[Rolanzo]

Hello everyone,

I am trying to solve system of non-linear equations for two phase flow (oil and water without capillary effect) for reservoir simulation application. Primitive variables are pressure (p) and water saturation (Sw). The equations are as follows

ddt(f(Sw),p) - laplacian(g(Sw),p) = A (Eqn.1)
ddt(k(Sw),p) + ddt(B,Sw) - laplacian(m(Sw),p) = C (Eqn.2)

where f, k, g, and m are functions of Sw and A, B, C are constant.

As far as solving system of equations in OpenFOAM, segregate method is the most common way, this can be done by solving each equation sequentially, which will automatically linearize both equations rather than solving both of them simutaneously. I have tried this method and the solution is not correct. So I am looking at the conventional way of solving this system of equations which is FIM (fully Implicit Method) since OpenFOAM-3.1 ext provides blockmatrix solver which allows users to solve coupled equations.

I have looked at example of this solver and it is quite different from my problem as the problem consists of two linear equations as below:

fvScalarMatrix TEqn
(
fvm : : div ( phi , T ) - fvm : : laplacian ( DT , T ) == alpha*Ts - fvm : : Sp ( alpha , T )
) ;

fvScalarMatrix TsEqn
(
- fvm : : laplacian ( DTs , Ts ) == alpha*T - fvm : : Sp ( alpha , Ts )
) ;

As I have gone thru the code, the key part is to discretize both equations then insert into blockMatrix as below

// Insert equations into blockMatrix
blockMatrixTools : : insertEquation ( 0 , TEqn , blockM , blockX , blockB ) ;
blockMatrixTools : : insertEquation ( 1 , TsEqn , blockM , blockX , blockB ) ;

Then you add the off-diagonal terms

forAll ( d , i )
d [ i ] ( 0 , 1 ) = alpha.value( )*mesh.V( )[ i ] ;
d [ i ] ( 1 , 0 ) = alpha.value( )*mesh.V( )[ i ] ;

The remove off-diagonal terms from RHS or source term

blockB [ i ] [ 0 ] = alpha.value ( )*blockX [ i ][ 1 ]*mesh.V( )[ i ] ;
blockB [ i ] [ 1 ] = alpha.value ( )*blockX [ i ][ 0 ]*mesh.V( )[ i ] ;

Since above system of equations are linear, the coupled part can be done easily. However, for my problem, I need to linearzed both equations first, so I use Newton Raphson method making both equations become

Let F = Eqn.1 and G = Eqn.2, R1(p,Sw) = RHS-LHS or residual of Eqn.1 , and R2 = RHS-LHS or residual of Eqn.2, then we have

Sp(dF(Sw)/dp,delp) + Sp(dF(p)/dSw,delSw) = R1(p,Sw) (Eqn.3)
Sp(dG(Sw)/dp,delp) + Sp(dG(p)/dSw,delSw) = R2(p,Sw) (Eqn.4)

where delp is p^n+1-p^n or difference between current and next iteration step and delSw is Sw^n+1-Sw^n. The initial guess for p^n+1 is 1.1*p^n and for Sw^n+1 is 1.1*Sw^n.

R1 and R2 are calculated using value of p and Sw from current iteration step. Same for, dF(Sw)/dp, dF(p)/dSw, dG(Sw)/dp, and dG(p)/dSw, all using value of p and Sw from current iteration step. Eqn.3 and 4 are to be solved using blockMatrix solver. Once p and Sw are updated, all other terms are recalculated and the iteration keeps going until delp and delSw become smaller than convergence value. The we move to next time step.

Now the problem is how to solve these coupled equations, this is what I have done:

fvScalarMatrix delpEqn
( fvm::Sp(dF(Sw)/dp,delp) == fvc(R1) );

fvScalarMatrix delSwEqn
( fvm::Sp(dG(p)/dSw,delSw) == fvc(R2) );

fvScalarMatrix SwinPEqn
( fvm::Sp(dF(p)/dSw,delSw) );

fvScalarMatrix PinSwEqn
( fvm::Sp(dG(Sw)/dp,delp) );

Now, assemble the blockMatrix
blockMatrixTools : : insertEquation ( 0 , delpEqn , blockM , blockX , blockB );
blockMatrixTools : : insertEquation ( 1 , delSwEqn , blockM , blockX , blockB);

As to my understanding, above commands put delpEqn and delSwEqn in diagonal part of the matrix and put fvc(R1) and fvc(R2) in block source; however, for off-diagonal terms (SwinPEqn and PinSwEqn), I am still not sure how to do this part, can I do something similar to example problem as below

forAll ( d , i )
d [ i ] ( 0 , 1 ) = alpha.value( )*mesh.V( )[ i ] ;
d [ i ] ( 1 , 0 ) = alpha.value( )*mesh.V( )[ i ] ;

Please advise? Also, is there a way to solve this coupled equations without using this Newton-Raphson method in OpenFOAM using fully implicit method? Any advices welcome!

Thank you!