Dear All, I am solving 1D multiphase flow in a pipe. The system of equations looks as follows:

d rho_i/dt + d rho_i V_j/dx =G_i (i=1..4, j=1,2) – discontinuity equation

0 = d p(rho_i)/dx +d (mu(rho_i) d V1/dx)/dx +Fmas – momentum for main phase, inertia is small.

0= d p(rho_i)/dx –K(V2-V1) - momentum for dispersed phase.

Boundary conditions: rho_i (0) are given, V1(0)=V2(0), -p(L)+2*(mu*dV1/dx) (L)=-pa.

Without viscous terms the system could be discretized and linearized on 2 points. This brings to equation:

A_i*deltaU_i –B_i*deltaU_i-1=F_i (1)

Also substituting deltaU_i with deltaU_i-1 repetitively one can get

deltaU_n=Phi*deltaU_0+theta.

These condition multiplied by discretized BC gives additional condition at x=0. Then using (1) U_i could be found. Without viscous terms (both in equations and BC) the method I use gives good convergence. If I add viscous term to BC convergence slows down dramatically and sometimes it diverges. In the method I use it's not possible to include viscous terms in the equations.

The questions is how to include viscous term (implicitly!) and how to force the convergence.

Best regards Oleg

Have you tryed using a deferred correction for your viscous terms. (or part of your viscous term to keep your stencil compact.)

Dear John, I am not strong in CFD terminology. Do you mean taking viscous term from previous time-step or previous iteration. Yes I have done both. It leads to divergence in the computation. The problem also is that viscous term is small in most of the flow and becomes dominant near the boundary x=L. This makes(?) possible to include it only into a boundary conditions. But having derivative there slows convergence dramatically (from 2 to 30 iterations). Could it be so because equations are not properly linked with boundary conditions?

Thank you, Oleg

What is the solver you are using??

Dear John,

I have briefly discussed the method that I use to solve the problem in my first message. I have an equation in the form: G(U_i,U_i-1)=0 without viscous terms. After linearisation it comes to

A_i*deltaU_i –B_i*deltaU_i-1=F_i or deltaU_i = P_i* deltaU_i-1-Q_i (P=A^-1B,Q=A^-1F) (1)

Removing intermediate terms I can calculate Phi and theta so that:

deltaU_n=Phi*deltaU_0+theta. (2)

I also have a boundary condition which after linearisation comes to

alpha* deltaU_n+beta=0. (3)

in both cases – pressure or stress is equal to atmospheric. (For the last case I also use (1) to calculate deltaU_n-1 as a function of delta_U_n). I multiply then (2) by (3) and get additional relation for deltaU_0. Then I solve (1) with all parameters known at the left boundary.

That's the method I use. It converges very fast in case of given pressure and slow in case of stress biundary condition.

Oleg

If your using an implicit method then you give the final linear equation to a linear solver to be solved. What is this linear solver??? CGSTAB ?? SIP ??? SOR??? If your solution is diverging then it may be a good idea to check whether your matrix is positive definite and diagonally dominate, but these types of things really depend on the solver you are using.

1) Check what solver you are using 2) find out what conditions this solver requires

for convergence 3) check to see if these conditions are met.

Dear John,

I just do exactly what I describe in previous message – no special solver. Just compute matrixes, inverse it and multiply as described. Possibly it's better to use standard solver for linear equations but I do it strait forward. I'll check diagonal dominance and positive definition of matrixes.

Thank you for comments

Oleg