
[Sponsors] 
June 10, 2009, 14:45 
Modifying the laplacian operator

#1 
New Member
Michael Lawson
Join Date: Apr 2009
Location: NREL  Boulder, CO
Posts: 11
Rep Power: 10 
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 xdirection. 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 

August 3, 2011, 16:27 

#2 
Senior Member
Mieszko Młody
Join Date: Mar 2009
Location: POLAND, USA
Posts: 142
Rep Power: 10 
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 

August 3, 2011, 16:54 

#3 
New Member
Michael Lawson
Join Date: Apr 2009
Location: NREL  Boulder, CO
Posts: 11
Rep Power: 10 
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 nonempty directions, the ydir 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 

August 3, 2011, 17:23 

#4 
Senior Member
Mieszko Młody
Join Date: Mar 2009
Location: POLAND, USA
Posts: 142
Rep Power: 10 
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 

August 4, 2011, 03:32 

#5 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,225
Rep Power: 23 
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).


August 4, 2011, 10:32 

#6 
New Member
Michael Lawson
Join Date: Apr 2009
Location: NREL  Boulder, CO
Posts: 11
Rep Power: 10 
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 

August 4, 2011, 10:49 

#7 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,225
Rep Power: 23 
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.


August 4, 2011, 12:13 

#8 
Senior Member
Mieszko Młody
Join Date: Mar 2009
Location: POLAND, USA
Posts: 142
Rep Power: 10 
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 ] (3e02 0 0 0 0 0 0 0 0 ); x,y and  directions diffusion: DTe DTe [ 0 2 1 0 0 0 0 ] (3e02 0 0 0 3e02 0 0 0 3e02 ); Thanks again! 

September 7, 2011, 11:51 

#9 
Senior Member
Matthias Voß
Join Date: Mar 2009
Location: Berlin, Germany
Posts: 446
Rep Power: 13 
hey there.
could anybody please give me a hint on where to find some sort of documentation for a part of the FSIsolver (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 

September 7, 2011, 12:05 

#10 
Senior Member
Mieszko Młody
Join Date: Mar 2009
Location: POLAND, USA
Posts: 142
Rep Power: 10 
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 

September 7, 2011, 12:06 

#11 
Senior Member
Matthias Voß
Join Date: Mar 2009
Location: Berlin, Germany
Posts: 446
Rep Power: 13 
okay.
found it. i guess. it´s a label for that part of the eq. to hook it up with the schemes in fvSchemes. 

September 7, 2011, 12:30 

#12 
Senior Member
Mieszko Młody
Join Date: Mar 2009
Location: POLAND, USA
Posts: 142
Rep Power: 10 
I think that all needed information you can find in the solver description:l
inearelastic, smallstrain 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 linearelastic 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) 

May 19, 2017, 10:35 

#13 
Member
Xinguang Wang
Join Date: Feb 2015
Posts: 42
Rep Power: 4 
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. 

May 19, 2017, 13:43 

#14 
Senior Member
Mieszko Młody
Join Date: Mar 2009
Location: POLAND, USA
Posts: 142
Rep Power: 10 
Hi Xinguang,
What do you need exactly ? Can you explain it with more details ? best ZMM 

May 20, 2017, 08:00 

#15 
Member
Xinguang Wang
Join Date: Feb 2015
Posts: 42
Rep Power: 4 
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 

May 20, 2017, 08:22 

#16 
Member
Xinguang Wang
Join Date: Feb 2015
Posts: 42
Rep Power: 4 
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)


May 20, 2017, 10:16 

#17 
Senior Member
Mieszko Młody
Join Date: Mar 2009
Location: POLAND, USA
Posts: 142
Rep Power: 10 
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) 

May 21, 2017, 07:00 

#18 
Member
Xinguang Wang
Join Date: Feb 2015
Posts: 42
Rep Power: 4 
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? 

May 21, 2017, 08:15 

#19 
Senior Member
Mieszko Młody
Join Date: Mar 2009
Location: POLAND, USA
Posts: 142
Rep Power: 10 
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. 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Elementwise multiplication operator  johndeas  OpenFOAM Running, Solving & CFD  2  May 14, 2012 09:45 
Question about the fvmatrix and Laplacian operator  liuhuafei  OpenFOAM Running, Solving & CFD  6  October 3, 2009 06:58 
Operator declaration in Thermophysical library  lena  OpenFOAM Running, Solving & CFD  0  March 12, 2009 10:47 
Material interfaces and the laplacian operator  cliffoi  OpenFOAM Running, Solving & CFD  8  November 8, 2006 09:57 
Material interfaces using the laplacian operator  cliffoi  OpenFOAM  0  November 6, 2006 11:42 