[mlawson]

Hi everyone,

I would like to modify the laplacian so the derivative is only calculated in one direction - that is I would like to compute only d^2T/dx^2 and not d^2T/dx^2 + d^2T/dy^2 + d^2T/dz^2

Then I want to solve a system of equations using this modified laplacian. I tried to do this by dotting the result of the laplacian operator with a unit normal vector in the x-direction. A simple example would be:

solve
(
fvm::ddt(T)
- fvm::laplacian(DT, T) & xDir // where xDir is a dimentionedVector with the units [0 1 -1 0 0 0 0] (1 0 0)
);


Unfortunately, I get the following error when trying to compile my solver.

---------------------------------------
1dScalarTransportFoam.C: In function ‘int main(int, char**)’:
1dScalarTransportFoam.C:67: error: no match for ‘operator&’ in ‘Foam:perator-(const Foam::tmp<Foam::fvMatrix<Type> >&, const Foam::tmp<Foam::fvMatrix<Type> >&) [with Type = double](((const Foam::tmp<Foam::fvMatrix<double> >&)((const Foam::tmp<Foam::fvMatrix<double> >*)(& Foam::fvm::laplacian(const Foam::dimensioned<Type2>&, Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&) [with Type = double, GType = double](((Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&)(& T))))))) & xDir’
---------------------------------------

I assume I get this error because the & operator is not defined for a fvScalarMatrix.

Does anyone know a method modifying the laplacian operator so that the derivative is only calculated in one direction?

-Mike

[ziemowitzima]

Hi,
I am trying to do the same thing.
Did you manage to do this ?
Method you described is wrong, because laplacian(DT,T) is not a vector , but scalar:
Txx + Tyy + Tzz, so you cant make dot product operation for scalar and vector...
I did it in explicit way, but this not the best way...

If you have some updates how to solve this problem, pleas let me know.

best
-ZM

[mlawson]

You are correct the result of the laplacian operator is a scalar in this case. What should have asked is, is there a way to dot the laplacian operator itself with a unit vector before it is applied to the variable T (i.e (d/dx^2 + d/dy^2 + d/dz^2) dot (1 0 0) * T).

Either way, I never did solve this problem in a robust way. I ended up tricking the solver into only computing diffusion in one direction by making my two dimensional domain very long in one of the non-empty directions, the y-dir in my case, so the d/dy^2 term was negligible with respect to the d/dx^2 term. I know this is not a good way to solve the problem, but it ended up working for the specific case I was working on.

Good luck and please let me know if you find a solution to the problem.

-Mike

[ziemowitzima]

in explicit way u can do this like that:
volVectorField gradT=fvc::grad(T);
volScalarField gradTx = gradT & vector(1,0,0);
volScalarField Tx = gradT.component(0); // T_x
volVectorField gradT2 = fvc::grad(Tx); // (T_xx, T_xy)
volScalarField Txx = gradT2.component(0); // T_xx

[akidess]

Michael, it sounds like you want anisotropic diffusion - why not just pass a diffusion coefficient tensor to laplacian instead of a scalar? I don't think you need any modifications to the code to do that (but I haven't tried it myself).

[mlawson]

Hi Anton,

Mathematically that makes sense, but I'm unsure how to convert the scalar field T into a tensor in OpenFOAM. Do you have an idean about how this would be done?

-Mike

[akidess]

