Does subroutine overwrite correctly?

[Tsubasa]

Hello,
Now I am trying 2D simulation of a building's room by k-omega SST.
Here, I am doing pressure-driven flow. (I imposed Dirichlet pressure value at inlet and outlet.) When I assigned 10m/s as inlet velocity in GUI, 4m as hydraulic diameter in the GUI, and pressure by subroutine:

Code: Select all

call getfbr('window1_extruded', nlelt, lstelt)
do ilelt = 1, nlelt
   ifac = lstelt(ilelt)
   itypfb(ifac) = ientre

!Dirichlet B.C for pressure:
   icodcl(ifac,ipr) = 1
   rcodcl(ifac,ipr,1) = 101325.633

!Neumann B.C for velocity:
   icodcl(ifac,iu)= 3
   icodcl(ifac,iv)= 3
   icodcl(ifac,iw)= 3
enddo

I supposed that "Neumann B.C for velocity" in subroutine overwrite "10m velocity" in GUI. Even though I imposed 10m inlet velocity, velocity magnitude at inlet is less than 2m or 1m when I visualized result. Therefore, we do not have to care about velocity in GUI. (E.g. we can assign every velocity magnitude like 0m/s 100m/s, 1000m/s and so on.)

However, when I changed 4m/s to 1m/s as a velocity in the GUI, this solution became divergence even though I used completely same files and computational condition.
If my subroutine overwrite correctly, velocity in GUI never affects its simulation result.
However, this cannot be true.

For another reason for this, is the possibility that my subroutine code is wrong?

Let me post succeeded case (velocity:10m/s, hydraulic diameter:4m) and divergent case (velocity:1m/s, hydraulic diameter: 4m) and could someone test this simulation?

Succeeded case:

succeeded_case.zip

(4.68 MiB) Downloaded 117 times

Divergent case:

divergence_case.zip

(4.68 MiB) Downloaded 116 times

The difference between them is only inlet velicuty (10m/s or 1m/s) in GUI.

Best regards,
Tsubasa

[Yvan Fournier]

Hello,

The "rcodcl" values were set by the GUI for velocity (iu, iv, iw); so setting them to the actual Neumann values (rcodcl(ifac, iu, 3) = 0 for example) might make a difference. Id would not expect the Dirichlet value (rcodcl(ifac, iu, 3)) to be used with a Neumann condition, but just in case, you can set it to 0.

Also, regarding the stability of the computation (from you other post on this subject), purely pressure-driven flows might be a bit more unstable, so several factors may help:
- extrude the boundary over 1 cell layer at open windows (for a better mesh quality in that sensistive region)
- use a lower time step (check the Courant number), or switch to a "local in space" time step.

Running your case with the 7.0-beta version, without extrusion, I reproduced your crash with the constant time step you chose, but ran fine with local time step.

Best regards,

Yvan

[Tsubasa]

Hello Yvan,

Thank you for your help.

As you said, simulation did not become divergence if I used "Time varying(adaptive)" of time setting.
However looking at the result, flow is reversed (from outlet to inlet.)
Could you check whether reverse flow happened or not.

Or could you post the setup case without result (RESU), which was done by varying time step?

Best regards.
Hamada

[Tsubasa]

Hello,

Even though I used this subroutine,

------------------------------------------------------------------
call getfbr('window1_extruded', nlelt, lstelt)
do ilelt = 1, nlelt
ifac = lstelt(ilelt)
itypfb(ifac) = ientre

!Dirichlet B.C for pressure:
icodcl(ifac,ipr) = 1
rcodcl(ifac,ipr,1) = 101325.633

!Neumann B.C for velocity:
icodcl(ifac,iu)= 3
rcodcl(ifac, iu, 3) = 0
icodcl(ifac,iv)= 3
rcodcl(ifac, iv, 3) = 0
icodcl(ifac,iw)= 3
rcodcl(ifac, iw, 3) = 0
enddo
---------------------------------------------------------------------

