Hi folks,

I am a bit new to CFD, and I am running into a problem that I wouldn't expect with a code I'm using. The code is 3D FV, using backward euler in time, power law for the advection-diffusion terms, and LBL TDMA as a solver.

I am simply trying to make it run a a few test cases on a uniform plug of chemical concentration in a square duct - observing the chemical's diffusion with no advecting velocity, or observing advection with diffusion turned off and an imposed uniform flow profile. In both cases, I am observing non-conservation. While LBL TDMA is converging at each time step, the total concentration begins to decrease slightly (at about 10^-5 per iteration).

I wouldn't expect to see this at all, since FV is a conservative method, my concentration plug is far from the inlet and outlet, and I have coded no flux BC's at the walls.

Two Questions -

1. What could be potential causes? Anything that is easily fixed?

2. If this non-conservation is unavoidable, and is simply a consequence of the discretization, are there tricks that can be played to minimize it while waiting for the code to march through time?

(3. As a followup regarding numerics in general - is it possible to write any code that is 100% conservative, or are they all like this, the only difference being that the non-conservation can be pushed to smaller and smaller magnitudes with improved schemes?)

Any insight is appreciated.

Thanks! Jon

A "thought problem" in 1 dimension is what's needed.

Consider volumes labeled V1, V2, ..., Vj, Vj+1, ..., Vn, stretching from one boundary to the opposite. The mass contained in V1 is D1 x V1, where D1 is the mass density in V1. F1 is the flux of D across the left face of V1 and F2 is flux of D across the right face of V1, etc.

For Vj, the balance equation is

Mass accumulation: Vj x [Dj(t+dt) - Dj(t)]/dt = (Fj x Aj) - (Fj+1 x Aj+1).

For Vj+1, MAC: Vj+1 x [Dj+1(t+dt) - Dj+1(t)]/dt = (Fj+1 x Aj+1) - (Fj+2 x Aj+2),

where Aj is the flow area of the left face of Vj.

F1 and Fn+1 are used to set boundary conditions.

Note the last term in the first equation and the next-to-last term in the 2nd equation are identical. In my experience, it's very easy to code these terms just slightly differently, making the mass leaving the left (jth) cell a little different from that entering the right (j+1st) cell. Of course, you need to do a similar hard-nosed comparison for the other two directions in your 3-d code.

Actually one trick (if you have sufficient memory available), is to calculate F1 x A1, F2 x A2, etc. in a setup loop, then go through the equations in a separate loop, just applying the pre-calculated fluxes. This just about halves the processor work to calculate the fluxes and guarantees that the flux entering the left side of the volume equals that leaving the right side of the left-hand neighbor.

Good luck!

Jim, thanks for your insight. I did investigate what you said, and unfortunately all of my fluxes are the same at the different interfaces. Following SV Patankar, to me this means that ae(i,j,k)=aw(i+1,j,k), for instance. And that's the case for all of them. I have to wonder what else could be wrong. I have a thought, and I'm wondering what you and others think -

For no velocity wall penetration, what is the proper way to implement a no flux condition on concentration? I simply solve the internal points (i.e. - x = (2:end-1)) in each direction, and then set the end points c(1)=c(2) and c(end)=c(end-1). However, I don't think it should matter anyway, because it can be shown (I think) that for no flux, aw(2) =0 and ae(end-1)=0. In other words, the points adjacent to the walls feel no influence from the wall, if there is no source generation term by the wall and we simply have no velocity penetration and no diffusive flux.

Thoughts are appreciated. Thanks again. -Jon

Jon,

There's another non-conservation issue that I overlooked.

If the mass conservation for each control volume is not exact, that imbalance introduces a small source or sink of the scalar in the scalar transport equation. For a Pantakar scheme, that means that the pressure equation needs to be pretty well converged at each time step in order to keep the scalar source-sink small. To test this, try tightening the convergence on the pressure equation by a factor of 10 and observe what changes (if anything) in your scalar distribution. This can get a bit difficult if your code is writen in real*4 (single precision) and you're already grinding things down to roundoff.

If your incompressible scalar equation is

c_t + (u x c)_x + (v x c)_y = RHS,

and you expand that to

c_t + c x (u_x + v_y) + u x c_x + v x c_y = RHS,

the term in ( ...) is the residual on the continuity equation. If it's not numerically zero, that term is a source or sink of the scalar. The source/sink is strictly the result of insufficient convergence of the pressure equation - resulting in violation of mass conservation.

Of course, if you already have good convergence of the pressure equation, I've wasted your time : )