CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Solving Coupled Set of PDEs (http://www.cfd-online.com/Forums/openfoam-solving/88729-solving-coupled-set-pdes.html)

 nickmai123 May 24, 2011 15:16

Solving Coupled Set of PDEs

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) }
However, the inner product is not defined between a primitive tensor and a fvMatrix. Does anyone here have any recommendations? I need to use the fvm:: namespace because it does implicit discretization of the Laplacian.

EDIT: I've found coupledFvScalarMatrix in the OpenFOAM 1.6-ext 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

 chegdan May 24, 2011 16:05

try blockCoupledScalarTransportFoam

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)             );
With some additional methods to form a block matrix. Maybe that will help.

Dan

 alberto May 25, 2011 02:15

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))  );
where "dimless" should become

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,

 Tron November 24, 2011 09:09

blockCoupledScalarTransportFoam - coupling via operator

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 off-diagonal 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

 tiat July 26, 2012 10:41

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 :(

 Tron July 27, 2012 01:20

Hi Tian,

Quote:
 Originally Posted by tiat (Post 373728) is there anybody still following this topic?
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

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;

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 off-diagonal 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();
}

# 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

 tiat July 27, 2012 02:39

hi, christian,

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!).

 Tron July 27, 2012 10:10

Quote:
Unfortunately not. It is quite a while ago and I didn't work on this project since then. Sorry.

 tiat July 29, 2012 04:53

hi, chritian,

thanks all the same :)

Tian

 yhaomin2007 September 5, 2012 16:05

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?

 alisina-s March 27, 2014 10:05

Quote:
Hi christan
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

 Tron March 28, 2014 03:18

used version

Hello!

Well, that was a long long time ago when I worked on this solver :) I just checked my source code archive: I used OpenFoam-1.6-ext for this project.

Christian

 alisina-s March 30, 2014 14:19

Thanks christian:)

 alisina-s April 6, 2014 08:35

Quote:
 Originally Posted by alberto (Post 309128) 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))  ); where "dimless" should become 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,
Hi
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?

 alisina-s April 24, 2014 08:01

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_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::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)

 All times are GMT -4. The time now is 11:33.