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)

I know there is an older version of this Solid Mechanics addon for OpenFOAM 2.0-2.3, but has anyone ported it to the ESI version of OpenFOAM (OpenFOAM-v1706 and such - the ever changing names becomes confusing)?

Quote:

Originally Posted by bigphil (Post 651072)

Hi Bartosz,

The problem is with the globalFaceZones: after a topological change, these get messed up and it is not trivially to fix them; consequently, as you have seen the GGI will complain. In your case, is there a topo change at the patch/faceZone? If so, have you tried a case when the topo change is away from the patch? Just an idea.

Similar issue discussed here: https://sourceforge.net/p/openfoam-e...ndrelease/280/

Philip

Dear Prof. Cardiff

I faced the same problem trying to use a dynamic motion solver coupled with adaptive mesh refinement with fsiFoam (dynamicMotionSolverRefineFvMesh). Is it possible to fix the mess after the topological change? Or is there another way to implement this kind of solver in fsiFoam?

Best regards
Fábio Fonseca

Hey All,
I am new to OpenFoam and found the solidMechanics solvers to be very helpful.
So now to my Problem:D:
I want to use the Gent elasticity model with the Terzaghi stress given by:

Code:

1/det(F)*[I_lim/(I_lim-(Tr(F*F.T)-3))*F.T*F-I]

The model is nonlinear and I am not sure if there is an easy way to use this
model with the solidMechanics solvers. All the solvers and most of the BCs
seem to need the Lamé coefficients. Has someone already done something
similar or knows how this could be done?

Quote:

Originally Posted by Tshiro (Post 687237)

Hey All,
I am new to OpenFoam and found the solidMechanics solvers to be very helpful.
So now to my Problem:D:
I want to use the Gent elasticity model with the Terzaghi stress given by:

Code:

1/det(F)*[I_lim/(I_lim-(Tr(F*F.T)-3))*F.T*F-I]

Hi Tshiro,

I have been working on extending the solidMechanics code to make it easier to implement new constitutive laws such as this.

I have some other hyperelastic laws implemented and you are welcome to try it out, just PM me your email address and I will send you a link.

Philip

Hi, we are trying to run a FSI simulation (OF-ext 4) using the PronyViscoelastic rheology model. We are clueless on how to define the parameters (k, kDimensions, etc) of this model and will thank any help.

Could any of you give us an example?

Or, for example, how the model g(t)=[1-g1*(1-exp(-t/tau1))-g2*(1-exp(-t/tau2))] could be implemented???

Thanks a lot!!!

Quote:

Originally Posted by cguti (Post 688315)

Hi cguti,

Separately from the PronySeries law in OF-4.0-EXT, I have also implemented the viscoelastic approach described in "Simo and Hughes, 2000, Computational Inelasticity", which uses a relaxation function just like you have shown.

PM me your email and I will share it with you.

Philip

elasticOrthoSolidFoam

Hello All,

Currently I'm trying to familiarize myself with OpenFOAM, particularly the elasticOrthoSolidFoam solver. (Foam Extend 4.0)
First I just wanted to understand the basics so I took a simple cantilever, and tried to see how the simulation goes for normal and shear (bending) force.

After I've read through this thread - if I understand correctly - I see that this shear case might be a little problematic with this solver.
I thought about stress singularity so I've made a "base plate" with a fillet but I still failed to get good results. Setting nCorrectors to 20000 in fvSolutions did not help.

What should I do in order to get close to the analytical results? With the normal force case I get perfect results.

Here is a link to the case folder, if needed:
https://www.dropbox.com/sh/m9roulbqe...grmz2VQNa?dl=0

Thanks,
Bence

Quote:

Originally Posted by szbence (Post 689241)

Hi Bence,

Yes, you are right that the segregated algorithm employed by the solver is not so efficient in a bending cantilever case. However, the stress singularity does not cause a problem (for convergence or efficiency).
From checkMesh, I can see that there are a number of severely non-orthogonal faces, these are not good for convergence (though the solver does seem to be converging).

I suggest you try the case without the base plate and increase the nCorrectors further until convergence is achieved at least for a residual of 1e-6: this should produce predictions close to the analytical results; if it slightly underestimates then you will need to tighten this tolerance.

Separately, if you PM me your email, I can share a case using the solid block coupled solver, which will be much faster for this case (but orthotropy is not implemented yet...).

Philip

Quote:

Originally Posted by bigphil (Post 689244)

Dear Philip,

