CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   Modifying the laplacian operator (

mlawson June 10, 2009 14:45

Modifying the laplacian operator
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:

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


ziemowitzima August 3, 2011 16:27

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.


mlawson August 3, 2011 16:54

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.


ziemowitzima August 3, 2011 17:23

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 August 4, 2011 03:32

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 August 4, 2011 10:32

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?


akidess August 4, 2011 10:49

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. and recompile.

ziemowitzima August 4, 2011 12:13

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 September 7, 2011 11:51

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,

ziemowitzima September 7, 2011 12:05

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.


mvoss September 7, 2011 12:06

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

ziemowitzima September 7, 2011 12:30

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)


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:


if (thermalStress)

JasonWang3 May 19, 2017 10:35

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 May 19, 2017 13:43

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


JasonWang3 May 20, 2017 08:00

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.


JasonWang3 May 20, 2017 08:22

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 May 20, 2017 10:16

Laplacian does not equal grad(grad),
the true statment is:
\nu\nabla^2(T) = \nabla\cdot(\nu\nabla(T))
it OF reads:
Laplacian(T) = div(grad(T))

so you could:
gradT = grad(T)
T_xx = div(gradT)

gradT = grad(T)
T_yy = div(gradT)

gradT = grad(T)
T_zz = div(gradT)

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:

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:

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:

JasonWang3 May 21, 2017 07:00

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=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 May 21, 2017 08:15

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.

All times are GMT -4. The time now is 15:47.