CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

div(grad) not equal to laplacian close to a discontinuity

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   September 12, 2012, 05:54
Default div(grad) not equal to laplacian close to a discontinuity
  #1
New Member
 
Join Date: Jun 2011
Posts: 13
Rep Power: 5
gtg258f is on a distinguished road
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) )
using PCG with a tolerance and relative tolerance of 1e-12 and 0, respectively. In my problem, the stressSource has a discontinuity.
I included a pdf description of the set-up for more details.
http://dl.dropbox.com/u/34770732/201...gence-free.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
gtg258f is offline   Reply With Quote

Old   September 12, 2012, 06:48
Default
  #2
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 531
Rep Power: 17
bigphil will become famous soon enough
Hi Diane,

Does your mesh contain non-orthogonality and/or skewness? And do the boundary cells contain non-orthogonality?

By default OpenFOAM sets boundary non-orthogonal 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 non-orthogonal correction, you can get them in OpenFOAM-1.6-ext 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
bigphil is offline   Reply With Quote

Old   September 12, 2012, 08:40
Default
  #3
New Member
 
Join Date: Jun 2011
Posts: 13
Rep Power: 5
gtg258f is on a distinguished road
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
gtg258f is offline   Reply With Quote

Old   September 12, 2012, 08:55
Default
  #4
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 531
Rep Power: 17
bigphil will become famous soon enough
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 neighbour-neighbours and then takes the normal component.

Philip
bigphil is offline   Reply With Quote

Old   September 13, 2012, 03:54
Default
  #5
New Member
 
Join Date: Jun 2011
Posts: 13
Rep Power: 5
gtg258f is on a distinguished road
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
Attached Images
File Type: jpg Laplacian-vs-DivGrad.jpg (52.7 KB, 49 views)
gtg258f is offline   Reply With Quote

Old   September 13, 2012, 15:32
Default
  #6
New Member
 
Join Date: Jun 2011
Posts: 13
Rep Power: 5
gtg258f is on a distinguished road
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*gradU-alpha*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
gtg258f is offline   Reply With Quote

Old   September 14, 2012, 11:48
Default
  #7
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 531
Rep Power: 17
bigphil will become famous soon enough
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
bigphil is offline   Reply With Quote

Old   September 18, 2012, 06:23
Default
  #8
New Member
 
Join Date: Jun 2011
Posts: 13
Rep Power: 5
gtg258f is on a distinguished road
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
gtg258f is offline   Reply With Quote

Old   September 18, 2012, 07:51
Default
  #9
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 531
Rep Power: 17
bigphil will become famous soon enough
OK great,

Let us know if you figure something out!

Philip
bigphil is offline   Reply With Quote

Old   September 19, 2012, 09:30
Default
  #10
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 5
lindstroem is on a distinguished road
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)
lindstroem is offline   Reply With Quote

Old   December 6, 2012, 05:58
Default
  #11
New Member
 
Dennis Arreborg Hansen
Join Date: Mar 2012
Location: Denmark
Posts: 1
Rep Power: 0
Dennis A is on a distinguished road
Quote:
Originally Posted by gtg258f View Post
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*gradU-alpha*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
Dear Diane

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
Dennis A is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 11:22.