[tomdylan]

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

[bigphil]

Hi Tom,

I am glad the code is now compiling for you.

Is there any documentation out there of the syntax changes during the evolution of o30ext - o40ext?

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/.

Originally Posted by tomdylan

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.

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

[Horst]

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.

[ilhado]

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

[Horst]

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

[tomdylan]

Hi Philip,

Originally Posted by bigphil

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

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

fvm::ddt(p) == fvm::laplacian(coeff_Dp, p) - fvc::div(fvc::ddt(coeff_Dp2,U))
            );
        solverPerfP = pEqn.solve(); 
            p.relax();
            dp = p - p.oldTime();
            dp.correctBoundaryConditions();   

          Kfluid = 1.0/(Sr/K_ + (1.0-Sr)/(p0+psi*gamma+dp));       
          coeff_Dp = (kSat/gamma*Kfluid/n);
          coeff_Dp2 = Kfluid/n;

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:

laplacian(Dp,p) Gauss harmonic corrected;

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

[bigphil]

Hi Tom,

Originally Posted by tomdylan

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

Interesting; is the code shared somewhere?

Originally Posted by tomdylan

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

fvm::ddt(p) == fvm::laplacian(coeff_Dp, p) - fvc::div(fvc::ddt(coeff_Dp2,U))
            );
        solverPerfP = pEqn.solve(); 
            p.relax();
            dp = p - p.oldTime();
            dp.correctBoundaryConditions();   

          Kfluid = 1.0/(Sr/K_ + (1.0-Sr)/(p0+psi*gamma+dp));       
          coeff_Dp = (kSat/gamma*Kfluid/n);
          coeff_Dp2 = Kfluid/n;

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.

You can relax Kfluid just like p:

...
          Kfluid.storePrevIter();
          Kfluid = 1.0/(Sr/K_ + (1.0-Sr)/(p0+psi*gamma+dp));       
          Kfluid.relax();
...

The add the name of Kfluid to the relaxationFactors dictionary within fvSolution.

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:

laplacian(Dp,p) Gauss harmonic corrected;

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.

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

Philip

[ilhado]

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

[shpalman]

Hello,

Originally Posted by will.logie

I've not yet managed to validate thermalStressFoam with the hotCylinder case. Bug fix done.

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

[Thaw Tar]

Originally Posted by Sourabhsup

Hello all,

I am using the icoFsiElasticNonLinULSolidFoam solver. I cannot understand how the traction force (tForce) calculated in setInterfaceForce.H is applied in the solveSolid.H. The solveSolidEuler.H has the governing equation implemented, but i dont see the force term in the governing equation. Can anyone please explain the steps in which solid is solved? How is that traction force applied? and where?

Please guide me on this. Following code is solveSolidEuler.H

{
# include "readSolidMechanicsControls.H"

int iCorr = 0;
lduMatrix::solverPerformance solverPerf;
scalar initialResidual = 0;

lduMatrix::debug = 0;

# include "EulerCoeffs.H"

do
{
DU.storePrevIter();

fvVectorMatrix DUEqn
(
Cn*rho*fvm::ddt(DU)
- Co*rho*DV.oldTime()
==
fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
- fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
+ fvc::div
(
mu*gradDU.T()
+ lambda*(I*tr(gradDU))
+ mu*(gradDU&gradDU.T())
+ 0.5*lambda*(I*tr(gradDU & gradDU.T()))
+ (sigma & DF.T())
+ (DSigma & DF.T()),
"div(sigma)"
)
);

solverPerf = DUEqn.solve();

DU.relax();

if(iCorr == 0)
{
initialResidual = solverPerf.initialResidual();
}

gradDU = fvc::grad(DU);

DF = gradDU.T();

# include "calculateDSigma.H"
}
while
(
solverPerf.initialResidual() > convergenceTolerance
&& ++iCorr < nCorr
);

Info << "Solving for " << DU.name()
<< ", Initial residual = " << initialResidual
<< ", Final residual = " << solverPerf.initialResidual()
<< ", No outer iterations " << iCorr << endl;

DV = fvc::ddt(DU);

lduMatrix::debug = 1;
}

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

[bigphil]

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

[Thaw Tar]

Originally Posted by bigphil

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

[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

[RLFerreira]

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.

[wjw]

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

[Tobi]

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.

[wjw]

Originally Posted by Tobi

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, Phillip 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

[bigphil]

Originally Posted by wjw

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

Originally Posted by wjw

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

[Thaw Tar]

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

[bigphil]

Originally Posted by RLFerreira

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.

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.