Thank you for your quick response, I will PM you!
As I'm using a virtualbox the 30000 iterations took 2 hours, (without the base plate) but the residuals "look nice" even though it would need some more time to reach the target value. (https://www.dropbox.com/s/nlt3x48kju...uals3.JPG?dl=0) The results are getting closer to the analytical one, maybe I'll try a longer simulation during night, but I believe the results would be all right.

Bence

Quote:

Originally Posted by rege (Post 649932)

Hi all,

I've been trying to implement LES to the "beamInCrossFlow" tutorial for the FSI toolkit for foam-extend 4.0.
My problem is that a heavily oscillating pressure field appears in the fluid field, eventually causing extreme velocities, causing the simulation to crash.

My changes with respect to the original tutorial are as follows:
- Rewritten blockMeshDict to a complete 3D domain for both the fluid and structure.
- Increased the max speed in setInletVelocity to 10, in order to get turbulence.
- Changed the properties of the solid to the following: rho = 8000, E = 4e6, nu = 0.3
- Implemented PISO solver with LES model oneEqEddy, with vanDriest smoothing.
- Decreased FSI tolerance to 1e-5 and increased nOuterCorr to 200.
- ddtScheme: CrankNicolson 0.8
- div(phi,U) Gauss filteredLinear2 1 0
- fvSolution changed to common settings for LES simulations.
- The Courant number has been kept at 0.4-0.5.

The ddtScheme was chosen by trying out all the ddtSchemes without FSI.
The use of other ddtSchemes, including CrankNicolson 0.7 or 0.9, causes the simulations without FSI to blow up, too.
Simulations with CrankNicolson 0.8 run smoothly without FSI however.

In all cases, the problem has been that a field of high pressure at the lower edge of the inlet has merged with the field of high pressure in front of the beam. When these two fields of high pressure merge, they start to "pump" the fluid over the beam, resulting in pressure oscillations for the entire field.
The pressure typically ends up oscillating from +200 Pa to -200 Pa for 80% of the internal field within 0.05 seconds.
I therefore assumed that the inlet was too close to the beam, and increased the length from the inlet by 0.5m. The pressure oscillations did still occur.

I have also tried the following adjustments:
- Max speed 5 m/s.
- Initial internal fields calculated by consistentIcoFluid with FSI, by LES simulation without FSI, and by increasing the velocity from zero with LES and FSI.
- Reduced cell lengths to the half in all directions.
- Adjusted time step to Co = 0.2 and Co = 0.8.

In all these configurations, the pressure oscillations do still occur.
The pressure oscillations do not occur when doing a LES simulation of the undeformed geometry without FSI, and neither if the stiffness (E) of the solid is set to be very high in the FSI simulation (causing it to be stationary).

Has anyone else tried to implement LES to this case?
I would be very grateful for any suggestions.

Best regards,
KR

Hi Kristen,

Have you tried using Euler D2dt2 (D_) scheme? it could be more stable. the CrankNicolson also may be better to be replaced with more stable Euler scheme.

Regards
Ali

Quick question about the hronTurek FSI function object. I'm using the FSI toolkit for foam-extend-4.0, however I've run into an issue with the hronTurek and pointHistory function objects.

Initially they would not compile, giving me errors about the virtual functions being pure. I fixed this by adding in public member functions for "updateMesh" and "movePoints" in the header file, and changing virtual bool execute() to virtual bool execute(const bool forceWrite);

They now compile successfully, however they are only outputting force and position data for the very first time-step. After this they no longer are writing any output.

Have the function object definitions changed in the updated foam-extend version, or perhaps my modifications to fix the compiling issues were in error? :confused:

Hi all,

I hope this is the right place to ask for help. I just want to understand the math theory and mechanics of a very simple solidMechanics solver, solidDisplacementFoam.

I have initial OF coding experience but in fluid FVM. I hope someone here could offer me some basic materials and tutorials for solidDisplacementFoam, and much better if they can explain to me the code structure details corresponding to the full mechanics including thermal stress implementation in the solver.

Many thanks in advance.

Quote:

Originally Posted by cfdopenfoam (Post 699344)

Hi cfdopenfoam,

The solidsDisplacementFoam solver is described in the following articles:

Weller et al. (1998) A tensorial approach to computational continuum mechanics using object-oriented techniques link
Jasak and Weller (2000) Application of the finite volume method and unstructured meshes to linear elasticity link

Also, I have provided descriptions of the background theory in some of my conference presentations (see my ResearchGate), such as:

Introduction to Solid Mechanics with OpenFOAM: Basic Training link
The finite volume method for solid mechanics in OpenFOAM link
Towards the Development of an Extendable Solid Mechanics and Fluid-Solid Interaction Toolbox for OpenFOAM link

