
[Sponsors] 
May 24, 2011, 15:16 
Solving Coupled Set of PDEs

#1 
New Member
Nick Mai
Join Date: Feb 2011
Posts: 5
Rep Power: 6 
Hi Everyone,
I'm trying to solve the following coupled system with OpenFOAM: A) B) I've realized that you can use matrices to represent the switch and declare a tensor as follows: Code:
tensor j(0,1,0,1,0,0,0,0,0); solve { fvm::ddt(U) == j & fvm::laplacian(U) } Thanks for reading! EDIT: I've found coupledFvScalarMatrix in the OpenFOAM 1.6ext version. However, from the examples that I see, it seems as though it is intended for multiple domains/meshes. Is it possible to use the coupledFvScalarMatrix for a single mesh? If so, then my problem is solved. Nick 

May 24, 2011, 16:05 
try blockCoupledScalarTransportFoam

#2 
Senior Member
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 543
Rep Power: 18 
What about the blockCoupledScalarTransportFoam? I know that is on the same mesh, with two equations coupled to eachother through a source term
Code:
fvScalarMatrix TEqn ( fvm::div(phi, T)  fvm::laplacian(DT, T) == alpha*Ts  fvm::Sp(alpha, T) ); TEqn.relax(); fvScalarMatrix TsEqn (  fvm::laplacian(DTs, Ts) == alpha*T  fvm::Sp(alpha, Ts) ); Dan 

May 25, 2011, 02:15 

