Internal coupling spurious velocity

[Luciano Garelli]

Hello,

I'm running in CS 6.0.4 a simple example using the internal coupling feature, but when setting different densities (rho_air=1, rho_solid=100) spurious velocities appears at the solid region as seen in the picture.
The densities and temperature are constant and the gravity is turned on. If gravity is turned off this behavior disappears.

I upload the test case. maybe something is not correctly set.

Thanks in advance.

Regards,
Luciano

[Yvan Fournier]

Hello Luciano,

Can you try activating the "fluid_solid" option for this ? :

After the call to cs_internal_coupling_add_volume (in csèuser_parameters.c), try

Code: Select all

  cs_stokes_model_t *stokes = cs_get_glob_stokes_model();
  stokes->fluid_solid = true;

Best regards,

Yvan

[Luciano Garelli]

Hello Yvan,

Thanks for your quick reply. I have tried your suggestion but I get the same behavior with some iterations in the velocity field in the solid part.

In the setup.log appears the text "Fluid-solid mode (disable dynamics in the solid part)"

Regards,
Luciano

setup.log

(24.02 KiB) Downloaded 305 times

run_solver.log

(269.24 KiB) Downloaded 307 times

[Yvan Fournier]

Hello Luciano,

OK, I'm not sure whether it is necessary to specify which cells are "solid" when using this, so I'll check with colleagues.

I'll keep you informed.

Best regards,

Yvan

[rodion]

Hello Luciano,

Another option to stop any momentum in the solid is to impose a high negative implicit source term for the momentum conservation equation as it is described in [Patankar, S., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation, Washington, D.C.]
The part of the book where this concept is given is in the attached figure.

Here is an axample of the code:

Code: Select all

call getcel('SOLID',nlelt,lstelt) 
do ilelt = 1, nlelt
  iel = lstelt(ilelt)
  crvimp(1, 1, iel) = -1.d+30
  crvimp(2, 2, iel) = -1.d+30 
  crvimp(3, 3, iel) = -1.d+30 
  crvexp(1, iel) = 0.d0
  crvexp(2, iel) = 0.d0
  crvexp(3, iel) = 0.d0  
enddo

You can embed this code in the subroutine ustsnv, just use the right name instead of the SOLID.

However, in CS6 the value -1E30 could result in some instabilities of the solver, giving you 10000 iterations in velocity from time to time. In that case just put something like -1E4, it should do the job as well.

Best regards,
Rodion

[Luciano Garelli]

Hello Yvan,

Thanks for the suggestion, I will try it. It is like to set a penalization for this term In the solid region.

Another workaround that I found is to set the same density in both, fluid and solid, and in order to maintain the volumetric heat capacity in the solid part I have set the heat capacity in the solid .

Regards,
Luciano

[Yvan Fournier]

Hello Luciano,

The previous suggestion was by Rodion (thanks for him).

It can be interesting as a workaround, as I guess it can be implemented as a user setup, but should not be necessary in an automatic setup, so we still need to look at that very soon.

Best regards,

Yvan

[Yvan Fournier]

Hello Luciano,

I do not reproduce any issues using the latest 6.0 branch (but I would expect 6.0.4 should be the same).

Just in case, I attach the updated cs_user_parameter.c file, as you were using a deprecated function (which should still be OK, but I updated it to be sure).

I ran on a debug build, which requires an update from master branch in src/base/post_util.f90 (to avoid division by 0 at the first time step when computing the boundary flux and Nusselt for postprocessing, which is a completely separate issue). L'll commit the patch tonight or tomorrow.

Best regards,

Yvan