Thank you very much for this hint, Philip!

Now the solver is compiled without any errors.

And that was the reason for confusion:

lduMatrix::solverPerformance -> lduSolverPerformance

Is there any documentation out there of the syntax changes during the evolution of o30ext - o40ext?
Next time I might check the errors myself.

The modifications I am trying to implement is to consider the dependence of the fluid (gas & water mixture) stiffness K' (equivalent to the fluid compressibility 1/K') from the actual pore pressure. This of course enhances the non-linearity since the pore pressure equation is affected by the volumetric strain increment *and* the change of the fluid stiffness K' in response to pore pressure changes.

First tests indicate numerous iterations and non-convergence after a couple of time steps. However, a "lagged scheme" in which K' from the previous time step is considered (instead of simultaneous evaluation of K' within the iteration loop) yields reasonable results. These of course depend on the time step size selection...

Strongly non linear systems (e. g. unsaturated flow in porous media) sometimes are addressed through Picard-iteration approaches. I wonder if a Picard scheme might improve the convergence of the current problem.

Anyway, I succeeded installing o40x in opensuse leap 42.2. thanks again for your support!

Tom

Hi Tom,

I am glad the code is now compiling for you.

Not that I know of, apart from the code itself; there is some info here: http://foam-extend.fsb.hr/release-of-foam-extend-4-0/.

It makes sense that the lagged scheme is more stable; as regards coupling K changes within the time-step, I suppose you are already using Picard/Fixed-point iterations i.e. successive substitution.
Under-relaxation of the K field will probably help, or more under-relaxation of the p and U fields (equations and fields).

Philip

fsiFoam in extend4.0 / running in parallel
 

Hello fsiFoamers,

is anyone running fsiFoam on ext-4.0 in parallel? In the ext-3.1 version there was this decomposeParFsi-tool, which was easy to use and worked fine, but on ext-4.0 there is nothing like that. The Allrun skript suggest the use of standard decomposePar, but all my parallel jobs break while calling the first solid solver.
It could depend on our machine here so it would be nice, if some of you could give me some short information about your experiences in running the ext-4.0 version of fsiFoam in parallel.

Thanks a lot in advance
Georg

BTW: the robin-bc approach in the ext-4.0 version is much more stable, even for high fluid density! Thanks to everyone who took part in developing, i highly appreciate your work.

Hello Georg,

I've been using fsiFoam for a couple months now and also used version 3.2 before. While for version 3.2 I was only able to decompose my cases with decomposeParFsi, in version 4.0 decomposePar worked great. There are only two issues I found so far:

1) I was not able to visualize the partitions on ParaView and

2) If I decompose from the latest time-step which is different from zero, the decomposition does not create the directory "processor*/constant" (instead, the utility name it "processor*/0")

However, since you say your solid solver is not called properly in run time, did you link the fluid and solid directories? I think this already happened to me, and I remember that I missed the links between fluid and solid. Is there any error message ?

Iago

iago

Hello Iago,

thanks for your reply. Now i know that the problem is somewhere on our machine or local extend installation.

I was succesfully viszualizing the partitions with foamToVTK. It works both in 3.1 and 4.0.

Georg

implementing spatial heterogeneity of the saturated hydraulic conductivity
 

Hi Philip,

I found a quite stable Picard approach in suGWFoam, a saturated unsaturated solver by Liu which could serve as template for a future implementation of coupled deformation saturated/unsaturated flow. I tested Liu's solver under quite mean conditions (flow with large gradients into extremely dry soil and it yielded excellent results (even commercial software would not resolve this problem). To be honest I presently do not feel experienced enough to effectively implement the Picard scheme (which is a time adaptive approach) into elastoPlastoBiotFoam. In view of the high non-linearity of the conductivity and storage properties in unsaturated flow considered by the Picard approach a segregated scheme to consider deformation should not represent a big issue.

I played a little bit around trying to implement the pressure dependence of the Kfluid-field (stiffness of the water gas mixture within the pores) and I came to stable results with the fixed point approach (relaxing only p and U):

This way the pressure changes dp are considered inside the iteration loop. The absolute pressure consists of the atmospheric (p0) and the hydrostatic (psi=pressure head below phreatic surface) components.

Probably relaxing Kfluid (how ?), as you suggest, might enhance convergence.

Now I tried to implement heterogeneity of the saturated hydraulic conductivity, which in natural settings is always an issue, by introducing the field kSat. in open water environments it is not uncommon to find conductivity contrasts spanning 4 orders of magnitude: sand 1^-04 m/s, silt 1^-8 m/s. The coeff_Dp field varies correspondingly over orders of magnitude.

