CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Solving non-linear coupled equations using blockmatrix solver (OpenFOAM-3.1ext)

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 2 Post By Rolanzo

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 12, 2018, 21:49
Default Solving non-linear coupled equations using blockmatrix solver (OpenFOAM-3.1ext)
  #1
New Member
 
Rolanzo
Join Date: Jan 2016
Location: Texas, USA
Posts: 18
Rep Power: 10
Rolanzo is on a distinguished road
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!
Kummi and jennyjian like this.
Rolanzo is offline   Reply With Quote

Reply

Tags
---

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
conjugate heat transfer in OpenFOAM skuznet OpenFOAM Running, Solving & CFD 99 March 16, 2017 06:07
simpleFoam error - "Floating point exception" mbcx4jc2 OpenFOAM Running, Solving & CFD 12 August 4, 2015 03:20
Compressor Simulation using rhoPimpleDyMFoam Jetfire OpenFOAM Running, Solving & CFD 107 December 9, 2014 14:38
SLTS+rhoPisoFoam: what is rDeltaT??? nileshjrane OpenFOAM Running, Solving & CFD 4 February 25, 2013 05:13
Error while running rhoPisoFoam.. nileshjrane OpenFOAM Running, Solving & CFD 8 August 26, 2010 13:50


All times are GMT -4. The time now is 15:27.