I hope it helps,
Philip

Quote:

Originally Posted by bigphil (Post 699357)

Hi cfdopenfoam,

The solidsDisplacementFoam solver is described in the following articles:

Weller et al. (1998) A tensorial approach to computational continuum mechanics using object-oriented techniques link
Jasak and Weller (2000) Application of the finite volume method and unstructured meshes to linear elasticity link

Also, I have provided descriptions of the background theory in some of my conference presentations (see my ResearchGate), such as:

Introduction to Solid Mechanics with OpenFOAM: Basic Training link
The finite volume method for solid mechanics in OpenFOAM link
Towards the Development of an Extendable Solid Mechanics and Fluid-Solid Interaction Toolbox for OpenFOAM link

I hope it helps,
Philip

Thank you, Philip. This is great help. I appreciate. But I still cannot understand the mechanics of the thermal stress included in the solver. Any math/theory background on this?

Hooke's law can be extended to the so-called Duhamel-Neumann form of Hooke's law by including a thermal stress term e.g. books like Mase and Mase, Continuum Mechanics for Engineers discuss linear thermoelasticity.

So the solid stress is a function of temperature and you end up with a source term in the momentum equation that is a function of temperature. For cases where the kinetic energy is small, this temperature equation is decoupled for momentum so you can solve temperature first and then solve momentum (only one-way coupling).

If you PM your email, I can send you some slides I have that might help explain its OpenFOAM implementation.

Philip

Hi all,

just one more question here,

is mesh updated in solidDisplacementFoam (I guess not)? and how?

if I want to implement mesh motion/displacement due to solid mechanics, is there a basic solverFoam in OF and related tutorials?

Thanks in advance.

Hi,

short: No. The basic stress solvers in the Foundation version do not support mesh motion such as plastic movement of the solid. The solvers only handle linear-elastic problems. If you are interested into problems which include plasticity, you should check the foamExtend version. Phil made huge workarounds here. I also refer to his Ph.D. if you are interested in the theory.

Hello,

First off, hopefully this is the right spot for this question -- if not, sorry!

I've been struggling with an issue trying to run an FSI solver in parallel for a while now. It's based off of the standard icoFsiElasticNonLinULSolidFoam, but modified to work in a single rotating reference frame and to use the PISO algorithm for fluid solving. All of the solid mechanics stuff was untouched. For what it's worth, I tried running the standard solver (icoFsi...Foam) and received the same error. The case runs perfectly serially but refuses to run in parallel :mad:

The mesh was created in blockMesh and builds without errors. It is then decomposed using the patchConstrained method, keeping the entire fluid/solid interface patch on one processor:

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM Extend Project: Open Source CFD        |

|  \\    /   O peration     | Version:  1.6-ext                               |

|   \\  /    A nd           | Web:      www.extend-project.de                 |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       dictionary;

    location    "system";

    object      decomposeParDict;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



numberOfSubdomains   4;

method                               patchConstrained;



preservePatches ( AF side1 side2 );



patchConstrainedCoeffs

{

        method                                metis;

        numberOfSubdomains        4;



        patchConstraints

        (

                (side1 0) // cyclicGgi sides

                (side2 0)

                (AF 1) // this is the fluid patch that corresponds to the solid interface

        );

}



metisCoeffs

{

    processorWeights ( 1 1 1 1 );

}



distributed     no;



roots           ( );



// ************************************************************************* //

And a similar file for the solid region, assigning the interface to processor 1 again. However, when trying to run the case I get the following error:

Code:

kellis@desktop:~/foam/kellis-3.2/run/wellsTurbine/FSITests/300FSITest_1022/caseDir/fluid$ mpirun -np 4 SRFPisoFSIFoam -parallel

/*---------------------------------------------------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | foam-extend: Open Source CFD                    |

|  \\    /   O peration     | Version:     3.2                                |

|   \\  /    A nd           | Web:         http://www.foam-extend.org         |

|    \\/     M anipulation  | For copyright notice see file Copyright         |

\*---------------------------------------------------------------------------*/

Build    : 3.2-a4e7bf3b15c0

Exec     : SRFPisoFSIFoam -parallel

Date     : Oct 22 2018

Time     : 22:14:36

Host     : desktop

PID      : 13246

CtrlDict : "/home/kellis/foam/kellis-3.2/run/wellsTurbine/FSITests/300FSITest_1022/caseDir/fluid/system/controlDict"

Case     : /home/kellis/foam/kellis-3.2/run/wellsTurbine/FSITests/300FSITest_1022/caseDir/fluid

nProcs   : 4

Slaves : 

3