simulation became divergent if I set a velocity value as 1m in GUI.
But If I set it as 10m, simulation is OK.

I am wondering why this subroutine cannot overwrite the velocity of GUI.

Best regards,
Hamada

[Yvan Fournier]

Hello,

The velocity is icodcl(iac, ivar, 1) for Dirichlet (set by the GUI), and icodcl(iac, ivar, 3) for Neumann in you case, so you do not overwrite the Dirichlet value if you do not set it.

I am not sure where it is used (perhaps mixed conditions), but you could overwrite it more completely.

Regards,

Yvan

[Tsubasa]

Hello Yvan and everyone,

Yvan, do you mean velocity set by GUI is reflected in the simulation result even if I used the above subrouine code?
However, when I tested purely pressure-driven-2D- Hagen–Poiseuille flow with same subroutine code above, the reuslt is almost same to analytical value even if I set 10m/s, 100m/s and 1000m/s for velocity.

Also, whatever I set velocity magnitude in GUI, its velocity cannot appear in postprocessing. For example, when I set 100m/s for velocity for 2D- Hagen–Poiseuille flow by GUI, 100m/s never happened. This is subroutine code:

cs_user_boundary_conditions.f90

(21.49 KiB) Downloaded 119 times

This figure shows comparison between CFD and analytical value. In GUI, I set 10000m/s, but it is not reflected in simulation.

Therefore, I think I do not have to set velocity in cubroutine code if I used above subroutine code for pure pressure-driven flow. However, is it wrong? Do I have to set something by subroutine code to remove velocity component compeletely?

Could someone tell me where is there example or reference to do it in Doxygen, or could someone give me an example to overwrite it correctly?

Best regards,
Tsubasa

[Yvan Fournier]

Hello,

As I said in my previous posts, I do not know exactly how the previous GUI setting manages to have an influence on the computation. I suspect if is used to compute some reference value bat am not sure

I did not have the time to check that in detail, so making sure you overwrite the theoretically unused values is a first step.

Best regards,

Yvan

[Tsubasa]

Hello,

So far I did 2 dimentional simulation with 2D-mesh, which has single layer in Z direction.
As a cause of a lot of problem, I suspected that this problem is cased by 2D-mesh.
Then, I used a new mesh, which is same geometry but has 2 layer in Z direction.
After that, I was suprised that simulations became stable, and I can set any hydraulic diameter and velocity magnitude.

Yvan, are there some limitations for 2 dimentional simulation in Code_saturne?
I used k-omega SST model.
Do the developer team assume that simulation should be done by 3 dimentional mesh for k-omega SST or something?

Best regards,
Tsubasa

[Yvan Fournier]

Hello,

code_saturne is usually more robust/stable using "orthogonal" or nearly "orthogonal" meshes, especially at outlets, buf if possible at all boundaries, and even better, in the domain. This means the vectors joining cell centers and face centers are as well aligned as possible with face normals.

Using a tetrahedral mesh, mesh quality in this regard is lower. If you use a single cell layer (2D mes), it is much better if this is an extruded 2D mesh (with prisms) than if you use tetrahedra (which will have non orthogonalty at boundaries). In this case, a single layer (with symmetry boundary conditions on top and bottom layers) should be fine.

I do not remember how your mesh is in this aspect.,

Regards,

Yvan

[Tsubasa]

Hello Yvan,

I have a good news.
So far I used 2D simulation with prism mesh, but its simulation did not work.
However, after I changed mesh from prism to hexa, everything is fine.



Also, I used three types of mesh; Tetra, Prism and Hexa, for pressure-driven flow,
and let me share my experiences.

Tetra; calculation never started. Error occurred at first step.
Prism; calculation start sometimes, and it depends on subroutine and GUI setting.
Hexa; calculation is very stable and fast even by any subroutine and GUI setting.

Best regards,
Tsubasa