|
[Sponsors] |
Solving non-linear coupled equations using blockmatrix solver (OpenFOAM-3.1ext) |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
April 12, 2018, 21:49 |
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 |
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! |
|
Tags |
--- |
Thread Tools | Search this Thread |
Display Modes | |
|
|
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 |