
[Sponsors] 
September 12, 2012, 05:54 
div(grad) not equal to laplacian close to a discontinuity

#1 
New Member
Join Date: Jun 2011
Posts: 13
Rep Power: 7 
Dear all,
I am solving a stress equilibrium equation for a solid, which can be written as: Code:
fvVectorMatrix UEqn ( fvm::laplacian( etaS*(2*mu + lambda), U) ==  fvc::div( stressSource) ) I included a pdf description of the setup for more details. http://dl.dropbox.com/u/34770732/201...gencefree.pdf Results converge very well. However, when I plot the resultant stress (stress = etaS*(2*mu + lambda)*fvc::grad(U) + stressSource), it is appparent that it is not divergence free. Comparing the residual of the equation vs. the divergence of the stress (which should be equal) I obtain very different behaviors, especially in the vicinity of the discontinuity, reflecting the fact that in that case Code:
fvc::laplacian( etaS*(2*mu + lambda), U) != fvc::div(etaS*(2mu+lambda)*fvc::grad(U)) Code:
*(2*mu + lambda)*fvc::grad(U)) while I did not expect an exact match of the discrete operators, I believe that the error I obtain is huge. On the opposite, on a smooth field (such as that obtained with a pure diffusion, e.g. basic laplacian tutorial) I do have fvc::laplacian = fvc::div(fvc::grad). I thus believe that my problem comes from the discontinuity at the boundary. I am using: gradSchemes: Gauss skewCorrected linear; divSchemes: Gauss skewCorrected linear; laplacianSchemes: Gauss linear corrected; I tried a few other schemes but haven't seen much difference. I am open to any hints and suggestions.  Should I use other spatial differentiation schemes (which one would you recommend)?  Are there some settings that should be modified in the vicinity of the boundary for a proper calculation of the stress? I can include the source code if it might be more helpful than the attached case description. Cheers, Diane 

September 12, 2012, 06:48 

#2 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 621
Rep Power: 22 
Hi Diane,
Does your mesh contain nonorthogonality and/or skewness? And do the boundary cells contain nonorthogonality? By default OpenFOAM sets boundary nonorthogonal corrections to zero, so stresses and displacements will be wrong in boundary cells unless you apply the appropriate corrections in your boundary conditions. I have developed displacement and traction BCs with nonorthogonal correction, you can get them in OpenFOAM1.6ext in the solidMechanics feature branch (see here). Also "extendedLeastSquares 0" will be much more accurate for the gradient scheme. Philip Last edited by bigphil; September 12, 2012 at 06:48. Reason: typo 

September 12, 2012, 08:40 

#3 
New Member
Join Date: Jun 2011
Posts: 13
Rep Power: 7 
Hi Philip,
Thank you for your super fast answer as always! It is true that my figures do not show the mesh. It is a basic Cartesian mesh generated with blockMesh, as such it should be perfectly orthogonal. Following your advice, I reverted my setting to extendedLeastSquare 0, but results are the same. Currently downloading the new extend solvers to see if I can find a solution there, but I believe that my traction boundary condition should actually be the same  as it was taken from the updatedLagrangian code before your new release). Do you have any other suggestion? in case it makes things easier, I uploaded the initialization code and a case file at http://dl.dropbox.com/u/34770732/201...continuity.zip Thank you again! Diane 

September 12, 2012, 08:55 

#4 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 621
Rep Power: 22 
Hi Diane,
OK, since your mesh is perfectly orthogonal, you don't need to use correct boundary conditions. Also the skewCorrected interpolation won't do anything and Gauss will be perfect for grad, div and laplacian. I don't entirely understand your problem, essentially you are saying that fvc::laplacian( etaS*(2*mu + lambda), U) should equal fvc::div(etaS*2mu+lambda)*fvc::grad(U)) but it doesn't..? Is this only near the boundary or everywhere? Mathematically both expressions are the same but numerically the laplacian is calculated using a compact molecule (i.e. each face uses just the cell centre either side) whereas the second expression calculates the full gradient using all the surrounding face neighbourneighbours and then takes the normal component. Philip 

September 13, 2012, 03:54 

