[gtg258f]

Dear all,

I am solving a stress equilibrium equation for a solid, which can be written as:

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

fvc::laplacian( etaS*(2*mu + lambda), U) 
!=
fvc::div(etaS*(2mu+lambda)*fvc::grad(U))

*(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

[bigphil]

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

[gtg258f]

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

[bigphil]

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

[gtg258f]

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

[gtg258f]

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

[bigphil]

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

[gtg258f]

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

[bigphil]

OK great,

Let us know if you figure something out!

Philip

[lindstroem]

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)

[Dennis A]

Originally Posted by gtg258f

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