code_saturne User's Forum

Hello,
i am preparing an easy test file for analyze the usage of mass source terms (with only one processor) in code saturne.
To get a thermal balance i am using the routine "usproj.f90" which i have bounded to only balance at inlet and outlet.
For the mass balance i have integrated a function like m = A*v*rho
 
For the inlet:
massin = 0.0d0 CALL GETFBR('in-hot or in-cold',NLELT,LSTELT) !========== do ilelt = 1, nlelt ifac = lstelt(ilelt) ! --- Element de bord iel = ifabor(ifac) betrag = sqrt(surfbo(1,ifac)**2 + surfbo(2,ifac)**2 + surfbo(3,ifac)**2 ) normal(1) = -surfbo(1,ifac)/betrag normal(2) = -surfbo(2,ifac)/betrag normal(3) = -surfbo(3,ifac)/betrag massin = massin + propce(iel,ipproc(irom(iphas)))*betrag*normal(1)*rtp(iel,iu(iphas)) & + propce(iel,ipproc(irom(iphas)))*betrag*normal(2)*rtp(iel,iv(iphas)) & + propce(iel,ipproc(irom(iphas)))*betrag*normal(3)*rtp(iel,iw(iphas))
 
For the outlet:
massout = 0.0d0 CALL GETFBR('out-cold or out-hot',NLELT,LSTELT) !========== do ilelt = 1, nlelt ifac = lstelt(ilelt) ! --- Element de bord iel = ifabor(ifac) betrag = sqrt(surfbo(1,ifac)**2 + surfbo(2,ifac)**2 + surfbo(3,ifac)**2 ) normal(1) = surfbo(1,ifac)/betrag normal(2) = surfbo(2,ifac)/betrag normal(3) = surfbo(3,ifac)/betrag massout = massout + propce(iel,ipproc(irom(iphas)))*betrag*normal(1)*rtp(iel,iu(iphas)) & + propce(iel,ipproc(irom(iphas)))*betrag*normal(2)*rtp(iel,iv(iphas)) & + propce(iel,ipproc(irom(iphas)))*betrag*normal(3)*rtp(iel,iw(iphas))
 
For the balance:
massbal = massin - massout
 
While for the energy balance for timesteps < 100 i get an average about zero, the average for the mass balance is about 2.4e-6.
In ensight 9.1.1 i did made an clip at the inlet and one at the outlet to calculate the mass balance and i got a average about 2e-6, too.
In the listing of the calculation the inlet and outlet have the same mass flow ( 0.760000e-2 ).
I do not use mass source terms yet and my time steps are uniform and constant. My grid has an constant hexa cell size of 0.001m and the maximum velocity of 50 m/s. The time step is 0.0002 s. That is coarse, but it should work and the listing shows the right mass balance.
The calculation is incompressible.
 
Does anyone know what is wrong with my mass balance in "usproj.f90" and ensight 9.1.1 ?
I would be glad if i could fix it.
 
Best regards,
Patrick

I am wondering if is not just a round-off error?
Perhaps you should try to use directly the mass flux stored in propfb(ifabor, ipprob(ifluma(iu(1)))) and so on.

Don't you think 2.6e-6 is a bit big for a round-off error?
And should a round-off error not have a average around zero?
I did tested the the mass balance with a flow through a cube and the error is bit greater (something about 3e-6).
Strange is, that the error is still the same in every time step. It only varies at the 8th or 9th decimal place.
 
If i only make a balance of the velocity (the flow is incompressible and the inlet and outlet has the same space) i see, that the velocity has the same error. It is a bit taller, but the same error.
 
Best regards,
Patrick

Hello,
As mass fluxes are defined at faces and not cell centers, computing a correct balance requires using the values in propfb(ifabor, ipprob(ifluma(iu(1)))) and not directly cell values (otherwise, you add an interpolation error).
The balance that is output in the listing file is generated in src/base/typecl.f90, so you might want to look at that file to see how it is done: basically, we have something like this:
do ii = 1, ntypmx
flumty(ii) = 0.d0
enddo

