Mass source term and periodicity

[mdegonville]

Hello everybody,

I'm using the brand new CS 3.0.0, with GUI, and I want to simulate a simple tube with periodicity between inlet and outlet to get a fully developed flow in the whole domain.
The problem is that the mass flow decreases until it reaches 0: so I wrote a subroutine using the ustsma.f90 examples, with a gamma calculated from the difference between the reference velocity (or a given mass flow) and the actual mass flow, inspired by answers that I found on this forum (I'm a beginner in FORTRAN 90).

But nothing happens, I mean, the file is compilated without error and the values are computed (I did some print within the file to check it out), but the mass flow keeps decreasing whatever the value of gamma.

Any idea about what is wrong in my code?

Thanks a lot,

Maximilien

[Yvan Fournier]

Hello,

I have'nt checked your ustsma.f90 file in detail, but mass source terms are useful for representation of mass sources or sinks only. So I suspect the flow simply "slows down" in your case (you can check looking at the maximum velcoity). For a periodic tube, you probably do not need to inject mass, but need to keep the flow moving, using a velocity source term, which means programmng ustsnv (or ustsns if using the ivelco = 0 keyword value).

Also, for developed flow, without need of periodicity, we have a new example of an inlet conditions, in the cs_user_boundary_conditions-auto_inlet_profile.f90 example (bascically, pseudo-periodicity is prescribed between the inlet and first cell; this works only if the first cell layer is orthogonal, and may require more iterations to set up than a purely periodic case, but may be simpler, and adaptable more easily to more complex geometries).

Regards,

Yvan

[mdegonville]

Thanks a lot ! It was actually quite simple to do it with the cs_user_boundary_conditions-auto_inlet_profile.f90 file...

But actually, if you could get a look to my file and explain me why it didn't modify anything, I would be thankful, as it could be useful in the future to do such modifications..

[mdegonville]

Hi again,

the auto-inlet subroutine was actually not what I searched for, as the cylinder was only a first step before simulating more complex profiled pipes. I tried to programm some loop in the cs_user_source_terms.f90, but I've got some problem with it:

I get the mass flow rate from the "inlet" (now an interior face thanks to periodicity) and the total surface using

Code: Select all

call getfac('inlet', ifac, lstfac)

do ielt = 1, ifac
  iifac = lstfac(ielt)
  iel = ifacel(1, iifac)
  flux = flux + propce(iel, ipproc(ifluma(iu)))
  surf = surf + surfan(iifac)
enddo

then I compute the bulk velocity from it using a value of propce(iel, ipproc(irom)) and the calculated surface, and then compute the MMT coefficient as:

mmt = rho * (uref - ubulk) / dtref

and finally apply the source term at each cell next to my inlet using

crvexp(1, iel) = volume(iel) * mmt


Then I got some problems:

- my density is absolutely not the one that I defined using the GUI (it is much smaller)
- the mass flow rate calculated is bigger than the one I expect
- the ifacel(1, iifac) works, but not the ifacel(2, iifac) that give me a 0 value for my velocity and my mass flow rate, which I don't understand as the inlet, as an interior face, must now have 2 neighboring cells with "normal" values
- finally, I don't even reach a stabilization of my bulk velocity (divergence or convergence to 0)

I also tried calculating the area-weighted average of the mass flow rate, and the velocity using rtp(iel, iu) (which also returns incoherent values with the results from the monitoring points), without success.

Can somebody help me?

Thanks a lot,

Maximilien.

[Yvan Fournier]

Hello,

Are you also using cs_user_physical_properties.f90 ? If you are, the explanation for the density should be there (it takes priority over the GUI). If not, you probably have a array overwrite somewhere in your code (which can be very bad)...

In any case, for mass flux, you should use the face-based values instead of cell values.

Here is an example for ustsnv with periodic conditions, which I'll let you adapt. Careful, I use nlfac where you use ifac, so as to keep ifac for the loop (instead of iifac). This example also only has a source term for U. You may need to adapt it for u,v,w depending on the orientation of your tube.

Code: Select all

! Local variables
integer ipcrom
integer ii, ielt, ifac, nlfac, cptfac, iflmas
double precision flowint, flowref, surf
double precision flown, flown1, dp

save flown, flown1, dp

!==================
! 1. Initialization
!==================

uref = 1.0d0 ! reference velocity: set your value here

if (ivar.ne.iu) return ! adapt if you need also iv and iw)