T stays a scalar field, you only have to change the diffusion coefficient into a tensor (you wouldn't do anything else on pen and paper). Redefine DT from dimensionedScalar to dimensionedTensor (see e.g. http://www.foamcfd.org/Nabla/guides/...sGuidese5.html) and recompile.

[ziemowitzima]

Hi Akidess,
Thanks for your suggestion.
It works just as it should.
Human being is learning all its life...

To Mlawson:
x - direction diffusion:
DTe DTe [ 0 2 -1 0 0 0 0 ] (3e-02 0 0 0 0 0 0 0 0 );

x,y and - directions diffusion:
DTe DTe [ 0 2 -1 0 0 0 0 ] (3e-02 0 0 0 3e-02 0 0 0 3e-02 );

Thanks again!

[mvoss]

hey there.

could anybody please give me a hint on where to find some sort of documentation for a part of the FSI-solver (icoFSIFoam??) like mentioned above. i am espacially interessted in how the pde´s have been derived for the solid/stress part.

e.g. fvm::laplacian (2*mu, + lambda, Usolid, "laplacian (DU,U)")

what is that third part good for? I wasn´t able to get any doc. for a third part in laplacian(..).
Is this some sort of comment??

Thanks in advance,
neewbie

[ziemowitzima]

Third part "laplacian(DU,U)" is just needed for fvScheme file, where laplacian(DU,U) is defined. Here it just says that variable DU = 2*mu + lambda.

ZMM

[mvoss]

okay.
found it. i guess.
it´s a label for that part of the eq. to hook it up with the schemes in fvSchemes.

[ziemowitzima]

I think that all needed information you can find in the solver description:l
inear-elastic, small-strain deformation of a solid body, with optional thermal
diffusion and thermal stresses.
For example in solver solidDisplacementFoamfact there are only two PDE solved (I dont know where FSIsolver is):

fvm::ddt(T) == fvm::laplacian(DT, T)

and

fvm::d2dt2(D)
==
fvm::laplacian(2*mu + lambda, D, "laplacian(DD,D)")
+ divSigmaExp

first equation is just standart heat equation,
second equation is just linear-elastic equation with optional thermat stress (+ divSigmaExp) influence (calculated by solving first equation).
This are standart equations so they derivation can be easly found in basic PDE books, or online:
http://en.wikipedia.org/wiki/Linear_elasticity

ZMM



if (thermalStress)

[JasonWang3]

I still confused after reading all the thread. I need to separate Txx, Tyy and Tzz as mentioned.

Is there anyone give some example to do it?

Thanks a lot.

[ziemowitzima]

Hi Xinguang,
What do you need exactly ?
Can you explain it with more details ?

best
ZMM

[JasonWang3]

hi ZMM

I want to calculate the second derivatives of a scalar field such as Temperature separately. For example, d^2T/d^2x, d^2T/d^2y and d^2T/d^2z, and then these three terms sum up as d^2T/d^2x + d^2T/d^2y + d^2T/d^2z = laplacian(T). I need to use each term for the calculation.

Many thanks for your reply.

Jaosn

[JasonWang3]

I also tried to use grad(grad) to calculate d^2T/d^2x d^2T/d^2y d^2T/d^2z , but when I add up these three terms, it's not equal to laplacian(T)

[ziemowitzima]

Laplacian does not equal grad(grad),
the true statment is:

it OF reads:
Laplacian(T) = div(grad(T))

so you could:
gradT = grad(T)
gradT.component(1)=0
gradT.component(2)=0
T_xx = div(gradT)

gradT = grad(T)
gradT.component(0)=0
gradT.component(2)=0
T_yy = div(gradT)

gradT = grad(T)
gradT.component(0)=0
gradT.component(1)=0
T_zz = div(gradT)

OR:
you could do as it is explained above:
and define the diffusivity coeff as a tensor:

for x - direction diffusion only:
DTx DTx [ 0 2 -1 0 0 0 0 ] (1 0 0 0 0 0 0 0 0 );

in the code:
T_xx=fvc::Laplacian(DTx,T)

for y - direction diffusion only:
DTy DTy [ 0 2 -1 0 0 0 0 ] (0 0 0 1 0 0 0 0 0 );

in the code:
T_yy=fvc::Laplacian(DTy,T)

for z - direction diffusion only:
DTz DTz [ 0 2 -1 0 0 0 0 ] (0 0 0 0 0 0 1 0 0 );

in the code:
T_zz=fvc::Laplacian(DTz,T)

[JasonWang3]

In my code Ux is a scalar field denotes the streamwise flow velocity.


const volScalarField du2dy2=fvc::laplacian(Eiy,Ux);
where Eiy is (0,0,0,0,1,0,0,0,0)

const volVectorField gradUx=fvc::grad(Ux);
const volScalarField dudy=gradUx.component(vector::Y);
const volVectorField graddudy=fvc::grad(dudy);

const volScalarField divgradUx=fvc::div(gradUx);
(for 2D flat plate flow dudx is very small)

when I compare:

term1= du2dy2
term2=graddudy.component(vector::Y)
term3=divgradUx

term2=term3(nearly the same), but term1 is not equal to term2 or term3 (30% less)

In fvSchemes, the default scheme are all Gauss linear.

Do I miss something to make such difference happens?

[ziemowitzima]

in your case :
term3 = Ux_xx + Uy_yy + Ux_yz,

so it wont be same as term1,

to have ti equal to term1, you should first:
const volVectorField gradUx=fvc::grad(Ux);
gradUx.component(vector::X) = 0;
gradUx.component(vector::Z) = 0;

then you would have:
gradUx = (0,Ux_y,0),

an then applying:
const volScalarField divgradUx=fvc::div(gradUx);
you wpould have:
(_x,_y,_z)*(0,Ux_y,0) = Ux_yy

and term3 sould be same as term1

I am not sure about term2.

[Elham]

Dear Foamers,

In a multiphase fluid type, I have the following form of TEqn

    fvScalarMatrix TEqn
    (
          fvm::ddt(T)
        + fvm::div(phi,T)
        - Foam::fvm::laplacian( k/rhoCp , T)
    );


    TEqn.relax();

        solve
         (
            TEqn
            ==
            fvc::ddt(p)*(1/rhoCp)
            +Foam::fvc::laplacian( D23*T , XV)
         );

where XV is vapor mass fraction as following:

rhogp= (rho2*limitedAlpha2+rho3*limitedAlpha3)/(limitedAlpha2+limitedAlpha3);
XV= rho2/rhogp

with the above equation I got too unbounded alpha fields at very first step .If I manipulate TEqn as the following:

fvScalarMatrix TEqn
    (
          fvm::ddt(T)
        + fvm::div(phi,T)
        - Foam::fvm::laplacian( k/rhoCp , T)
        - Foam::fvm::laplacian( D23*T , XV)
    );


    TEqn.relax();

        solve
         (
            TEqn
            ==
            fvc::ddt(p)*(1/rhoCp)

         );

It will blow up at the first time step with the following error:

gradientInternalCoeficients cannot be called for a calculatedFvPatchField

I will apprciate any hint for fixing.

Cheers,

Elham