First I decided to adapt the laplacian scheme to harmonic, which makes more sense than a linear one from an geohydraulic point of view:

Unfortunately the solver now diverges during the first iterations.

From the (very illustrative) convergence analysis in "On finite volume method implementation of poro-elastic-plasticity soil model" by Tang, Heledal & Cardiff it is clear, that a very large Kfluid (due to almost compressed gas in 20 m below phreatic surface) becomes stiffer than the bulk modulus of the soil and therefore mesh and time step sizes are determined by the bulk modulus, the water unit weight and the hydraulic conductivity! It turned out that a mesh size of approx. 0.2 m is necessary. I reduced the model to 2D (app. 120 m x 20 m) but could not reach convergence yet.

I could imagine that my problem is not the first one to exhibit strongly varying laplacian terms. So I would appreciate any hints on strategies to overcome these kind of difficulties.

Tom

Hi Tom,

Interesting; is the code shared somewhere?

You can relax Kfluid just like p:
The add the name of Kfluid to the relaxationFactors dictionary within fvSolution.

Hmnn relaxing the p field should help with convergence but you are already doing this; you could try more.

Philip

Hello,

In the new version of the FSI toolkit, a great number of new boundary conditions (BC) was implemented for the fluid part. One of them is called outflowPressure, which I suppose is intended for outlets of the flow domain with a specified pressure on it.

However, I am trying to understand what exactly this BC performs. I looked in its .C and .H files, but so far I only understood that the BC accounts for a correction/modification of the pressure involving a relation with the face normal velocity gradient through a variable called pp:

pp = nu.value()*(n & pU.snGrad())

where

pU = patch().lookupPatchField<volVectorField, vector>("U")

I could not guess more (my C++ is very limited, unfortunately) and, so far, I did not find any information about this specific implementation of this boundary condition in the forums or even comparing with the implemented BC's in other versions of foam-extend.
Does anyone knows what exactly this BC implements?

Thanks in advance for any help :)
Iago

elasticThermalSolidFoam
 

Hello,

I'm trying to use foam-extend-4.0, but it seemed like tractionBoundaryGradient.C still had "if (!incremental)" here.

What I'm actually trying to do is model steady-state elastic strain in semiconductor structures.[1] There's a "built-in" strain which comes from a "misfit" between two crystals with different atomic spacings. I had been using solidDisplacementFoam for this, and it basically worked, although I was faking the misfit strain with a thermal strain (setting the thermal conductivity to zero).

Basically, elasticThermalSolidFoam does not give the expected result (the solidDisplacementFoam result is correct) which is why I'm curious about the status of that !incremental bug, or the outcome of this thread.

Ideally I'd like to be able to properly specific a misfit strain somehow, i.e. that two different (orthotropic) materials are obliged to maintain the same lattice spacing at the interface despite having different natural lattice spacings, but for now the fake thermal strain thing would be ok.

Thanks!

1. http://aip.scitation.org/doi/abs/10.1063/1.4765009

Hi, could anybody please explain me about this question, too. I also have this very confusion. I cannot figure out how the computed fluid force (traction force) is used in the solid solver.

Thaw Tar

Hi Thaw Tar,

The fluid force/traction is applied to the solid as a traction boundary condition in the "setInterfaceForce.H" header file (just before "solveSolid.H" in the "icoFsiElasticNonLinULSolidFoam.C" file).

Philip

Thank you very much Dr. Cardiff for your very fast reply and for your contribution to the solid solver in OpenFOAM.

Now I understood how the fluid force was used as a boundary condition for the solid solver.

And also for the advice to increase the solid solver iteration count and to control the convergence with residuals. It worked very well!!!

Thaw Tar

Dear Dr. Cardiff and others,

I am sorry for disturbing again. But I have one problem with the fsiFoam solver.

I ran the simulation for some time steps and wrote out the data. But when I restart the simulation from the last time step, it showed this following error:

Create extended GGI zone-to-zone interpolator
Checking fluid-to-solid face interpolator
Fluid-to-solid face interpolation error: 0.0126103
Checking solid-to-fluid point interpolator (GGI)
Extended GGI, master point distance, max: 1e+15, avg: 2.30143e+14, min: -0.002
Extended GGI, master point orientation (<0), max: 1, min: -1, nIncorrectPoints: 638/982


--> FOAM FATAL ERROR:
Master point addressing is not correct

From function ExtendedGGIInterpolation::masterToSlavePointInterp olate(const Field<Type> pf)
in file ./numerics/ggi/ExtendedGGIInterpolation/ExtendedGGIInterpolation.C at line 703.