#3 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26 
Hi,
there is a coupled solver under development in extend, as pointed out in comment 2 (warning: it is a bit complicated at a first glance). However, you could define j as a volTensorField, and initialise it with your value of j. Code:
volTensorField j ( IOobject ( "j", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionedTensor("jTensor", dimless, Foam::Tensor(0,1,0,1,0,0,0,0,0)) ); dimensionSet(.....) if the tensor has dimensional units. Keep in mind however that the solution won't be coupled but sequential: meaning component by component. Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image. OpenQBMM  An opensource implementation of quadraturebased moment methods Last edited by alberto; May 25, 2011 at 02:16. Reason: Added comment 

November 24, 2011, 09:09 
blockCoupledScalarTransportFoam  coupling via operator

#4 
New Member
Christian Jungreuthmayer
Join Date: Mar 2009
Posts: 9
Rep Power: 8 
Hi!
I am trying to implement a solver for the following two coupled equations: Equation 1: laplacian(D1, phi_real)  laplacian(D2, phi_img) = 0 Equation 2: laplacian(D1, phi_img) + laplacian(D2, phi_real) = 0 First, I implemented a segregated version: solve( fvm::laplacian(D1, phi_real)  fvc::laplacian(D2, phi_img)); solve( fvm::laplacian(D1, phi_img) + fvc::laplacian(D2, phi_real)); which seems to work for weakly couplings. Naturally, this approach is rather useless if the coupling grows stronger. Hence, I had a look at 1.6.ext's solver "blockCoupledScalarTransportFoam". I guess I understand the main idea of it. I think one of the critical parts is the manipulation of the offdiagonal elements and the block source in order to create the correctly coupled system which looks like that in blockCoupledScalarTransportFoam.C: forAll(d, i) { d[i](0,1) = alpha.value()*mesh.V()[i]; d[i](1,0) = alpha.value()*mesh.V()[i]; blockB[i][0] = alpha.value()*blockX[i][1]*mesh.V()[i]; blockB[i][1] = alpha.value()*blockX[i][0]*mesh.V()[i]; } The above approach is rather simple for a coupling of the equations via the scalar "alpha". However, can something similar also be done if the coupling is done via an operator such as laplacian()? I tried to figure that out, but I am stuck. Any helpful comments would be very appreciated. Cheers, Christian 

July 26, 2012, 10:41 

#5 
New Member
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5 
hi,
is there anybody still following this topic? i'm having exactly the same interests as christian has pointed out. The blockCoupledScalarTransportFoam is very understandable for the idea of block matrix construction. However, when it comes to the coupling of transport terms (like grad(p) in the momentum equation, not as simple as the source term coupling in the blockCoupledScalarTransportFoam case), it is hard to construct the block because no implict gradient scheme is available does anyone have idea about a fvm::grad() scheme? 

July 27, 2012, 01:20 

#6 
New Member
Christian Jungreuthmayer
Join Date: Mar 2009
Posts: 9
Rep Power: 8 
Hi Tian,
Yes, I do. At least I did. I managed to solve my problem. If my memory servers me right then I wrote a solver that could solve the following coupled equations grad * (sigma * grad(phi_real))  grad * ((E0*Er*w) * grad(phi_imag)) = 0 grad * (sigma * grad(phi_imag)) + grad * ((E0*Er*w) * grad(phi_real)) = 0 using the following code: #include "fvCFD.H" #include "fieldTypes.H" #include "Time.H" #include "fvMesh.H" #include "blockLduSolvers.H" #include "VectorNFieldTypes.H" #include "volVectorNFields.H" #include "blockVectorNMatrices.H" #include "blockMatrixTools.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" # include "createFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // # include "CourantNo.H" for (runTime++; !runTime.end(); runTime++) { Info<< "Time = " << runTime.timeName() << nl << endl; # include "readSIMPLEControls.H" for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix phi_realEqn (  fvm::laplacian(EpsilonTimesAngVel, phi_real) ); phi_realEqn.relax(); fvScalarMatrix phi_realEqn2 (  fvm::laplacian(Sigma, phi_real) ); phi_realEqn2.relax(); fvScalarMatrix phi_imgEqn (  fvm::laplacian(EpsilonTimesAngVel, phi_img) ); phi_imgEqn.relax(); fvScalarMatrix phi_imgEqn2 (  fvm::laplacian(Sigma, phi_img) ); phi_imgEqn2.relax(); // Prepare block system BlockLduMatrix<vector2> blockM(mesh); BlockLduMatrix<vector2> blockM2(mesh); // Grab block diagonal and set it to zero Field<tensor2>& d = blockM.diag().asSquare(); d = tensor2::zero; Field<tensor2>& d2 = blockM2.diag().asSquare(); d2 = tensor2::zero; // Grab linear offdiagonal and set it to zero Field<tensor2>& u = blockM.upper().asSquare(); Field<tensor2>& l = blockM.lower().asSquare(); u = tensor2::zero; l = tensor2::zero; Field<tensor2>& u2 = blockM2.upper().asSquare(); Field<tensor2>& l2 = blockM2.lower().asSquare(); u2 = tensor2::zero; l2 = tensor2::zero; // forAll(u, j) // { // u[j](0,0) = 0.0; // u[j](1,0) = 0.0; // u[j](0,1) = 0.0; // u[j](1,1) = 0.0; // u2[j](0,0) = 0.0; // u2[j](1,0) = 0.0; // u2[j](0,1) = 0.0; // u2[j](1,1) = 0.0; // } // forAll(l, j) // { // l[j](0,0) = 0.0; // l[j](1,0) = 0.0; // l[j](0,1) = 0.0; // l[j](1,1) = 0.0; // l2[j](0,0) = 0.0; // l2[j](1,0) = 0.0; // l2[j](0,1) = 0.0; // l2[j](1,1) = 0.0; // } vector2Field& blockX = blockT.internalField(); vector2Field& blockX2 = blockT2.internalField(); vector2Field blockB(mesh.nCells(), vector2::zero); vector2Field blockB2(mesh.nCells(), vector2::zero); // Inset equations into block Matrix blockMatrixTools::insertEquation(0, phi_realEqn, blockM, blockX, blockB); blockMatrixTools::insertEquation(1, phi_imgEqn, blockM, blockX, blockB); blockMatrixTools::insertEquation(0, phi_realEqn2, blockM2, blockX2, blockB2); blockMatrixTools::insertEquation(1, phi_imgEqn2, blockM2, blockX2, blockB2); forAll(d, i) { d[i](0,1) = d2[i](1,1); d[i](1,0) = d2[i](0,0); blockB[i][0] += blockB2[i][1]; blockB[i][1] = blockB2[i][0]; } forAll(u, j ) { u[j](0,1) = u2[j](1,1); u[j](1,0) = u2[j](0,0); } forAll(l, j ) { l[j](0,1) = l2[j](1,1); l[j](1,0) = l2[j](0,0); } // Block coupled solver call BlockSolverPerformance<vector2> solverPerf = BlockLduSolver<vector2>::New ( word("blockVar"), blockM, mesh.solver("blockVar") )>solve(blockX, blockB); solverPerf.print(); // Retrieve solution blockMatrixTools::blockRetrieve(0, phi_real.internalField(), blockX); blockMatrixTools::blockRetrieve(1, phi_img.internalField(), blockX); phi_real.correctBoundaryConditions(); phi_img.correctBoundaryConditions(); } volVectorField gradPhi_real = fvc::grad(phi_real); volVectorField gradPhi_img = fvc::grad(phi_img); # include "write.H" } Info<< "End\n" << endl; return(0); } I do not remember any details anymore. However, I remember that the solver worked nicely. I hope this helps, Christian 

July 27, 2012, 02:39 

#7 
New Member
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5 
hi, christian,
thanks soooo much for your kind reply! i understand your code generally, but i just have a question, in my equation, the coupling term is the gradient of p (in your equation, the coupling is the laplacian, and luckily OpenFOAM has the fvm::laplacian(), but they don't have fvm::grad()...that's my big problem!). do you also have any idea about this? 

July 27, 2012, 10:10 

#8 
New Member
Christian Jungreuthmayer
Join Date: Mar 2009
Posts: 9
Rep Power: 8 

July 29, 2012, 04:53 

#9 
New Member
Tian Tang
Join Date: Jun 2012
Location: Copenhagen, Denmark
Posts: 16
Rep Power: 5 
hi, chritian,
thanks all the same Tian 

September 5, 2012, 16:05 

#10 
Member
Haomin Yuan
Join Date: Jan 2012
Location: Madison, Wisconsin, USA
Posts: 54
Rep Power: 5 
HI, All
As I read your guys' reply, you are couple two equations maximum. Is there any example that we can couple more equations? like 3 or 4? thank you in advance 

March 27, 2014, 10:05 

#11  
New Member
Join Date: Mar 2014
Posts: 20
Rep Power: 3 
Quote:
In which version of openfoam you compiled your solver? I want to solve some coupled transport equation and some one suggested me the blockCoupledScalarTransportFoam.C that is similar to your solver and I know this just work in extended version 3 and higher and not in version 2.2 

March 28, 2014, 03:18 
used version

#12 
New Member
Christian Jungreuthmayer
Join Date: Mar 2009
Posts: 9
Rep Power: 8 
Hello!
Well, that was a long long time ago when I worked on this solver I just checked my source code archive: I used OpenFoam1.6ext for this project. Christian 

March 30, 2014, 14:19 

#13 
New Member
Join Date: Mar 2014
Posts: 20
Rep Power: 3 
Thanks christian


April 6, 2014, 08:35 

#14  
New Member
Join Date: Mar 2014
Posts: 20
Rep Power: 3 
Quote:
I have the same problem somehow I want to solve dT/dt +d(Uk)/dx =0 dk/dt +d(UT)/dx =0 that are coupled and I have to use some tensors also thanks for the reply but in dimensionedTensor("jTensor", dimless, Foam::Tensor(0,1,0,1,0,0,0,0,0)) what is jTensor? because i have seen "0" somewhere I sthe "0" means 0 directory? 

April 24, 2014, 08:01 

#15 
New Member
Join Date: Mar 2014
Posts: 20
Rep Power: 3 
Does anyone still following this topic?
I want to solve a coupled equation dT/dt + d/dx(Uk) = 0 dk/dt + d/dx(UT) = 0 for solving that volTensorField j ( IOobject ( "j", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionedTensor( "jTensor", dimless,tensor(0,1,0,1,0,0,0,0,0) ) ); and solving fvVectorMatrix TEqn ( fvm::ddt(f) == ( j & fvc::div(phi,f) ) ); TEqn.solve(); runTime.write(); and defined f as volVectorField f ( IOobject ( "f", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); forAll(f,I) { f[I].component(0)=T[I]; f[I].component(1)=kappa[I]; f[I].component(2)=0; } my code compiled correctly but my results are unexpected I am in doubt with f,I mean I am not sure what kind of boundary conditions I have to implement for a vectorfield that including two scalarfield(in this case T and kappa) 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
icoLagrangianFoam OF1.6 myNewParticleSolver  heavy_user  OpenFOAM  16  February 11, 2012 06:15 
Error while running rhoPisoFoam..  nileshjrane  OpenFOAM Running, Solving & CFD  8  August 26, 2010 12:50 
ForcesCoeffs  ronaldo  OpenFOAM  4  September 14, 2009 07:11 
Problems with simulating TurbFOAM  barath.ezhilan  OpenFOAM  13  July 16, 2009 05:55 
MRFSimpleFoam amp cyclic patches  david  OpenFOAM Running, Solving & CFD  36  October 21, 2008 21:55 