CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)
-   -   Laplacian(A*B) (https://www.cfd-online.com/Forums/openfoam/122138-laplacian-b.html)

Ricardo August 13, 2013 02:28

Laplacian(A*B)
 
Hello, is possible solve in openFoam a function laplacian(k,A*B) where k is constant, A is known volScalarField and B is searched unknown (volScalarField)?

Thanks
Richard

Bernhard August 13, 2013 03:42

Did you consider fvm::laplacian(k*A,B)?

Ricardo August 13, 2013 06:22

Yes I did, but I don't know if it is correct from mathematical point of view.
I need solve this:

dF/dt=laplacian(delta,P)
P is partial pressure and could be rewrite F = P / Psat
Psat is partial saturated pressure (function of Temperature)

I tried this, but without succes... :(
fvm::ddt(dw,F)=fvm::laplacian(delta*Psat,F)
fvm::ddt(dw,F)=fvc::laplacian(delta,Psat*F)

Thanks
Richard

niklas August 13, 2013 06:50

when you say tried, but without success....what does that mean?

Could you not compile it?
Could you not run it?
Did it not yield the desired result?

fvm::ddt(dw,F)=fvm::laplacian(delta*Psat,F)

this should work and be your best bet...assuming Psat is just a constant.

just a note, but remember that
laplacian(k,AB) = div(k * grad(AB))
and is not, in general, equal to
laplacian(kA,B) = div(kA * grad(B))
only if A = const.

niklas August 13, 2013 06:53

also
= is an assignement
== is checking if its equal

so in that implementation you should use ==

Bernhard August 13, 2013 07:00

Quote:

Originally Posted by niklas (Post 445406)
== is checking if its equal

Not always: http://www.cfd-online.com/Forums/ope...-operator.html

niklas August 13, 2013 07:05

Quote:

Originally Posted by Bernhard (Post 445410)

true :)
sorry, should have clarified that I ment in this example

Ricardo August 13, 2013 07:11

I have operator "==" in code, it is right, compile and solution done, but result is wrong, I think.

I write solver for Heat and Moisture transport in porous materials.
1 step - Temperature calculation (changing thermal parameters with moisture)
2 step - Humidity calculation (changing moisture parameters with temperature)
next Iterations...

Saturated Pressure is not a constant, it depends on temerature and the temperature is changing in time. Transport parameters are calculated in time and stored in volScalarFields.

Thanks
Richard

immortality August 13, 2013 07:32

Hi Ricardo
how you have solved this equation?! with which solver of OF?

Ricardo August 13, 2013 07:36

I started from basic Laplacian example, and still using PCG solver with DIC preconditioner, is it wrong?

immortality August 13, 2013 07:46

No I just wanted to know.you are modifying the laplacianFoam,right? could you send me your case?

Ricardo August 13, 2013 07:51

Quote:

Originally Posted by immortality (Post 445425)
No I just wanted to know.you are modifying the laplacianFoam,right? could you send me your case?

It is not problem, but it doesn't work correctly now, ...

immortality August 13, 2013 07:56

I just like to know the way you are doing it.no problem.maybe an idea occurred to me.
my email is: force.of.love@gmail.com if you want to send it.

Ricardo August 13, 2013 08:10

Quote:

Originally Posted by immortality (Post 445429)
I just like to know the way you are doing it.no problem.maybe an idea occurred to me.
my email is: force.of.love@gmail.com if you want to send it.

It has been send.

I am Ph.D. student, I would like use this program for simulations in my thesis. I wrote similar algorithm in 1D in Visual Basic, it works without any problem. OpenFOAM and C++ is new for me...
Richard

fabian_roesler August 13, 2013 16:45

Dont forget the chain rule
 
Hi,

when psat is a function of temperature, you are not allowed to pull it out of the gradient in the laplacian function. But you could use the chain rule to get an implicit term of F and an explicit term of psat. Here we go:

Eqn (1): This is your starting point
Code:

    fvm::ddt(dw,F)
  + fvm::laplacian(delta,psat*F)

Eqn (2): Split the gradient term
Code:

    fvm::ddt(dw,F)
  + fvc::div(delta*F*fvc::grad(psat) + delta*psat*fvc::grad(F))

Eqn (3): Split the divergence term
Code:

    fvm::ddt(dw,F)
  + fvc::div(delta*F*fvc::grad(psat))
  + fvm::div(delta*psat*fvc::grad(F))

Eqn (4): form the last two terms into laplacians again
Code:

    fvm::ddt(dw,F)
  + fvc::laplacian(delta*F,psat)
  + fvm::laplacian(delta*psat,F)

Eqn 4 should work for you. Try it out.

Regards

Fabian

immortality August 13, 2013 18:02

Code:

dF/dt=laplacian(delta,P)
P is partial pressure and could be rewrite F = P / Psat
Psat is partial saturated pressure (function of Temperature)

I tried this, but without succes...
fvm::ddt(dw,F)=fvm::laplacian(delta*Psat,F)
fvm::ddt(dw,F)=fvc::laplacian(delta,Psat*F)

whats dw in these equations while it isn't in dF/dt?
Code:

Eqn (1): This is your starting point
Code:
    fvm::ddt(dw,F)
  + fvm::laplacian(delta,psat*F)
Eqn (2): Split the gradient term
Code:
    fvm::ddt(dw,F)
  + fvc::div(delta*F*fvc::grad(psat) + delta*psat*fvc::grad(F))
Eqn (3): Split the divergence term
Code:
    fvm::ddt(dw,F)
  + fvc::div(delta*F*fvc::grad(psat))
  + fvm::div(delta*psat*fvc::grad(F))

how and why implicit terms convert to explicit ones and vice versa?

fabian_roesler August 14, 2013 02:43

Hi

The first laplacian

Code:

  + fvc::laplacian(delta*F,psat)
accounts for the change in saturation pressure and would be zero for constant psat. It has to be treated explicit (fvc), as the conservation equation is for F. The second term

Code:

  + fvm::laplacian(delta*psat,F)
can be treated implicit (fvm) now due to the chain rule.

Regards

Fabian

immortality August 14, 2013 05:07

Quote:

laplacian(k,AB) = div(k * grad(AB))
Hi Niklas, laplacian(A)=d^2A/dx^2+d^2A/dy^2 in 2D but what means above laplacian(k,AB)?if they are multiply why laplacian(k,AB) is different from laplacian(kA,B)?

Ricardo August 14, 2013 09:56

Quote:

Originally Posted by fabian_roesler (Post 445542)
Hi,
Code:

    fvm::ddt(dw,F)
  + fvc::laplacian(delta*F,psat)
  + fvm::laplacian(delta*psat,F)

Eqn 4 should work for you. Try it out.

Regards

Fabian


Thank you very much, I am going to try it. :)
Yesterday I used other approach (similar which I using in Visual Basic) for implicit scheme and I think that results was correct.