FOAM aborting
Aborted (core dumped)

If anyone has the solution and could explain to me how to fix it, I would be very grateful.

Thaw Tar

Dear Formers, I trying simulate a 3D valve compressors. Initially, we created a 2D simple model with a cantilever beam, which will be a valve, as showed by mesh picture attached. At a certain time of simulation, the FSI interface doesn't match with solid mesh anymore (Uncoupled_meshes.png). In addition, the simulation does not diverge and the solution proceeds with those uncoupling meshes. Can someone help me? I know that if I try to restart the simulation from the last time step the error "Master point addressing is not correct" occurs.

Thanks.

How to get the contact area and normal contact stress?
 

Hello Everybody,
Nice to be here!
I am a phD student in Germany facing with a task of simulating and calculating dry sliding wear in a total hip replacement in a long term. As is all known, FEM method in Abaqus is a common and relative simple approach to perform this task. So far, I have not found any paper which used FVM method to simulate dry sliding wear. But with solidmechanics solver in combination with dynamical mesh technology in foam extend can the same goal be theoretically achieved. My plan is finally to publish a paper with both fem and fvm approach, and show that fvm can be a good alternative approach.
Now, I am at the beginning stage and already tried NonLinUL solvers to build some simple contact models. The results are pretty good compared with abaqus. The problem I am confonted with is that i am not albe to get the contact area and normal contact stress from the soildmechanics solver in foam-extend , but I really need these two for simulating wear. Can I get some help or cooperation from here? Since these two results are very common and important for contact problems and can usually be directly obtained from some fem softwares, I am looking forward to seeing these two can also be included in foam-extend 5.0.

Best regards,
Jingwei

Hi,

I guess Phil can help you more based on the fact that he developed it but since I am using an own solid stress solver, it should be not a big deal to get the contact area. In the contact class there should be some information about the cells / faces which are in contact and thus, you should easily be able to calculate the area. The normal contact stress should be similar. However, I do not know how the normal contact stress is given in terms of mathematical expressions. Again, Philip should be the one who could give more fundamental insight.

I just wanted to let you know that it is possible.

Hi Tobias,
Thanks for your answer. Actually my German name is also Tobias. I agree with you that the key point is to get the information about which cell/surfaces are really in contact. The rest things are pretty easy. By the way, If you are interested in the mathematical expression of normal conact stress from FEM method which is vertex-based, you can send me a PM, I can send you a link or a paper.

Jingwei

Hi Jingwei,

Further to Tobias's response:

The topic of my PhD was on modelling a hip joint in OpenFOAM, so some of the following articles may be useful:

Hi Jingwei,

Further to Tobias's response:

The topic of my PhD was on modelling a hip joint in OpenFOAM, so some of the following articles may be useful:

• Development of a Hip Joint Model for Finite Volume Simulations
• Development of mapped stress-field boundary conditions based on a Hill-type muscle model
• TOWARDS IN SILICO ANALYSIS OF TOTAL HIP ARTHOPLASTY MECHANICS (recent work)

Also, you can find my PhD on my ResearchGate, which may be useful.

As regards calculating the normal contact stress (), it is not written by default but you can write a utility to do this (send me a PM with your email and I will send you a utility that will work).

For the contact area, it is possible to do it in a similar way.

Philip

VIV of a finite length vertical cylinder by fsiFoam
 

Dear all,

I tried to simulate the fluid structure interaction of a finite length vertical cylinder with a fixed bottom using fsiFoam and succeeded. It could simulate the occurrence of VIV and lock-in.

https://www.youtube.com/watch?v=NnFrh9ygWos

Special thanks to Dr. Cardiff for developing the solid solver and his advice to tighten the solid solver by increasing the iteration count.

I will test on more cases of reduced velocities.

Thaw Tar

Hi RLFerreira,

If the solid and fluid interface do not match then this typically means that the FSI loop did not converge to a tight enough tolerance (or did not converge at all) in each time-step.

I would suggest looking at the FSI residual and check if it is converging in each time-step; if it is not then try a larger maximum number of FSI correctors or a different FSI coupling method or smaller time-steps; if it is converging, then try a tighter FSI tolerance.

Separately, the segregated solid solver is typically not very efficient for slender geometries (like your thin beam): it much prefers "blocky" geometry. I understand you probably can't change your geometry, but if you PM me your email address then I can send you a work-in-progress coupled solid solver that may be more efficient/robust in this case.

Best,
Philip

Also, nice simulation Thaw Tar.

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)?

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:

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?

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!!!

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

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

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

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.

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:

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:

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

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