Convergence issues while using temperature-dependent density and user drag function

I am trying to simulate the particle sedimentation in water using CFD-DEM approach. I have defined a user-defined function for density of water (varies with temperature as in Boussinesq approximation). If I use constant density of water (i.e. 1000 kg/m^3), then the simulation runs. However, if I use the user-defined density, I get an error of DT < DTmin (recovery not possible).
I have attached the files. Kindly help me resolve the issue.

Thank You

fluid_bed_dem_2d - water - ht - udf.zip (30.0 MB)

I’m not sure if you attached the right files, but the usr_drag.f in the uploaded ZIP file is invalid. I see an error message immediately when I start your simulation, due to the following code in usr_drag.f:

! The following error message is used to make sure that if a user
! defined drag law is invoked, that this routine has been modified.

!- REMOVE THE FOLLOWING ---------------------------------------------->>

      lDgA = 3*3.14159265359*Mug*DPM*VREL

      WRITE(ERR_MSG,9999)
      CALL LOG_ERROR()

 9999 FORMAT('ERROR 9999: The user-defined drag routine was invoked ', &
         'but this',/'generic error message exits. Either choose a ',  &
         'different drag law',/'or correct mfix/model/usr_drag.f')

!- END REMOVE --------------------------------------------------------<<

However, after changing the drag model to the default SYAM_OBRIEN I see a problem with your case. Setting fluid density to a constant value of 1000 kg/m³ the job runs to completion in a little over 13 minutes - there is a big dropoff in DT after about 0.3 seconds of simulation time - I don’t know exactly what is causing this:


Also noteworthy is that the gas temperature is between 0 and 1K in the VTK output (?)

However, when I enable the user-provided density model (usr_rog) and make the simplest possible UDF:

      SUBROUTINE USR_PROP_ROg(IJK)
      INTEGER, INTENT(IN) :: IJK
      RO_G(IJK) = 1000.0D0
      RETURN
      END SUBROUTINE USR_PROP_ROg

the simulation does not get past t=0, as can be seen in the MFiX status tab

{'crashed': False,
 'dt': 0.00014088441290426774,
 'elapsed_time': 46.77308797836304,
 'error': None,
 'finished': False,
 'io_time': 0.0032379627227783203,
 'nit': 14,
 'paused': False,
 'paused_time': 0.0,
 'pid': 11409,
 'remaining_time': 331949.44329041266,
 'residuals': [['HYDRO', 1.9945984832906412],
               ['THETA', 0.0],
               ['ENERGY', 1.5812548440622666e-16],
               ['SPECIES', 0.0],
               ['SCALAR', 0.0],
               ['K-EPS.', 0.0]],
 'run_name': 'fluid_bed_dem_2d',
 'running': False,
 'time': 0.0,
 'tstop': 1.0,
 'version': '95'}

The HYDRO residual is dominating.

Setting RO_G to 1000.0 in the UDF should be the same as using a constant density model, but it 's not. I do not know why this is. I’m working with the rest of the MFiX team to figure this out, I will update you when we have an answer.

Thank you for this bug report, and for your patience.

– Charles

The difference between constant density (ro_g0) and a UDF (usr_rog) that sets ro_g0 to 1000 may be that ro_g0 code treats the fluid as incompressible while the usr_rog code treats it as compressible. Stay tuned for further developments.