## Configuration Files

### couplingProperties (CFD/constant/)

```cpp
modelType B;

couplingInterval 40;  // CFD_dt / DEM_dt = 0.001 / 2.5e-5 = 40

voidFractionModel divided;
locateModel engine;
meshMotionModel noMeshMotion;
IOModel basicIO;
dataExchangeModel twoWayMPI;
averagingModel dense;
smoothingModel constDiffSmoothing;

forceModels
(
    KochHillDrag
    LaEuScalarTemp
    Archimedes
);

momCoupleModels
(
    implicitCouple
);

constDiffSmoothingProps
{
    lowerLimit 0.0;
    upperLimit 1e99;
    smoothingLength 0.012;  // 4× particle diameter
}

implicitCoupleProps
{
    velFieldName           "U";
    granVelFieldName       "Us";
    voidfractionFieldName  "voidfraction";
    KslLimit               30;
}

KochHillDragProps
{
    velFieldName           "U";
    densityFieldName       "rho";
    voidfractionFieldName  "voidfraction";
}

dividedProps
{
    alphaMin   0.3;
    volScale (1.0);
}
```

### controlDict (CFD/system/)

```cpp
application     cfdemSolverPisoScalar;
deltaT          0.001;
writeControl    runTime;
writeInterval   0.05;  // Problem occurs with 0.05, NOT with 1.0
```

### in_run.Drum (DEM script) - Key sections

```liggghts
units           si
atom_style      granular
boundary        f f f
timestep        2.5e-5

# Rotation
variable rpm equal 2.0
variable D_rot equal 60/${rpm}  // 30 s period

fix movecad1 all move/mesh mesh cad1 rotate origin 0.0 0.0 0.0 axis 0.0 1.0 0.0 period ${D_rot}
fix movecad2 all move/mesh mesh cad2 rotate origin 0.0 0.0 0.0 axis 0.0 1.0 0.0 period ${D_rot}
fix movecad3 all move/mesh mesh cad3 rotate origin 0.0 0.0 0.0 axis 0.0 1.0 0.0 period ${D_rot}

# CFDEM coupling
fix cfd all couple/cfd couple_every 10 mpi
fix cfd2 all couple/cfd/force/implicit
fix cfd3 all couple/cfd/convection T0 ${start_temp}

# Output - Problem occurs with dumptime = 0.05
variable dumptime equal 0.05
variable dumpsteps equal ${dumptime}/${dt}
dump dmp all custom/vtk ${dumpsteps} post/particle_*.vtk id type x y z vx vy vz fx fy fz omegax omegay omegaz radius mass f_Temp[0] f_heatFlux[0]
```

---

## What I Have Already Tried

### 1. Drag Model Variations
| Model | Result |
|-------|--------|
| GidaspowDrag | Unstable (particles ejected) |
| **KochHillDrag** | More stable, but still ejections |
| DiFeliceDrag | Not yet tested |

### 2. Coupling Parameters
| Parameter | Tested Values | Result |
|-----------|---------------|--------|
| KslLimit | 100 → 50 → 30 | Reduced frequency, not eliminated |
| alphaMin | 0.1 → 0.3 → 0.4 | Minor improvement |
| smoothingLength | 0.009 → 0.012 → 0.018 | Minor improvement |
| couplingInterval | 40 (correct) | Verified calculation |

### 3. Solver Settings
- Added solvers for sSmoothField and vSmoothField (PCG/DIC)
- Verified all residuals are converging properly

### 4. Output Intervals (KEY FINDING)
| writeInterval / dumptime | Result |
|--------------------------|--------|
| **0.05 s** | **Instability occurs** |
| **1.0 s** | **Stable simulation!** |


## Questions

1. **Is there a known issue with I/O synchronization between CFD writeInterval and DEM dump when using frequent output intervals?**

2. **Could there be a buffer overflow or memory issue when writing many output files in quick succession?**

3. **Is there a recommended ratio between writeInterval/dumptime and couplingInterval?**

4. **Are there any special considerations for rotating geometries with implicit coupling?**

5. **Should I use explicitCouple instead of implicitCouple for rotating systems?**

6. **Is snappyHexMesh recommended over blockMesh for this type of geometry?**