Test cases

[philipcardiff]

@iBatistic : this seems like a good place to discuss the test cases.

Verification case

I looked at implementing the analytical solution from Timoshenko and Goodier for the stress around a spherical cavity in uniaxial tension; it is straightforward for the stresses along a line, but one terms is not clear to me for the 3-D solution. I will try figure it out but it looks like some information is missing or there is an error. So, using an MMS solution might be best. In which case, we can select a 1x1x1 cube and impose a known nonlinear solution. I will see if I have one already that we can reuse.

Other cases

• ellipticPlate: already set up and easy to run
• narrowTmember: already set up and easy to run
• ideal ventricle: this is set up and the results seem reasonable but I have not check the answers are fully correct. We can do a mesh study and check the predictions converge
• ...

[iBatistic]

@philipcardiff, the solution in Timoshenko's book is not clear to me either. Maybe it could be found in another book, although I doubt it would be more detailed/accurate if it's "wrong" here.

I'm curious, how did you make the mesh for the ideal ventricle case?

[philipcardiff]

I found the spherical cavity solution described in https://apps.dtic.mil/sti/tr/pdf/ADA185392.pdf. I have gone through the equations (they give displacement and stress). There are quite a few terms, but I think I understand how every term is calculated. I will see if I can implement it.

Yep, I have the ideal ventricle case; it is here: https://github.com/solids4foam/cardiacFoam/tree/main/tutorials/LandEtAl2015/problem2/implicitCellCentredSnes However, I plan to copy this case to solid-benchmarks too.

[philipcardiff]

After reading many articles, I finally found the analytical solution for the stress distribution around a spherical cavity in an infinite elastic medium under uniaxial tension, in R.V. Southwell M.A. F.R.S. & H.J. Gough M.B.E. B.Sc. (1926) VI. On the concentration of stress in the neighbourhood of a small spherical flaw; and on the propagation of fatigue fractures in “Statistically Isotropic” materials, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science: Series 7, 1:1, 71-97, DOI: 10.1080/14786442608633614.

I have implemented these as a function object and a traction boundary condition and you can find the case at https://github.com/solids4foam/solids4foam/tree/feature-petsc-snes/tutorials/solids/linearElasticity/sphericalCavity. From some quick checks, the numerical stresses appear to be approaching the analytical stresses so the analytical stresses are probably correct. Also, if you plug in z = 0, it gives the same formula as in Timoshenko and Goodier (which I used in my thesis). We need to move the case to the solid-benchmarks repo and set up the scripts to check the order of accuracy for the stress.