#5 
New Member
Join Date: Jun 2011
Posts: 13
Rep Power: 7 
Hi Philip,
Thank you for the mathematical insights. I was thinking that there was a difference in the molecules used, because the difference is indeed located close to the boundary where the discontinuity is, which is not the case when the field is smooth (and I guess this stems from the fact that for smooth fields, gradients are continuous as well, so that the second derivative is more or less the same irrespective of the molecule used). I tried to attach a figure with more details (using a much coarser grid to better show the difference). I am plotting fvc::laplacian(DU) vs. fvc::div(fvc::grad(DU)) This run uses a simpler formulation, so that the stress is equal to: Dstress = 2*mu*strain(DU) + lambda*div(DU)*I  alpha*P At t=0, the material is suddenly loaded with pressure, while at the boundary fluid is free to flow (in the porous media) so that the pressure at the boundary is 0. In order to have a divergence free stress, we would thus need a DU field that has a constant gradient except at the boudary where it should compensate for the fact that P=0. As you can see the gradient is constant throughout most of the domain, but I guess it starts to change too far away from the boundary. Laplacian(DU) on the other hand does make sense. I agree that this must be due to the difference in computational molecule and the presence of the discontinuity at the boundary. Are there some cases where you had discontinuous fields (shocks?), are there some special treatments of the derivatives I should include to improve the coherence of my spatial operators when I get close to the discontinuity? Thank you again for your time, Diane 

September 13, 2012, 15:32 

#6 
New Member
Join Date: Jun 2011
Posts: 13
Rep Power: 7 
Hi Philip,
Thinking back to the computational molecule, it completely makes sense that this results in the observed error in the vicinity of the pressure discontinuity. However, since I am doing an updated Lagrangian approach I would need my stress to be more accurate than that. The place where both grad(U) and p are most accurate for me are on the cell faces. I am thus thinking that if I could:  compute grad(U) on the cell faces  have p on the cell faces (using either a staggered data storage,or interpolating the cell centered values)  compute my stress on the cell faces (something like mu*gradUalpha*p)  and then only interpolate back to the cell centers then I believe the accuracy of the computed stressed would actually be greater(especially since grad(U) would then be coherent with the laplacian molecule) Is there a way to perform these computations on the cell faces, without interpolating grad(U) back to the cell centers? Thank you again, Diane 

September 14, 2012, 11:48 

#7 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 621
Rep Power: 22 
Hi Diane,
Hmnnn, I am not sure how to solve your problem. I am still not entirely clear on your problem. I can see that you find "div(grad(DU)) != laplacian(DU)", did you try write out by hand why these are different in your case? Philip 

September 18, 2012, 06:23 

#8 
New Member
Join Date: Jun 2011
Posts: 13
Rep Power: 7 
Hi Philip,
Sorry for the late reply. Yes, it is through the comparison of the values I was getting to the ones I had expected from a computation by hand, that I started wondering about what was really happening. I think I have it figured out. It took a little while for me to figure out the proper functions to use to go from cell centers to cell faces and vice versa. But computing stresses on the faces and then div(stress) at the cell centers seems to work just fine. I also looked at your solid modules and what I am doing kind of goes along the same lines as your "decompose" approach for the source term in the stress equation. I will now run a few checks and cases to verify that I haven't spoken too quickly. Cheers, Diane 

September 18, 2012, 07:51 

#9 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 621
Rep Power: 22 
OK great,
Let us know if you figure something out! Philip 

September 19, 2012, 09:30 

#10 
Senior Member
Join Date: Nov 2010
Posts: 113
Rep Power: 8 
Hi,
if I may dig in here: Isn't it the same as the calculation of the curvature in interFoam? What we need is div(grad(alpha)) and, as you said, it is calculated on the faces. So maybe it will help you to compare the procedure to your case?! (You find it in interfaceProperties.C) 

December 6, 2012, 06:58 

#11  
New Member
Dennis Arreborg Hansen
Join Date: Mar 2012
Location: Denmark
Posts: 1
Rep Power: 0 
Quote:
Would you please elaborate a bit on why you belive that the gradients are best defined on the cell faces and subsequently clarify how you would define them on the cell faces. IMO, an interpolation will newer increase the accuracy of a value. Kind regards Dennis Arreborg 

Thread Tools  
Display Modes  