if (iwarni(ivar).ge.1) then
  write(nfecra,*) 'User source terms for variable ' , nomvar(ipprtp(ivar))
endif

ipcrom = ipproc(irom)
iflmas = ipprof(ifluma(iu))

!===============================
! 2. Source term for component u:
!===============================

flowref = uref*propce(1,ipcrom)
flowint = flown
flown = 0.0d0
surf  = 0.0d0
cptfac = 0

call getfac('inlet', nlfac, lstfac)

! First loop to determine surface and current flow

do ielt = 1, nlfac
  ifac = lstfac(ielt)
  surf = surf + surfan(ifac)
  flown = flown + propfa(ifac,iflmas)
  cptfac = cptfac + 1
enddo

if (irangp.ge.0) then
  call parsom(surf)
  call parsom(flown)
  call parcpt(cptfac)
endif

! convert flow to surface flow
! (if you need mass flow instead, for example with variable density, adapt it)

flown = abs(flown) / surf

if (ntcabs .eq. 1) then
  dp = 0.0d0
  flown  = flowref
  flown1 = flown
else
  flown1 = flowint
endif

! Recall force check

! Caution: this is valid only with a spatially constant time step
! (i.e. for an unsteady Code_Saturne calculation, not for a steady calculation)

dp = dp + (2.0d0*(flown-flowref)-(flown1-flowref))/(2.0d0*dt(1))

do ielt = 1, nlfac
  ifac = lstfac(ielt)
  do ii = 1, 2
    crvexp(1,ifacel(ii,ifac)) = - volume(iel) * dp
  enddo
enddo

write(nfecra,*) 'number of inlet faces = ', cptfac
write(nfecra,*) 'inlet surface = ', surf
write(nfecra,*) 'flow at time n at inlet = ', flown
write(nfecra,*) 'ref. flow', flowref
write(nfecra,*) 'flow at time n-1 at inlet = ', flown1
write(nfecra,*) 'flow forcing = ', dp

This is code adapted from a working example, but the adaptation is not tested, even for compilation, so I'll let you finish it.

Regards,

Yvan

[mdegonville]

Hello,

about the face-based mass flux, it was a mistake I corrected in the newer versions of my code.

Thank you a lot for your code, I didn't think of such an algorithm for the source term.

As I don't use the cs_user_physical_properties.f90 file (I'm using the GUI for most of the case configuration), I did some tests, and my propce(iel, ipcrom) actually returns... the time step used for the calculation of the cell.
Moreover, the dt() command doesn't work and returns me an error, not during the compilation but during the computation.

I tested my code on the 3.0.1 CS version on another computer, and the same thing appears. I join the code file to this post, hoping that someone could find the problem in the code (as I honestly don't know where such a mix between the density and the time step could occur in it...).

[Yvan Fournier]

Hello,

Did you use the GUI or cs_user_parameters.f90 to tell the code you would be using a variable density ? If not, density is assumed to be constant, so propce(iel, ipcrom) does not point to the correct array.

Regards,

Yvan

[mdegonville]

I work with a constant density, and the only subroutine is the cs_user_source_terms.f90. Is the density stored anywhere else in the case of a constant density?

[Yvan Fournier]

Hello,

Sorry, my mistake. The array for rho seems to always be present. I thought rho was optional, like the variable viscosity for scalars (ivisls), but it seems this is not the case.

I don't have much time to look at your user file this week, so I would suggest commenting most of that routine, adding your verifications (such as printing the value of rho or dt), and adding back (uncommenting) things progressively, to see where there might be a mistake. Assuming there is probably a "hidden" error in your code, once you find it, it should be easy to fix (bugs in the non-user parts of the code are not excluded, but rarer, so let's start with the more probable causes first).

Good luck !

Yvan

[mdegonville]

Hello,

I tested my code as you said, and it appears that the problem comes from "ipcrom = ipproc(irom)", as it returns "1" which doesn't seem to be the index of the density values at the cells.

All the code before is just declaring variables, initializing the reference velocity and stopping if the concerned variable isn't iu (I didn't modify much of this part in the code you gave to me, just adding my own variables and a counter).

I am wondering if, as the density is constant, it is not stored in the cell and boundary face arrays but just in the ro0 variable ? (and the irom values are just there in case of a variable density) Or is that what you meant in your previous post?

Thanks a lot!