(

desktop.13247

desktop.13248

desktop.13249

)



Pstream initialized with:

    nProcsSimpleSum   : 16

    commsType         : blocking

SigFpe   : Enabling floating point exception trapping (FOAM_SIGFPE).



// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Create time



Create dynamic mesh for time = 0



Selecting dynamicFvMesh dynamicMotionSolverFvMesh

Initializing the GGI interpolator between master/shadow patches: side1/side2

Selecting motion solver: laplace

Selecting motion diffusivity: quadratic



Reading transportProperties



Reading field p



Reading field Urel



Creating field Uabs



Creating SRF model



Selecting SRFModel rpm

Reading/calculating face flux field phi



Selecting incompressible transport model Newtonian

Selecting turbulence model type RASModel

Selecting RAS turbulence model kOmegaSST

kOmegaSSTCoeffs

{

    alphaK1         0.85;

    alphaK2         1;

    alphaOmega1     0.5;

    alphaOmega2     0.856;

    gamma1          0.555556;

    gamma2          0.44;

    beta1           0.075;

    beta2           0.0828;

    betaStar        0.09;

    a1              0.31;

    b1              1;

    c1              10;

    F3              false;

}



Reading incremental displacement field DU



Patch AFS        Traction boundary field: DU

        nonLinear set to updated Lagrangian

Patch AFS        Traction boundary field: Usolid

        nonLinear set to updated Lagrangian

Reading incremental displacement field DV



Reading accumulated velocity field V



Reading accumulated stress field sigma



Reading incremental stress field DSigma



Selecting rheology model linearElastic

Creating constitutive model



Reading coupling properties

Create fluid-to-solid and solid-to-fluid interpolators

Check fluid-to-solid and solid-to-fluid interpolators

[1] 

[1] 

[1] --> FOAM FATAL ERROR: 

[1] given field does not correspond to patch. Patch size: 2842 field size: 0

[1] 

[1]     From function PatchToPatchInterpolation::faceInterpolate(const Field<Type> ff)

[1]     in file /home/kellis/foam/foam-extend-3.2/src/foam/lnInclude/PatchToPatchInterpolate.C at line 132.

[1] 

FOAM parallel run aborting

[1] 

--------------------------------------------------------------------------

MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD 

with errorcode 1.



NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.

You may or may not see output from other processes, depending on

exactly when Open MPI kills them.

--------------------------------------------------------------------------

--------------------------------------------------------------------------

mpirun has exited due to process rank 1 with PID 13247 on

node desktop exiting improperly. There are two reasons this could occur:



1. this process did not call "init" before exiting, but others in

the job did. This can cause a job to hang indefinitely while it waits

for all processes to call "init". By rule, if one process calls "init",

then ALL processes must call "init" prior to termination.



2. this process called "init", but exited without calling "finalize".

By rule, all processes that call "init" MUST call "finalize" prior to

exiting or it will be considered an "abnormal termination"



This may have caused other processes in the application to be

terminated by signals sent by mpirun (as reported here).

--------------------------------------------------------------------------

You may have noticed there are cyclic GGI patches involved, but I don't think they are to blame. I swapped them to standard cyclics and the problem still persisted.

The patch size is consistent with the number of faces on the fluid side of the interface patch, and the error originates from the processor I assigned it to. However, I can't find anything wrong with the case files and am not sure what to try next. Has anyone encountered a similar issue before?

Thanks for any assistance!
-Kellis

Quote:

Originally Posted by Kellis (Post 712066)

Code:

/*--------------------------------*- C++ -*----------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | OpenFOAM Extend Project: Open Source CFD        |

|  \\    /   O peration     | Version:  1.6-ext                               |

|   \\  /    A nd           | Web:      www.extend-project.de                 |

|    \\/     M anipulation  |                                                 |

\*---------------------------------------------------------------------------*/

FoamFile

