CFD Online Discussion Forums - [solidMechanics] Support thread for "Solid Mechanics Solvers added to OpenFOAM Extend"

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM CC Toolkits for Fluid-Structure Interaction (https://www.cfd-online.com/Forums/openfoam-cc-toolkits-fluid-structure-interaction/)

- - [solidMechanics] Support thread for "Solid Mechanics Solvers added to OpenFOAM Extend" (https://www.cfd-online.com/Forums/openfoam-cc-toolkits-fluid-structure-interaction/126706-support-thread-solid-mechanics-solvers-added-openfoam-extend.html)

Quote:

Originally Posted by codder (Post 484047)

@ Dr. Cardiff -

I believe there exists a bug in the arbitraryCrack solid model. I have tested 2D trigonal meshes with the elasticAcpSolidFoam solver using meshes generated in both gmsh and ICEM (having initially suspected the error to be associated with a mesh import utility).

But I am now thinking the issue is:

$FOAM_SRC/solidModels/arbitraryCrack/faceCracker/detachFaceCracker.C

As part of this issue, I see under paraview warp-by-vector an unphysical separation of faces not a member of (but adjacent to) the crack patch; conjointly, elements with no opportunity to detach nodes, because no traction will develop on attached faces.

By thumbnail:

#1 - uncracked mesh - crackable boundary patch specified along the five triangular notches
#2 - cracked mesh with labeled displaced problem face (time = 54)
#3 - crack patch in slightly 3-D view, similarly labelled (time = 54)

For example, compilation using "detachFaceCracker_orig.C" worsens the issue.

Best, Eric

Hi Eric,

I was aware of this issue with certain meshes, but I finally got around to looking at it and hopefully it is now solved with the attached detachFaceCracker.C.

Try is out and let me know.

Best regards,
Philip

Stress & Displacement Errors

Hi,

A thin beam has been modeled by plane strain problem. When the number of cells along the thickness increases there is some error increasing in the simulation but by increasing number of cells along the length there is decreasing in the error, as shown in attached figures. Why does increasing cell numbers along the thickness cause an increasing in the error, please?

Thank you.

Quote:

Originally Posted by Brayanashel (Post 489143)

Hi,

Using the finite volume method in its current linear form, large aspect ratio cells act stiffer in bending. So if you have long skinny cells along the length of the beam then it will act more stiff in bending than it should be.
This behaviour can also occur in finite elements.

I think someone has noted this earlier in this thread.

Philip

Quote:

Originally Posted by bigphil (Post 488716)

Hi Eric,

Try is out and let me know.

Hi Phil -

To my eyes, the problem is gone.

A big thanks for your bug fix. A second big thanks to the continued generosity of the UCD lab for releasing the code.

Best, Eric

Quote:

Originally Posted by misklach (Post 488205)

Thank you very much Philip for your support.

Ok I fixed the issue with the tolerance (the tolerance for the linear solver was bigger than the convergenceTolerance ...shame on me!) ...but apparently it doesn't affect very much the result.

For a time step of 0.005 at time 10 the displacement is still the same. A zoom-in at the middle point shows that the difference is really negligible.
DU convergenceTolerances are 1e-7(blue), 1e-9(red), 1e-12(green)
Attachment 30377

Attachment 30378

Decreasing the timestep still makes the solution diverge at some point
https://drive.google.com/file/d/0B1d...it?usp=sharing
deltaT = 5e-3(blu), 3.125e-4(red), 7.8125e-5(green)
DU convergenceTolerance = 1e-9

Hi Giampaolo,

You could try this case with a total displacement solver such as elasticNonLinTLSolidFoam instead of the incremental displacement solver and see how it compares.

There are only a few changes to the case required:
change all BCs with "nonLinear updatedLagrangian" to "nonLinear totalLagrangian"; and also change all "DU" to "U" in fvSchemes/fvSolution.

Best regards,
Philip

Hi Philip,

unfortunately there are no improvements using the total lagrangian solver. For the same time-step TL and UL give different behaviour.
All the simulations use backward scheme and DU tolerance = 1e-09.
https://drive.google.com/file/d/0B1d...it?usp=sharing
https://drive.google.com/file/d/0B1d...it?usp=sharing

At least all the total lagrange simulations converge (with different modes) towards the steady state solution ...but now where the dissipation/dumping comes from?
https://drive.google.com/file/d/0B1d...it?usp=sharing

Best regards
Giamp

Quote:

Originally Posted by misklach (Post 491025)

HI Giamp,

for the TL simulations using the backward scheme, did you modify the temporal term in the solver to be rho*d2dt2(U) instead of d2dt2(rho, U)?

Philip

Hi Philip,

also with the temporal term rho*fvm::d2dt2(U) the solution is damped :(
https://drive.google.com/file/d/0B1d...it?usp=sharing

decreasing the time step, even if dumped, at least the solution seems to not change anymore
https://drive.google.com/file/d/0B1d...it?usp=sharing

but unfortunately at some some it crashes
these are the last lines from the log of DeltaT=3.125e-4

Code:

Time 1.87156, Corrector 9992, Solving for U using DICPCG, residual = 3.70566e-05, relative residual = 2.45566, inner iterations 9

        Time 1.87156, Corrector 9993, Solving for U using DICPCG, residual = 3.51757e-05, relative residual = 0.244511, inner iterations 4

        Time 1.87156, Corrector 9994, Solving for U using DICPCG, residual = 1.73622e-05, relative residual = 1.30172, inner iterations 18

        Time 1.87156, Corrector 9995, Solving for U using DICPCG, residual = 2.47067e-05, relative residual = 0.183089, inner iterations 1

        Time 1.87156, Corrector 9996, Solving for U using DICPCG, residual = 4.02368e-05, relative residual = 0.197761, inner iterations 1

        Time 1.87156, Corrector 9997, Solving for U using DICPCG, residual = 3.76838e-05, relative residual = 0.389101, inner iterations 1

        Time 1.87156, Corrector 9998, Solving for U using DICPCG, residual = 2.09545e-05, relative residual = 0.672903, inner iterations 2

        Time 1.87156, Corrector 9999, Solving for U using DICPCG, residual = 1.19139e-05, relative residual = 0.731561, inner iterations 6



Time 1.87156, Solving for U, Initial residual = 0.0049075, Final residual = 1.19139e-05, Relative residual = 0.731561, No outer iterations 10000

ExecutionTime = 2376.78 s  ClockTime = 2398 s

ExecutionTime = 2376.78 s





Time: 1.87187



        Time 1.87187, Corrector 0, Solving for U using DICPCG, residual = 0.00534617, relative residual = 1, inner iterations 1

        Time 1.87187, Corrector 1, Solving for U using DICPCG, residual = 1.96508e-05, relative residual = 1.01578, inner iterations 14

        Time 1.87187, Corrector 2, Solving for U using DICPCG, residual = 2.91777e-05, relative residual = 5.95844, inner iterations 5

        Time 1.87187, Corrector 3, Solving for U using DICPCG, residual = 2.74333e-05, relative residual = 1.4765, inner iterations 2

        Time 1.87187, Corrector 4, Solving for U using DICPCG, residual = 2.01264e-05, relative residual = 0.445946, inner iterations 1

        Time 1.87187, Corrector 5, Solving for U using DICPCG, residual = 1.45802e-05, relative residual = 4.99339, inner iterations 4

        Time 1.87187, Corrector 6, Solving for U using DICPCG, residual = 2.34268e-05, relative residual = 0.336645, inner iterations 1

        Time 1.87187, Corrector 7, Solving for U using DICPCG, residual = 2.84046e-05, relative residual = 0.698792, inner iterations 2

        Time 1.87187, Corrector 8, Solving for U using DICPCG, residual = 3.61341e-05, relative residual = 2.96854, inner iterations 5

        Time 1.87187, Corrector 9, Solving for U using DICPCG, residual = 3.78643e-05, relative residual = 2.09039, inner iterations 2

        Time 1.87187, Corrector 10, Solving for U using DICPCG, residual = 3.68786e-05, relative residual = 4.16278, inner iterations 7

        Time 1.87187, Corrector 11, Solving for U using DICPCG, residual = 3.29223e-05, relative residual = 1.70304, inner iterations 4

        Time 1.87187, Corrector 12, Solving for U using DICPCG, residual = 3.08959e-05, relative residual = 9.57103, inner iterations 8

        Time 1.87187, Corrector 13, Solving for U using DICPCG, residual = 4.94034e-05, relative residual = 0.192556, inner iterations 1

        Time 1.87187, Corrector 14, Solving for U using DICPCG, residual = 9.53101e-05, relative residual = 0.190417, inner iterations 1

        Time 1.87187, Corrector 15, Solving for U using DICPCG, residual = 0.000259233, relative residual = 0.18017, inner iterations 1

        Time 1.87187, Corrector 16, Solving for U using DICPCG, residual = 0.00113276, relative residual = 1.74945, inner iterations 1

        Time 1.87187, Corrector 17, Solving for U using DICPCG, residual = 0.0123064, relative residual = 1.9589, inner iterations 1

        Time 1.87187, Corrector 18, Solving for U using DICPCG, residual = 0.377465, relative residual = 1.48288, inner iterations 1

        Time 1.87187, Corrector 19, Solving for U using DICPCG, residual = 0.995577, relative residual = 1.00002, inner iterations 1

        Time 1.87187, Corrector 20, Solving for U using DICPCG, residual = 1, relative residual = 1, inner iterations 1

--> FOAM Warning :

    From function eigenValues(const tensor&)

    in file primitives/Tensor/tensor/tensor.C at line 170

    complex eigenvalues detected for tensor: (-5.24059e+18 1.52383e+20 6.46432e-26 378752 1 -4.87694e-11 653.646 -18991.2 1)

--> FOAM Warning :

    From function eigenValues(const tensor&)

    in file primitives/Tensor/tensor/tensor.C at line 170

    complex eigenvalues detected for tensor: (-5.24059e+18 1.52383e+20 6.46432e-26 378752 1 -4.87694e-11 653.646 -18991.2 1)

Here there is the case that I'm running if you want to check it...
https://drive.google.com/file/d/0B1d...it?usp=sharing

Best regards,
Giamp

Hi Giamp,

Hmnnn for transient analysis the beam should never reach steady-state, can you check if a small bug-fix in TL solver makes a difference:
in writeFields.H, change this:

Code:

rho = rho/J;

to this:

Code:

//rho = rho/J;

Also, it would be interesting to check if the small stain solver elasticSolidFoam has same issues.

Philip

good! with the bug-fix, the beam keeps oscillating
https://drive.google.com/file/d/0B1d...it?usp=sharing

but still, decreasing the time-step, at some point the solutions are different. ElasticSolidFoam keeps oscillating but crashes when decreasing the time-step
https://drive.google.com/file/d/0B1d...it?usp=sharing

Hi Dr. Cardiff

Have you ever tried the FSI solver using velocityLaplacian mesh motion solver in parallel?
I am facing an error:

Code:

Constructing processor meshes



Processor 0

    Number of cells = 10470

    Number of faces shared with processor 1 = 171

    Number of processor patches = 1

    Number of processor faces = 171

    Number of boundary faces = 21405



Processor 1

    Number of cells = 10470

    Number of faces shared with processor 0 = 171

    Number of processor patches = 1

    Number of processor faces = 171

    Number of boundary faces = 21395



Number of processor faces = 171

Max number of processor patches = 1

Max number of faces between processors = 171



Processor 0: field transfer





--> FOAM FATAL ERROR: 

size of field = 21744 is not the same as the size of mesh = 21578



    From function DimensionedField<Type, GeoMesh>::DimensionedField(const IOobject& io,const Mesh& mesh, const dimensionSet& dims, const Field<Type>& field)

    in file /opt/OpenFOAM/OpenFOAM-3.0-ext/src/OpenFOAM/lnInclude/DimensionedField.C at line 71.



FOAM aborting



Aborted (core dumped)

This problem does not happen with FEM solvers.

Best wishes

Hi Philip,
I'm trying to modify the solver elasticNonULSolidFoam using an incompressible Mooney-Rivlin hyperelastic model.

From the strain energy function

W = c1*(tr(C) - 3) + c2*(0.5*((tr(C))^2 - tr(C^2)) -3)

I derived the second Piola-Kirchhoff stress tensor in its incremental form

DSigma = 2*c1*I + 2*c2*tr(C)*I - 2*c2*C - p*inv(C)

where C is the right Cauchy-Green deformation tensor (C=F.T() & F), and c1+c2=mu/2

Following
http://powerlab.fsb.hr/ped/kturbo/Op...tressPaper.pdf
and
http://powerlab.fsb.hr/ped/kturbo/op...18-06-2007.pdf
I came up with the following Updated lagrangian momentum equation (neglecting the pressure term)

rho*fvm::d2dt2(DU)
==
- fvm::laplacian(2*c2, DU, "laplacian(DDU,DU)")
+ fvc::div
(
4*c2*tr(gradDU)*I
+ 2*c2*(gradDU.T()),
"div(sigma)"
)
+ fvc::div
(
2*c2*(I*tr(gradDU.T() & gradDU))
- 2*c2*(gradDU.T() & gradDU)
+ ((sigma + DSigma) & gradDU),
"div(sigma)"
)

Basically in elasticNonLinULSolidFoam I changed only the momentum equation and DSigma (neglecting the pressure term for now). I calculated C as C = (2*DEpsilon + I) and used c1=0.4375*mu and c2=0.0625*mu.

I tried to run a steady state simulation with this solver with a case which works fine with elasticNonLinULSolidFoam, but unfortunately the simulation crashes

Code:

        Time 0.0001, Corrector 61, Solving for DU using GAMG, res = 0.786476, rel res = 0.117896, inner iters 1

        Time 0.0001, Corrector 62, Solving for DU using GAMG, res = 0.783442, rel res = 0.116527, inner iters 1

        Time 0.0001, Corrector 63, Solving for DU using GAMG, res = 0.780196, rel res = 0.115416, inner iters 1

        Time 0.0001, Corrector 64, Solving for DU using GAMG, res = 0.776698, rel res = 0.116016, inner iters 1

        Time 0.0001, Corrector 65, Solving for DU using GAMG, res = 0.77302, rel res = 0.156058, inner iters 1

        Time 0.0001, Corrector 66, Solving for DU using GAMG, res = 0.799312, rel res = 0.628151, inner iters 1

        Time 0.0001, Corrector 67, Solving for DU using GAMG, res = 0.998106, rel res = 0.998607, inner iters 2

        Time 0.0001, Corrector 68, Solving for DU using GAMG, res = 1, rel res = 0.999996, inner iters 2

        Time 0.0001, Corrector 69, Solving for DU using GAMG, res = 1, rel res = 1, inner iters 2

        Time 0.0001, Corrector 70, Solving for DU using GAMG, res = 1, rel res = 1, inner iters 2

Floating point exception (core dumped)

I don't know if the problem comes from the equations (maybe because I neglected the pressure term?) or from the discretization. Apparently I don't know in my case how to do the trick with the implicit/explicit manipulation of the laplacian(DU).

Any suggestion to make it work?

About the pressure term, since that should be an arbitrary pressure, does assuming p=0 make sense? ...or the problem might come from this assumption?

Thank you very much!

Hi Giampaolo,

See this paper by Bijelonja et al, they develop the mathematical model for an incompressible Mooney-Rivlin material using a large strain total Lagrangian approach. They use the SIMPLE algorithm to couple the pressure term, just like in fluids.

So essentially merge elasticNonLinTLSolidFoam and simpleFoam and change calculation of stress in terms of Mooney-Rivlin.
Also traction boundary gradient method must be changed.

Let me know how it goes, I previously looked at creating a solid solver for incompressible materials using SIMPLE/PISO coupling and I had problems with traction boundary conditions, but it is possible as the authors above have achieved it using cell-centred Finite Volume method.

Good luck,

By the way, as regards the time-step defence of the oscillating membrane problem, I am not sure what the issue is, I suppose a more accurate time discretisation (e.g. Newmark, Crank-Nicholson, etc.) would help but I am not sure.

Philip

MesquiteMotionSolver with icoFsiElasticNonLinULSolidFoam

@Philip
I wanted to remesh the deformed part of the mesh while using icoFsiElasticNonLinULSolidFoam but when i used dynamicTopoFvMesh along with mesquiteMotionSolver I got this error msg!!!

--> FOAM FATAL ERROR:
Problem with mesh motion solver selection

From function icoFsiElasticNonLinULSolidFoam
in file moveFluidMesh.H at line 140.

FOAM aborting

Aborted (core dumped)

and when I forced the moveFluidMesh.H to use the given motion solver, it reported an error msg.:

--> FOAM FATAL ERROR:

request for tetPointVectorField motionU from objectRegistry region0 failed
available objects of type tetPointVectorField are

0
(
)

From function objectRegistry::lookupObject<Type>(const word&) const
in file /home/a/foam/foam-extend-3.0/src/foam/lnInclude/objectRegistryTemplates.C at line 139.

FOAM aborting

can I know where the problem is? and what should I do to solve this issue?

regards,

MA

Hi,

How can a traction and displacement BCs be applied simultaneously in a face, please? For example, having tractionDisplacement and slip types BC in the same time in a face.

Thanks.

Quote:

Originally Posted by Brayanashel (Post 494495)

Hi,

How can a traction and displacement BCs be applied simultaneously in a face, please? For example, having tractionDisplacement and slip types BC in the same time in a face.

Thanks.

If you mean that the normal displacement is specified and the tangential traction is zero then use the fixedDisplacementZeroShear boundary condition.

Philip

Actually there are shear and displacement in a single face, as shown in the left BC in attachment.

Quote:

Originally Posted by Brayanashel (Post 494529)

Actually there are shear and displacement in a single face, as shown in the left BC in attachment.

I am not sure I understand your schematic drawing;

there are 'frictionless' rollers on the left which implies there is zero shear?
And there is a fixed point at the centre of the left boundary.

If this is the case, then you would split the left boundary so as there is a thin fixedDisplacement boundary for the fixed part and fixedDisplacementZeroShear for the other parts.

Philip

Please look at the paper in attachment. I would like to know how beams D and E in Fig. 2 on page 2 can be modeled by stressedFoam. Applying the left dual BCs is my question.

Thanks.

Quote:

Originally Posted by Brayanashel (Post 494537)

Please look at the paper in attachment. I would like to know how beams D and E in Fig. 2 on page 2 can be modeled by stressedFoam. Applying the left dual BCs is my question.

Thanks.

attachment?