iuiph = iu(iphas)
iflmab = ipprob(ifluma(iuiph))

do ii = 1, ntypmx
istart = idebty(ii,iphas)
iend = ifinty(ii,iphas)
do jj = istart, iend
ifac = itrifb(jj,iphas)
flumty(ii) = flumty(ii) + propfb(ifac,iflmab)
enddo
if (irangp.ge.0) then
call parsom (flumty(ii))
endif
enddo

where ii loops on boundary condition types (the number matching each type is defined in include/base/paramx.h).
You may adapt this calculation to usproj.f90 if you want to generate a plottable file instead of just a printout in the  "listing" file.
Best regards,   Yvan

In addition to what Yvan wrote, I would say that even if you add no interpolation error when computing a "mass flux" at cell centers, you still wouldn't have the correct mass flux. The mass flux that respects the mass conservation [div(rho.u) = 0] is the mass flux defined at faces, that's why you need to use the code given in example in the previous post, or something like:
massin = 0.0d0
iflmab = ipprob(ifluma(iu(1)))

call getfbr('in-hot or in-cold', nlelt, lstelt)

do ilelt = 1, nlelt
  ifac = lstelt(ilelt)
massin = massin + propfb(ifac, iflmab)
enddo

if (irangp.ge.0) then
call parsom(massin)
endif

David

Thak you very much,
but now i would like to calculate the mass balance from the components velocity, surface area and density.
I think the velocity is stored in the cell centers, or?
 
best regards,
Patrick

Hello,
Yes, velocity is stored at cell centers, but if you calculate the mass balance from the cell-centered velocity, surface area and density, you will have a small consistency error as before. The mass fluxes stored in propfb are already based on the velocity, surface area, and density, so I see no reason not to use them.
Note that mass fluxes at boundary faces are in propfb, while mass fluxes at interior faces are in propfa (if you are using source terms and want to obtain a balance in a restricted part of the domain, you may need both; if you only need the balance for the full domain, boundary faces are enough).
Best regards,
  Yvan

Furthermore, the mass flux stored at face centers is the right mass flux. You won't have something correct at cell centers for the following reason, the Navier-Stokes resolution is done thanks to a segregated approach:
First, compute a predicted velocity u * , solution of the following equation (where q n is the mass flux at the previous time-step):
(u * - u n ) / dt + div[ q n .u * - nu.grad(u * ) ] = -grad(p n ) + div[ nu. t grad(u n ) ] + ...

Then, solve the pressure equation and get the correct mass flux (at faces) that balance the mass conservation equation:
div[ dt.grad(p) ] = div[ rho.u * ]
And the final velocity field is:
u (n+1) = u * - dt/rho . grad(p)

But, to be complete in my description of Code_Saturne algorithm, this is a bit more complex... In order to avoid potential decoupling in the calculation of the pressure gradient (due to the co-located formulation), we add a Rhie & Chow filter and actually solve:
div[ dt.grad(p) ] = div[ rho.u * + dt.grad cells (p n ) ] - div[ dt.grad faces (p n ) ]

To sum up, the only valid mass flux is the face-stored mass flux. Otherwise, you won't really account for everything.
Hope that clarifies a bit (and not the contrary...)
David
ps: I also hope not to have miss something in my equations (not easy to write it in the forum...) but if you see an error, let me know!

Thank you David and Yvan!
Looks like a kind of SIMPLE pressure correction procedure for me. The correct velocity is stored in the inner faces.  The filter is to much for me ;)
Sorry if that question is bit silly, but in the result files the velocity  and all other scalars are only available in the cell centers. Is that velocity the correct one interpolate to the cell centers?
I did try to calculate the mass flux in a incompressible flow by EnSight, but even the area mean velocity at the inlet and outlet have had an error like the mass flux.