[buaad635]

Dear foamers,

i am new to openfoam programming, and use a codedFixedValue to realize a user-defined pressure outlet. And, the value of the pressure on the outlet need to be solved from linear equations.

So is there any functions that can be used to solve linear equations in the codedFixedValue boudary?

Thanks a lot!

[jherb]

I did something similar: I have a volVectorField, so for each cell a vector, which is the right side of the equation. And a volTensorField, so for each cell a tensor. Then your code would look somehow like this:

volVectorField rightSide(vector(1, 2, 3));
volTensorField leftSide(I); // I is the identity matrix as example

volVectorField solution(inv(leftSide) & rightSide); // using the inverted tensor (in each cell) to get the solution

This is pseudo code not tested. I guess, it should be possible to do the same with the vector and tensor class, but again, not tested.

[buaad635]

Originally Posted by jherb

I did something similar: I have a volVectorField, so for each cell a vector, which is the right side of the equation. And a volTensorField, so for each cell a tensor. Then your code would look somehow like this:

volVectorField rightSide(vector(1, 2, 3));
volTensorField leftSide(I); // I is the identity matrix as example

volVectorField solution(inv(leftSide) & rightSide); // using the inverted tensor (in each cell) to get the solution

This is pseudo code not tested. I guess, it should be possible to do the same with the vector and tensor class, but again, not tested.

Thank you, jherb! I will try this way soon!

[buaad635]

dear foamers,

I have solved this problem, if someone have similar questions, please refer to the following codes,

outlet1 //outlet
{
//type fixedValue;
//value uniform 13.0;
type codedFixedValue;
value uniform 0;
name BcPressureLB1;

code
#{
const surfaceScalarField& phi = db().lookupObject<surfaceScalarField>("phi");
//reading the flowRate on the output1 patch
label output1PatchID = this->patch().boundaryMesh().findPatchID("output1");
scalar phip1 = gSum(phi.boundaryField()[output1PatchID]);
//reading the flowRate on the output2 patch
label output2PatchID = this->patch().boundaryMesh().findPatchID("output2");
scalar phip2 = gSum(phi.boundaryField()[output2PatchID]);

Info << "total flow rate in output1 " << phip1 << " total flow rate in output2 " << phip2 << " m3/s" << endl;

//static scalar phipp;
const scalar k=0.6;
const scalar R1=k*0.255; //Rpca1
const scalar R2=k*0.255; //Rpca1
const scalar R3=k*4.73; //Rpca2
const scalar R4=k*4.73; //Rpca2
const scalar R5=k*29.1; //Rp
const scalar R6=k*29.1; //Rp
const scalar R7=k*16.05; //Rpcoa
const scalar R8=k*16.05; //Rpcoa
const scalar R9=k*1.56; //Rica
const scalar R10=k*1.56; //Rica
const scalar R11=k*15.49; //Rmca
const scalar R12=k*15.49; //Rmca
const scalar R13=k*22.4; //Rm
const scalar R14=k*22.4; //Rm
const scalar R15=k*0.428; //Raca1
const scalar R16=k*0.428; //Raca1
const scalar R17=k*9.64; //Raca2
const scalar R18=k*9.64; //Raca2
const scalar R19=k*43.7; //Ra
const scalar R20=k*43.7; //Ra
const scalar R21=k*3.21; //Racoa
const scalar Pv=0.0; //central venous pressure

SquareMatrix<scalar> AA(19, 0.0);
AA(0,0) = 1.0;
AA(0,2) = -1.0;
AA(1,1) = 1.0;
AA(1,3) = -1.0;
AA(2,2) = 1.0;
AA(2,8) = -(R3+R5);
AA(3,3) = 1.0;
AA(3,9) = -(R4+R6);
AA(4,2) = 1.0;
AA(4,4) = -1.0;
AA(4,10) = -R7;
AA(5,3) = 1.0;
AA(5,5) = -1.0;
AA(5,11) = -R8;
AA(6,4) = 1.0;
AA(6,14) = -R9*(R11+R13)/(R9+R11+R13);
AA(7,5) = 1.0;
AA(7,15) = -R10*(R12+R14)/(R10+R12+R14);
AA(8,4) = 1.0;
AA(8,6) = -1.0;
AA(8,12) = -R15;
AA(9,5) = 1.0;
AA(9,7) = -1.0;
AA(9,13) = -R16;
AA(10,7) = 1.0;
AA(10,6) = -1.0;
AA(10,18) = -R21;
AA(11,6) = 1.0;
AA(11,16) = -(R17+R19);
AA(12,7) = 1.0;
AA(12,17) = -(R18+R20);
AA(13,8) = 1.0;
AA(13,10) =1.0;
AA(14,9) = 1.0;
AA(14,11) = 1.0;
AA(15,10) = 1.0;
AA(15,12) = -1.0;
AA(15,14) = -1.0;
AA(16,11) = 1.0;
AA(16,13) = -1.0;
AA(16,15) = -1.0;
AA(17,12) = 1.0;
AA(17,18) = 1.0;
AA(17,16) = -1.0;
AA(18,13) = 1.0;
AA(18,17) = -1.0;
AA(18,18) = -1.0;

scalarField bb(19,1);
bb[0] = phip1*R1*1000000.0; //m3 to ml
bb[1] = phip2*R2*1000000.0;
bb[2] = Pv;
bb[3] = Pv;
bb[4] = 0.0;
bb[5] = 0.0;
bb[6] = Pv;
bb[7] = Pv;
bb[8] = 0.0;
bb[9] = 0.0;
bb[10] = 0.0;
bb[11] = Pv;
bb[12] = Pv;
bb[13] = phip1*1000000.0;
bb[14] = phip2*1000000.0;
bb[15] = 0.0;
bb[16] = 0.0;
bb[17] = 0.0;
bb[18] = 0.0;
//bb = -bb;

//Info << bb <<endl;

QRMatrix<scalarSquareMatrix> QR(AA);
scalarField xx(QR.solve(bb));
//scalarField xx;
Info << "P1 " << xx[0] << " P2 " << xx[1] << " mmHg" << endl;

//****************************************//
//pp += h*sum(phip); //test
(*this) == (xx[0]*133.3/1000.0);
#};

codeInclude
#{
#include "SquareMatrix.H"
#include "RectangularMatrix.H"
#include "QRMatrix.H"
#};
}