{

    version     2.0;

    format      ascii;

    class       dictionary;

    location    "system";

    object      decomposeParDict;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



numberOfSubdomains   4;

method                               patchConstrained;



preservePatches ( AF side1 side2 );



patchConstrainedCoeffs

{

        method                                metis;

        numberOfSubdomains        4;



        patchConstraints

        (

                (side1 0) // cyclicGgi sides

                (side2 0)

                (AF 1) // this is the fluid patch that corresponds to the solid interface

        );

}



metisCoeffs

{

    processorWeights ( 1 1 1 1 );

}



distributed     no;



roots           ( );



// ************************************************************************* //

And a similar file for the solid region, assigning the interface to processor 1 again. However, when trying to run the case I get the following error:

Code:

kellis@desktop:~/foam/kellis-3.2/run/wellsTurbine/FSITests/300FSITest_1022/caseDir/fluid$ mpirun -np 4 SRFPisoFSIFoam -parallel

/*---------------------------------------------------------------------------*\

| =========                 |                                                 |

| \\      /  F ield         | foam-extend: Open Source CFD                    |

|  \\    /   O peration     | Version:     3.2                                |

|   \\  /    A nd           | Web:         http://www.foam-extend.org         |

|    \\/     M anipulation  | For copyright notice see file Copyright         |

\*---------------------------------------------------------------------------*/

Build    : 3.2-a4e7bf3b15c0

Exec     : SRFPisoFSIFoam -parallel

Date     : Oct 22 2018

Time     : 22:14:36

Host     : desktop

PID      : 13246

CtrlDict : "/home/kellis/foam/kellis-3.2/run/wellsTurbine/FSITests/300FSITest_1022/caseDir/fluid/system/controlDict"

Case     : /home/kellis/foam/kellis-3.2/run/wellsTurbine/FSITests/300FSITest_1022/caseDir/fluid

nProcs   : 4

Slaves : 

3

(

desktop.13247

desktop.13248

desktop.13249

)



Pstream initialized with:

    nProcsSimpleSum   : 16

    commsType         : blocking

SigFpe   : Enabling floating point exception trapping (FOAM_SIGFPE).



// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Create time



Create dynamic mesh for time = 0



Selecting dynamicFvMesh dynamicMotionSolverFvMesh

Initializing the GGI interpolator between master/shadow patches: side1/side2

Selecting motion solver: laplace

Selecting motion diffusivity: quadratic



Reading transportProperties



Reading field p



Reading field Urel



Creating field Uabs



Creating SRF model



Selecting SRFModel rpm

Reading/calculating face flux field phi



Selecting incompressible transport model Newtonian

Selecting turbulence model type RASModel

Selecting RAS turbulence model kOmegaSST

kOmegaSSTCoeffs

{

    alphaK1         0.85;

    alphaK2         1;

    alphaOmega1     0.5;

    alphaOmega2     0.856;

    gamma1          0.555556;

    gamma2          0.44;

    beta1           0.075;

    beta2           0.0828;

    betaStar        0.09;

    a1              0.31;

    b1              1;

    c1              10;

    F3              false;

}



Reading incremental displacement field DU



Patch AFS        Traction boundary field: DU

        nonLinear set to updated Lagrangian

Patch AFS        Traction boundary field: Usolid

        nonLinear set to updated Lagrangian

Reading incremental displacement field DV



Reading accumulated velocity field V



Reading accumulated stress field sigma



Reading incremental stress field DSigma



Selecting rheology model linearElastic

Creating constitutive model



Reading coupling properties

Create fluid-to-solid and solid-to-fluid interpolators

Check fluid-to-solid and solid-to-fluid interpolators

[1] 

[1] 

[1] --> FOAM FATAL ERROR: 

[1] given field does not correspond to patch. Patch size: 2842 field size: 0

[1] 

[1]     From function PatchToPatchInterpolation::faceInterpolate(const Field<Type> ff)

[1]     in file /home/kellis/foam/foam-extend-3.2/src/foam/lnInclude/PatchToPatchInterpolate.C at line 132.

[1] 

FOAM parallel run aborting

[1] 

--------------------------------------------------------------------------

MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD 

with errorcode 1.



NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.

You may or may not see output from other processes, depending on

exactly when Open MPI kills them.

--------------------------------------------------------------------------

--------------------------------------------------------------------------

mpirun has exited due to process rank 1 with PID 13247 on

node desktop exiting improperly. There are two reasons this could occur:



1. this process did not call "init" before exiting, but others in

the job did. This can cause a job to hang indefinitely while it waits

for all processes to call "init". By rule, if one process calls "init",

then ALL processes must call "init" prior to termination.



2. this process called "init", but exited without calling "finalize".

By rule, all processes that call "init" MUST call "finalize" prior to

exiting or it will be considered an "abnormal termination"



This may have caused other processes in the application to be

terminated by signals sent by mpirun (as reported here).

--------------------------------------------------------------------------

Kellis,

I have had a similar problem in the past, with no success finding a solution. I had to resort to performing a full simulation (in my case, 3 blades of a wind turbine). I only post here to add my name to the list of people who are looking for a solution..

Cheers,
Dave

Hi Kellis, Dave,

I suggest running the solver through a debugger (e.g. gdb) in debug mode and finding the lines of code where the solver stops. From there, it should be clear where the error is coming from and then hopefully a solution can be found.

Philip