- generate fvscalarmatrix "mat" from laplacian(delta,F)
- multiply each cell in fvscalarmatrix mat by Psat in point of this cell (using lduAddressing)
- solve (fvm::ddt(dw,F)==mat)

But this code is more complicated and your approach seems better.

dw is aproximation function of sorption curve (relation between water content and relative humidity of air) it is nonlinear function, I forgot write in equation:o

Richard

Ricardo August 14, 2013 10:47

The part "fvc::laplacian(delta*F,Psat)" causes stability problem,... :(

volScalarField F

internalField nonuniform List<scalar>
4000
(
0.503748
0.500001
0.5
0.5
0.5
0.5
0.49999
0.44179
-361.167
-2.31903e+06
0.500001
0.5
0.5
0.5
0.5
0.5
0.49999
0.44179
-361.165
-2.31902e+06

niklas August 14, 2013 12:22

you could try to implement it implicitly
fvc::laplacian(delta*F,Psat) = fvc::div( delta*F*grad(Psat) )

volVectorField dps = delta * fvc::grad(Psat);
surfaceScalarField dpsf = fvc::interpolate(dps) & mesh.Sf();

then you can view the dps field as a 'convective' change of F due to the Psat gradient
and rewrite it like this

fvc::laplacian(delta*F,Psat) = fvm::div(dpsf, F)


...should work

fabian_roesler August 14, 2013 13:35

Hi

That's a pitty. But the code of niklas seems to be promissing. Ricardo, I didn't get your approach in total. Could you explain it in a little more detail? Especially the part of the lduAddressing.

Regards

Fabian

Ricardo August 14, 2013 14:21

Quote:

Originally Posted by niklas (Post 445771)
fvc::laplacian(delta*F,Psat) = fvm::div(dpsf, F)
...should work

implementation to solver >> solve (fvm::ddt(dw,F)==fvm::laplacian(lambdaD*Psat,F)+fv m::div(dpsf, F)); doen't works, too :(

Quote:

Originally Posted by fabian_roesler (Post 445791)
Hi

That's a pitty. But the code of niklas seems to be promissing. Ricardo, I didn't get your approach in total. Could you explain it in a little more detail? Especially the part of the lduAddressing.

Regards

Fabian

This is my code, but it I am still working ...

a part of code ...
Code:

    float tdelta = 1.0;
    float tdeltaN = 1.0;
   
    float fdelta = 1.0;
    float fdeltaN = 1.0;

    double relFaktor = 1.0;
    int iter = 0;
   
    #include "boundaryCalculation.H"
   
    do
    {

    Tt = T;
   
    #include "heatParameters.H"

    solve
            (
              fvm::ddt(cap,T)==fvm::laplacian(lambda, T)//+hv*fvc::laplacian(lambdaD, Psat*F)
        );

    T=(relFaktor)*T+(1-relFaktor)*Tt;
   
    tdeltaN=diffFields(T,Tt);

    #include "saturatedpressure.H"
    #include "moistureParameters.H"


    fvScalarMatrix mat1
    (
          fvm::laplacian(lambdaD,F)
    );

   
    label Nc = mesh.nCells(); //Total number of cells
   
    // Diagonal coefficients calculation
        for(label i=0; i<Nc; i++)
        {
          mat1.diag()[i]*=Psat[i];
        }
       
    //Off-diagonal coefficients calculation
        for(label faceI=0; faceI<mat1.lduAddr().lowerAddr().size(); faceI++)
        {
    mat1.lower()[faceI]*=Psat[mesh.faceOwner()[faceI]];
      mat1.upper()[faceI]*=Psat[mesh.faceNeighbour()[faceI]];
        }


      dimensionedScalar Pa ("Pa",dimensionSet(1,-1,-2,0,0,0,0), 1.0 );

    mat1*=Pa;

    solve
        (
            fvm::ddt(dw,F)==mat1
    //fvm::laplacian(lambdaD*Psat,F)//+fvc::laplacian(dT,T)// // 
        );
   
    F=(relFaktor)*F + (1-relFaktor)*Ft;

    fdeltaN=diffFields(F,Ft);
   
    #include "w_prepocet.H"
    iter++;
   
    if (iter==1)
    {
    tdelta=tdeltaN;
    fdelta=fdeltaN;
    }

    if ((tdeltaN>tdelta) || (fdeltaN>fdelta))
    {
    relFaktor=relFaktor/2;
    T=(relFaktor)*T+(1-relFaktor)*Tt;
    F=(relFaktor)*F + (1-relFaktor)*Ft;
    }
   
    } while (((tdeltaN>maxDeltaT) || (fdeltaN>maxDeltaH)) && iter < maxIter);

my approach to laplacian(labmbdaD,F*Psat) is in code, maybe for expert it could be ridiculous :D , but I am working with openFoam and C++ 2 weeks, and this code gives some "expected" values. Other approach gives unstability or doesn't work...

I tried next approach, rewrite partial pressure by gradient of temperature, it seems, that gives similar values like lduaddressing script :)

I don't know where is problem, I have to check all functions, material parameters and other properties again,...

thanks to all for support and advices
Richard

niklas August 14, 2013 15:11

you REALLY need to be more specific when you say 'it doesnt work'.

Ricardo August 14, 2013 15:39

I was wrong :o, It works, but obtained values ...

dimensions [0 0 0 0 0 0 0];

internalField nonuniform List<scalar>
4000
(
-1.35814e-10
4.58031e-11
4.58031e-11
4.58031e-11
4.5803e-11
4.58031e-11
4.58031e-11
4.58031e-11
-1.82592e-10
-1.77054e-12
-2.23032e-11

So it seems (two scheme unstable), that a problem could be somewhere else,

Is right using parameters (lambda), stored as volScalaField and (ununiform field) calculated in each step equation by this way? (fvm::laplacian(lambda,T)) Is it interpreted correctly to scheme or is used average values?

Thanks
Richard

quarkz April 14, 2020 01:13

Quote:

Originally Posted by fabian_roesler (Post 445542)
Hi,

when psat is a function of temperature, you are not allowed to pull it out of the gradient in the laplacian function. But you could use the chain rule to get an implicit term of F and an explicit term of psat. Here we go:

Eqn (1): This is your starting point
Code:

    fvm::ddt(dw,F) + fvm::laplacian(delta,psat*F)
Eqn (2): Split the gradient term
Code:

    fvm::ddt(dw,F) + fvc::div(delta*F*fvc::grad(psat) + delta*psat*fvc::grad(F))
Eqn (3): Split the divergence term
Code:

    fvm::ddt(dw,F) + fvc::div(delta*F*fvc::grad(psat)) + fvm::div(delta*psat*fvc::grad(F))
Eqn (4): form the last two terms into laplacians again
Code:

    fvm::ddt(dw,F) + fvc::laplacian(delta*F,psat) + m::laplacian(delta*psat,F)
Eqn 4 should work for you. Try it out.

Regards

Fabian

Hi, anyone manage to solve the problem?

I face similar problem whereby I need to convert

Code:

fvm::laplacian((alpha*rho*DkEff(F1)), k_)
to
Code:

fvm::laplacian((alpha*DkEff(F1)), rho*k_)
I used eqn4 as recommended to get:
Code:

fvm::laplacian((alpha*rho*DkEff(F1)), k_) + fvc::laplacian((alpha*k_*DkEff(F1)), rho)
but seems like I can't place rho there in:
Code:

fvc::laplacian((alpha*k_*DkEff(F1)), rho)
The error is "error: no matching function for call to ‘laplacian(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >, const rhoField&)’"

Another error is "note: template argument deduction/substitution failed"

Any idea how to solver it? Must been trying to solve it for a month.

Thanks


All times are GMT -4. The time now is 23:13.