I also found the stress (and displacement) solution in the Bower book (see problem 5.4.8 at http://solidmechanics.org/Text/Chapter5_4/Chapter5_4.php), however, there seems to be errors in these formulae, as the zz stress does not reduce to the formula in Timoshenko when z = 0. You will see that I have implemented these formulae for the displacements in the function object but they are obviously wrong and do not agree with the numerical ones (although I could have made a mistake). In any case, for the paper it will be enough to examine the stress errors.

By the way, one of the seminal references for this problem seems to be Goodier, J. N., "Concentration of Stress Around Spherical and Cylindrical Inclusions and Flaws," J. Appl. Mech., Vol. 55, 1933, p. A39., but unfortunately, I could not find this anywhere online for free; I have asked for it via our library inter-library loan service, but in any case, we already have the stress solution now and I suspect Goodier does not give the displacement solution.

[iBatistic]

Wow, nice! I'm glad you found solutions to this problem. It's interesting how Southwell and Gough used the sigma symbol for Poisson's ratio in their paper, I've never seen that before.

I hope you didn’t spend too much time on Bower's solution. I had a similar experience with another analytical solution once; I only got to the correct equation after tracking down the original work. It seems like you went through something similar.

We need to move the case to the solid-benchmarks repo and set up the scripts to check the order of accuracy for the stress.

If you want, I can take care of that once I return from vacation. Since you've already added the calculations for the L1, L2, and Linf norms in a function object, all we need is a script to run everything and a gnuplot script to visualize

[philipcardiff]

I am surprised at how difficult it is to find this solution for what would appear to be a classic problem.

I am not sure where the error in the Bower solution comes from, although in theory, I could try to go through the derivation he described (tedious). However, in the meantime, I was able to access the original Goodier paper from 1933; here are the displacement and stress equations it gives for a spherical cavity: where A, B and C are defined as

These look straightforward to implement; I will implement them and see if they give the same stresses as Southwell's solution, and hopefully, the displacement fields are correct (as it would be good if we could show the displacement order of accuracy).

By the way, Section 4 of the paper is starting to take shape. Running the cases (or creating the scripts to do this) and generating the plots is the main thing to do.

[philipcardiff]

It seems that we are only getting 1st order accuracy for the displacements:

Either something is not quite right with the analytical solution (or error calculation) or the changes to the gradient calculation at the boundary are causing an issue. It will be useful to also run this case with the segregated solver (standard and "uns" version).

[philipcardiff]

The stress errors are worse (5 meshes):

The max error ($L_\infty$) seems to be zero-order accurate, i.e. not converging.

Here are the displacement errors with 5 meshes (additional finer mesh than above):

Something must be wrong...

Possibly, the analytical solution is slightly different than the numerical one, or the changes to the gradient calculations have introduced an error. Let's see what results the segregated solver gives.

[philipcardiff]

I tried linGeomTotalDispSolid in OFv2012 (Seg) and FE41 (Seg FE41) and unsLinGeomSolid in OFv2012 (Seg Uns) and they show similar or worse behaviour than the JFNK solver (although L2 for Seg FE41 does continue to converge for the finest mesh): Seg Uns in FE41 diverged for the third mesh so I did not include the results.

The max ($L_\infty$) and average ($L_2$) errors should be second order for the displacement so something is not right here.

Here are the displacements errors (JFNK):

and here is a slice through the cavity to the far field top corner (n = (1 -1 0))

I guess it is a bit surprising that there are some large errors near the far-field boundaries; it might suggest some issues there.

At this point, I will go back and run the plateHole verification case to check that the results from the https://doi.org/10.21278/TOF.42301 paper can be reproduced. After that, it may make sense to create a MMS cube case with a structured mesh and fixedValue boundaries before trying an unstructured mesh and traction boundaries.

By the way, I came across a nice book which focusses on code order-of-accuracy verification: https://www.taylorfrancis.com/books/mono/10.1201/9781420035421/verification-computer-codes-computational-science-engineering-patrick-knupp-kambiz-salari

[philipcardiff]

OK, I am making progress.

I went back to FE41 with Seg Uns on the plateHole case and am getting the following results:

#Mesh OrderOfAccuracyDisp(L2,Linf) OrderOfAccuracyStressXX(L2,Linf)
1 1.85562 1.86814 1.86521 1.52826
2 1.97409 1.97392 1.91893 1.74471
3 2.00422 2.00241 1.90945 1.81798
4 2.01881 2.01785 1.86501 1.08996

So, displacements are as expected, and stress is mostly as expected, although, in the paper, Zejko also showed 2nd order for stress LInf. It might be related to the skewCorrected scheme or maybe the solution tolerances need to be tighter.

Here are the XX stress errors on mesh 3 and 4:

where the max error can be seen to occur at the internal skewed cells, whereas we might have expected it to be at the boundary.

By the way, the plateHole analytical solution in solids4foam seems to be for plane strain, whereas in the paper, they used plane stress.

[iBatistic]

I have a question. When we calculate the analytical solution at the cell center and compare it with the numerical solution at the cell center, we're essentially comparing the analytical value at a single point with a numerical solution that represents the mean value over the entire cell. Wouldn't it be more accurate to integrate the analytical solution over the cell, average it, and then compare that averaged value with the numerical solution? This approach seems more consistent since both would represent an average over the cell, rather than comparing a point value with a mean. I understand this might not directly relate to your specific problem, but it seems to me that this method would be more precise, yet nobody seems to use it.

I once made a mistake by using the wrong slope to calculate first- and second-order analytical convergence. You probably didn't make the same mistake, but it might be worth double-checking just in case. The max error convergence seems a bit strange, and I suspect it could be due to a these few cells located far away from the cavity.

At this point, I will go back and run the plateHole verification case to check that the results So, displacements are as expected, and stress is mostly as expected, although, in the paper, Zejko also showed 2nd order for stress LInf. It might be related to the skewCorrected scheme or maybe the solution tolerances need to be tighter.

Hmm... spherical cavity have polyhedral cells, maybe the problem is somehow related to that.

By the way, I came across a nice book which focusses on code order-of-accuracy verification: https://www.taylorfrancis.com/books/mono/10.1201/9781420035421/verification-computer-codes-computational-science-engineering-patrick-knupp-kambiz-salari

thank you, it's great!

[iBatistic]

After that, it may make sense to create a MMS cube case with a structured mesh and fixedValue boundaries before trying an unstructured mesh and traction boundaries.

Do you have any ideas on how to implement the source? It can be easily hard-coded, but it might be more convenient to read it from a text file. A Python script could be used to generate the file. Just an idea... I did the same with the high-order code. Since plateHole is ok I guess you will not go to MMS case, but anyway, I was curious about this.

p.s. I'm returning from vacation this weekend. I was wondering, if I want to add something to the article, should I submit a pull request?

[philipcardiff]

I have a question. When we calculate the analytical solution at the cell center and compare it with the numerical solution at the cell center, we're essentially comparing the analytical value at a single point with a numerical solution that represents the mean value over the entire cell. Wouldn't it be more accurate to integrate the analytical solution over the cell, average it, and then compare that averaged value with the numerical solution? This approach seems more consistent since both would represent an average over the cell, rather than comparing a point value with a mean. I understand this might not directly relate to your specific problem, but it seems to me that this method would be more precise, yet nobody seems to use it.

In my understanding, this will not make a difference in our case as our flavour of second-order finite volume method can be equivalently derived in terms of mean cell values or cell-centre point values. In fact, the cell-centre point values is actually the approach we use in OpenFOAM when deriving the discretisation, e.g. snGrad uses the values at the cell-centre points not the cell mean value. And the volume integral discretisation uses the mid-point rule, which says the value at the centroid is exactly equal to the mean of the cell. So it should not matter. Having said that, currently, I am calculating the L2 norm as the numerical average for all cells, but it would probably make a bit more sense to take the volume-weighted average so that small cells do not overly dominate the average. However, this will not affect the maximum error (L infinity) or the order of accuracy for either.

I once made a mistake by using the wrong slope to calculate first- and second-order analytical convergence. You probably didn't make the same mistake, but it might be worth double-checking just in case. The max error convergence seems a bit strange, and I suspect it could be due to a these few cells located far away from the cavity.

I have also made these mistakes but I think it is correct here. For plateHole, I am using the following for displacement:

# Spacing 1st 2nd
6 6.125 0.005
12 12.5 0.02
24 25 0.08
48 50 0.32

At this point, I will go back and run the plateHole verification case to check that the results So, displacements are as expected, and stress is mostly as expected, although, in the paper, Zejko also showed 2nd order for stress LInf. It might be related to the skewCorrected scheme or maybe the solution tolerances need to be tighter.

Hmm... spherical cavity have polyhedral cells, maybe the problem is somehow related to that.

Yep, possibly. Unstructured meshes will definitely be more difficult than nice smoothly varying hexahedra.

[philipcardiff]

After that, it may make sense to create a MMS cube case with a structured mesh and fixedValue boundaries before trying an unstructured mesh and traction boundaries.

Do you have any ideas on how to implement the source? It can be easily hard-coded, but it might be more convenient to read it from a text file. A Python script could be used to generate the file. Just an idea... I did the same with the high-order code. Since plateHole is ok I guess you will not go to MMS case, but anyway, I was curious about this.

Federico hardcoded it into a copy of the solid model, e.g. see https://github.com/federico081997/solidBenchmarks/blob/main/smallStrainPaper/linearElasticity/cubeMMS/cubeMMS/solidModels/linGeomSolidMMS/linGeomSolidMMS.C where he uses this form:

In this case, we selected a manufactured solution that was zero on all the boundaries on a 1x1x1 cube with the origin at (0,0,0); this makes the boundary conditions easy (zero displacement). You will see that we also created a function object to calculate the errors (ideally there should be just one class, to avoid code duplication): https://github.com/federico081997/solidBenchmarks/blob/main/smallStrainPaper/linearElasticity/cubeMMS/vertexCentredLinearGeometry/structuredTriangles/base/system/controlDict

A more general solution is to use fvOptions but it is annoying that this is different between versions. Maybe we can add something like

#ifdef OPENFOAM_COM
    DEqn += fvOptions(D());
#endif

and then later we can get it working with FE, if needed.

p.s. I'm returning from vacation this weekend. I was wondering, if I want to add something to the article, should I submit a pull request?

[iBatistic]

In my understanding, this will not make a difference in our case as our flavour of second-order finite volume method can be equivalently derived in terms of mean cell values or cell-centre point values. In fact, the cell-centre point values is actually the approach we use in OpenFOAM when deriving the discretisation, e.g. snGrad uses the values at the cell-centre points not the cell mean value. And the volume integral discretisation uses the mid-point rule, which says the value at the centroid is exactly equal to the mean of the cell. So it should not matter.

Thanks for the explanation. I hadn't considered it from that point of view, I'm glad you cleared it up for me.

For plateHole, I am using the following for displacement:

# Spacing 1st 2nd
6 6.125 0.005
12 12.5 0.02
24 25 0.08
48 50 0.32

How do you calculate the average cell spacing in the case of a polyhedral mesh?

[philipcardiff]

How do you calculate the average cell spacing in the case of a polyhedral mesh?

I was using cbrt(V/nCells), where V is the total volume of the entire domain and nCells is the total number of cells. In the Gmsh script, I specify the cell size at the cavity and then a larger cell size at the far boundaries. I half both of these values in successive meshes.

I think this should be OK, but, yep, it is possible that this could distort the apparent order of accuracy.

You will have seen I emailed Zeljko about the plateHole case from the 2018 paper. Once I reproduce the behaviour there, then I will incrementally move back towards this case and/or maybe a MMS case.

[iBatistic]

I asked because sometimes the mesher does not strictly enforce the prescribed cell size values. I agree that cbrt(V/nCells) should be ok.

I double-checked the code for the stress analytical solution of the spherical cavity; found nothing wrong.

You will have seen I emailed Zeljko about the plateHole case from the 2018 paper. Once I reproduce the behaviour there, then I will incrementally move back towards this case and/or maybe a MMS case.

Understand;

[iBatistic]

I recall seeing a plateHole case where analytical displacement was prescribed instead of traction. Perhaps stress cannot be second order without this — just a thought. You might know better.

[philipcardiff]

Possibly. If the exact solution is prescribed at the boundary, then I guess the behaviour should be the same, assuming displacement and traction conditions are discretised to the same order of accuracy.

Zeljko is checking the extend bazaar code, which should be 2nd order for displacement and stress. Separately, I will use the plateHole (or MMS) solution to figure out which component of the discretisation is causing the issue, i.e. I will hard-code the exact displacement and see if the stress is 2nd order, etc.

I will keep you posted.

[iBatistic]

I agree. I didn't express myself well. I meant that we have finite dimensions of the domain, which might affect the stress convergence rate. I guess it is negligible, but maybe it is worth checking

I'm back to work, please let me know what I can do to help.

[philipcardiff]

Ah, ok. In my understanding, using the analytical displacement or analytical stress conditions should both eliminate the finite dimensions issues. I suspect that both approaches will give different absolute accuracy but the order of accuracy should be unaffected, assuming both displacement and traction conditions are implemented using the same order of accuracy.

In terms of progress on the paper, if you could help by setting up the cases and scripts (AllrunMeshStudy, gnuplot scripts, etc.) for any case in the paper that does not have it, that would be great. Feel free to add the graphs/data to the paper as you go; we can run all the cases again on the same platform for consistency before submission.

By the way, I set up a MMS case at https://github.com/solids4foam/solid-benchmarks/tree/main/linearElasticity/manufacturedSolution. Using a cube domain and a structured hexahedral mesh, the results are as expeced: 2nd order of average and max displacement errors, 2nd order for average stress errors, and 1st order for max stress errors.

Displacement errors:

Stress errors:

Above I used linGeomTotalDispSolid. Using Zeljko's "uns" approach should give 2nd order for the max errors too; I will check this later, but 1st order for max stress errors is not bad.

Next I will try an unstructured mesh. If that works, then we could use this case as the verification case in the paper. If it does not work as expected (lower order of accuracy) then I will perform checks to find the problem.

[philipcardiff]

Good news and bad news:

Good news: the polyhedral mesh shows the same behaviour on the MMS case as on the sphericalCavity case. Bad news: the displacement and stress errors approach zero order of accuracy as the mesh is refined.

Displacement errors:

Stress errors:

I will hard code various parts of the solution and see where the errors come from...

[iBatistic]

In terms of progress on the paper, if you could help by setting up the cases and scripts (AllrunMeshStudy, gnuplot scripts, etc.) for any case in the paper that does not have it, that would be great. Feel free to add the graphs/data to the paper as you go; we can run all the cases again on the same platform for consistency before submission.

Sure, I will start immediately!

Using a cube domain and a structured hexahedral mesh, the results are as expeced: 2nd order of average and max displacement errors, 2nd order for average stress errors, and 1st order for max stress errors. Good news: the polyhedral mesh shows the same behaviour on the MMS case as on the sphericalCavity case. Bad news: the displacement and stress errors approach zero order of accuracy as the mesh is refined.

Crazy. I guess we can say that we have mesh induced error. Hmm, even the 'skewCorrected` scheme does not help... My first guess is that the cause is LS gradient computation. I hope you find the cause soon.

[philipcardiff]

More progress...

I hard-coded and the exact (MMS) cell-centre displacements and then looked at the cell-centre stress errors. In addition, I hard-coded the stress values in cells adjacent to the boundaries to remove any boundary condition effects.

In this case (unstructured polyhedra), the max and average stress errors (interior cells) are close to 1st order. Looking again at the theory (see https://www.sciencedirect.com/science/article/pii/S0378475422003743), this seems to be the expected behaviour. Here is a quote from the Syrakos et al. paper: It is also noted that all of the presented gradient schemes are formally of first-order accuracy (in practice they are second-order accurate on smooth grids). In their study, they showed that the gradient order of accuracy is close to second on smooth quadrilateral grids, e.g. and drops to 1st order (or lower) on random triangular grids:

They discussed different linear approaches for the gradient calculation, but all are 1st order at best on random grids. They even note that it can be less than 1st if the grid refinement reduces the cell quality.

This does not explain the zero-order results I am getting but I suspect this is a problem at the boundary that can be fixed. For the paper, it may be easiest to proceed with a smooth grid for the verification study and simply note in the text that the stress order will drop for random grids. In any case, 1st order gradients should be sufficient to product 2nd order displacement errors.

By the way, I implemented a naive quadratic least squares gradient calculation using Eigen (similar to the higher-order approach) and this seems to approach 2nd order accuracy for the average and max stress errors (I have only checked the "internal" cells). It would be interesting to see if this would work with the overall solution procedure (i.e. use it in the momentum equation) but this may be too much work for this paper and it may detract from the main aim.

FYI: @ztukovic, in case you are interested.

[philipcardiff]

Here is the order of accuracy on the MMS case when I used a structured polyhedral mesh:

The mesh is created by creating a structured tetrahedral mesh in Gmsh and converting this to polyhedra using polyDualMesh. Here is a section through mesh3:

These results seem reasonable and are probably fine for the paper. Optionally, we could also add the structured hexahedra results. The "uns" solver will likely do better in terms of accuracy and order of accuracy, so we could make a JFNK version of it too.