thanks Dr. Cardiff

As for the last comment, yes you are right. Now I am clamping those patches with fixed displacement in every direction.

I will work on FSI patch boundary condition. Actually the solver does not allow any other boundary conditions exept fixedDisplacement types and solidTreaction.

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

edit: I have talked to the author of that paper, they have used prescribed pressure and traction for FSI interface and No load for outer wall.

How can I apply No load boundary condition in OpenFOAM?


regrads

Hi Philip,

Thanks for your answer. At this moment, i cannot able to plot the tip deflection from the solver. May I have the code for generate the file as "historySolid.dat"? Very sorry for keep troubling you and thanks a lot for your help.

:)

Regards,
Vincent

Hi Philip,

sorry for troubling you with questions again. :(
I tried to plot with smaller time step, but it still gave the same error message as:

ERROR: In /home/vin-82/OpenFOAM/OpenFOAM-1.6-ext/ThirdParty/rpmBuild/BUILD/ParaView-3.12.0/VTK/Filtering/vtkTable.cxx, line 353
vtkTable (0x9676cb0): Column "Probe Coordinates (Magnitude)" must have 1 rows, but has 2.

Please help me on plotting the tip deflection :)
Thanks.

Regards,
Vincent

Hello,

while trying to run the icoFsiElasticNonLinULSolidFoam solver in parallel as described in the info file of the tutorial I experience the exact same error as the error Daniel mentioned earlier.
After that I recompiled the solver with the moveSolidMeshLeastSquares.H file from Philip but the simulation keeps stopping at 0.065s.
As I'm new to OpenFoam I don't know how to fix this problem and would really appreciate some help.


Best regards

Michael

Hi Philip,

I'm currently validating the FSI solver, and as a first step I am trying to run the CSM3 test from the cylinder flap benchmark. So only a structure solver is used. I was wondering which settings to use, because with my current setup, the displacement of the flap is much larger than expected. I have used the elasticGravitySolidFoam solver.

What is the best approach for this problem? And I was wondering which divSigmaExpMethod to use. I'm currently using the standard approach, but I'm not aware of the differences between the different settings settings/surface/decompose/laplacian.

Best regards,
David

Hi David,

All divSigmaExp methods should produce very similar results; in theory "decompose" is the most accurate and I have found that "surface" is a bit more robust on bad meshes.

As regards gravity, the elasticGravitySolidFoam is a small strain solver so it's not appropriate for such large rotations as in the benchmark.
Essentially you should take elasticNonLinULSolidFoam and add the body force increment in the first time step (i.e. small modification to the solver).

Best regards,
Philip

EDIT:

Add something like this after DUEqn is created and before DUEqn.solve() is called:

Hi Philip,

Thanks for your quick response. This fix seems to have solved the problem. At first I also tried to add the gravity force at each time step before trying the elasticGravitySolidFoam solver. What is the reason that the gravity forces is only applied during the first time step? It seems counterintuitive.
Thanks alot.

Another question. Currently, a slightly different solid solver is used in the icoFsiElasticNonLinULSolidFoam fluid-structure interaction solver. Are there a lot of differences between the two solid solvers? I'm considering of replacing the structure solver of icoFsiElasticNonLinULSolidFoam with the elasticNonLinULSolidFoam solver.
Thanks alot!

Best regards,
David

The momentum equation in the elasticNonLinULSolidFoam/FSI solvers are cast in their incremental form i.e. the increment of momentum is conserved for each time-step, the primary variable is then the increment of displacement DU.
Therefore the gravity body force need only be added in the first increment and then you are only solving for changes in momentum/forces (gravity force stays constant).
Hope this makes sense.

The FSI solid solver uses the same approach so they are essentially the same, maybe just a few minor differences.

Philip

Hi Philip,

I've calculated the first 4 seconds of the CSM3 cylinder flap benchmark, and it seems that the displacement of the flap is growing over time. I've attached a plot of the displacement over time. Any idea what is causing the "growing/diverging"?
I have similar problems with the fluid structure interaction case where the displacement grows over time. Possibly because the mesh is deformed at every time step?
I've tried to use the elasticNonLinTLSolidFoam solver, but this solver crashes directly at the first time step.
Thanks alot for your help.

Attachment 27293

It looks like the amplitude of the displacement oscillation is getting smaller (not growing), is that what you meant?

It would make sense if the oscillation was getting smaller: this would be due to the use of Euler time discretisation which is first order accurate in time.
Two possible remedies:

• smaller time steps
• more accurate time discretisation, such as backward d2dt2 scheme.

Philip

Thanks. I was thinking that it was caused by the time integration scheme. Unfortunately, backward differencing d2dt scheme is not implemented in OpenFOAM.

I'm currently looking into the structure solver used by the fluid-structure interaction solver, and I've created a seperate solver with just the structure code. Can I add the same code snippet to this solver in order to simulate the same CSM3 problem?
Thanks in advance.

Hi David,

The backward d2dt2 scheme is implemented in the forthcoming next release of OpenFOAM Extend (the nextRelease git branch can be downloaded).
The code is found at:
OpenFOAM-1.6-ext/src/finiteVolume/finiteVolume/d2dt2Schemes/backwardD2dt2Scheme

Alternatively a smaller time-step should help.

Yes that should work, assuming you create all the relevant fields.

Philip

I have solved my previous problem but have a new one. I used the github for the official OpenFOAM 2.2x (I am using OpenFOAM-2.2.2). I tried to run the plate hole tutorial by using the blockMesh command then the elastricSolidFoam command, I then get the following error:

It cant seem to find the extendedLeastSquares grad scheme, is this something that needs to be changed or have I installed this feature incorrectly?

Thank you.

Hello Dr. Cardiff

I have a question about running a turbulence case with FSI solver. I have successfully changed the fluid solver according to pisoFoam solver. Now the FSI solver is able to model a turbulence flow.

However I have a little problem in running the tutorial case in turbulence regime. At the first time step the pressure equation leads to very wrong values for pressure distribution in the fluid domain. Also I have tested the fluid part of your case separately with pisoFoam solver and I can see the same problem there. So I think that the turbulence FSI solver is working fine (at-least working like pisoFoam solver).

This problem leads to the divergence of solid solver in FSI solver.

I have done a few tests with pisoFoam solver and it seems that the pressure distribution becomes reasonable after a few time steps ( 4 or more ). Also the steady state solution of the case is in well agreement with fluent.

So I am guessing that the solver just needs a few time step to get rid of that wrong prediction. This is not going to produce any problem in fluid domain simulation, but FSI simulation will definitely diverges because of this. So I have changed the strongly coupled algorithm as follow:

This will skip setting the pressure on solid side for a few time steps (4). By adding the above modification I am now able to run the tutorial case in turbulence regime.

Is there any side effects (skipping a few time steps) that I need to know?


regards

@ Daniel_Khazaei

Could you expand on this a bit, if possible? For instance, does that mean you were able to apply a fixedDisplacement boundary condition using a weak coupling on the FSI solver.

I'm looking into the feasibility (???) of using a dynamic mesh with this new class of solid domain solvers. Please excuse my ignorance I have attempted to put this question in the correct terms.

Hi Philip,

I'm currently testing foam extend 3.0. When running the CSM test of the cylinder flap benchmark with the elasticNonLinULSolidFoam solver, the correct solution is found. However, both the first order and second order time integration scheme give the same results. The order of the backward differencing scheme seems to be one.
Any idea what is causing this?

Also, I've tested the structure solver which is used by the fluid-structure interaction solver and looked at the two time integration schemes. The first order scheme gives the same results as the elasticNonLinULSolidFoam solver. But the second order backward differencing scheme does not give the correct result. Zero-th order behaviour is shown for the total displacement U, whereas the first order scheme shows the correct order.
Any ideas?

Many thanks.

Best,
David

Hi David,

Hmnn, I haven't played with the backwardD2Dt2 scheme much but it should work: it was developed and implemented by Željko Tukovic as outlined in
this paper. See Fig 5 and Fig 6 comparing Euler and backward schemes; the backward scheme certainly reduces the dissipation, although I don't thing that second order accuracy was actually shown. Did you examine the 2nd Order slope of the error plot?
Did you notice that the backward scheme is more accurate (irregardless of order)?

As regards the structural solver within icoFsiElasticNonLinULSolidFoam, it should be the same but it is fine to comment out all the time scheme stuff and just use the built-in backward/Euler schemes like elasticNonLinULSolidFoam.

Best regards,
Philip

Hi,

Thanks for your quick response. Since the first order and second order scheme are giving me the exact same response, I looked at the source.

This is the code in backwardD2dt2Scheme.C:

template<class Type>
tmp<fvMatrix<Type> >
backwardD2dt2Scheme<Type>::fvmD2dt2
(
const dimensionedScalar& rho,
GeometricField<Type, fvPatchField, volMesh>& vf
)
{
tmp<fvMatrix<Type> > tfvm
(
new fvMatrix<Type>
(
vf,
rho.dimensions()*vf.dimensions()*dimVol
/dimTime/dimTime
)
);

fvMatrix<Type>& fvm = tfvm();

scalar deltaT = mesh().time().deltaT().value();
scalar deltaT0 = mesh().time().deltaT0().value();

scalar coefft = (deltaT + deltaT0)/(2*deltaT);
scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);

scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);

if (mesh().moving())
{
notImplemented
(
type()
+ "::fvcD2dt2"
+ "("
+ "const dimensionedScalar& rho, "
+ "GeometricField<Type, fvPatchField, volMesh>& vf"
+ ")"
);
}
else
{
fvm.diag() = (coefft*rDeltaT2)*mesh().V()*rho.value();

fvm.source() = rDeltaT2*mesh().V()*rho.value()*
(
(coefft + coefft00)*vf.oldTime().internalField()
- coefft00*vf.oldTime().oldTime().internalField()
);
}

return tfvm;
}

This is the code in EulerD2dt2Scheme.C:

template<class Type>
tmp<fvMatrix<Type> >
EulerD2dt2Scheme<Type>::fvmD2dt2
(
const dimensionedScalar& rho,
GeometricField<Type, fvPatchField, volMesh>& vf
)
{
tmp<fvMatrix<Type> > tfvm
(
new fvMatrix<Type>
(
vf,
rho.dimensions()*vf.dimensions()*dimVol
/dimTime/dimTime
)
);

fvMatrix<Type>& fvm = tfvm();

scalar deltaT = mesh().time().deltaT().value();
scalar deltaT0 = mesh().time().deltaT0().value();

scalar coefft = (deltaT + deltaT0)/(2*deltaT);
scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);

scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);

if (mesh().moving())
{
scalar halfRdeltaT2 = 0.5*rDeltaT2;

scalarField VV0 = mesh().V() + mesh().V0();

scalarField V0V00 = mesh().V0() + mesh().V00();

fvm.diag() = rho.value()*(coefft*halfRdeltaT2)*VV0;

fvm.source() = halfRdeltaT2*rho.value()*
(
(coefft*VV0 + coefft00*V0V00)
*vf.oldTime().internalField()

- (coefft00*V0V00)*vf.oldTime().oldTime().internalFi eld()
);
}
else
{
fvm.diag() = (coefft*rDeltaT2)*mesh().V()*rho.value();

fvm.source() = rDeltaT2*mesh().V()*rho.value()*
(
(coefft + coefft00)*vf.oldTime().internalField()
- coefft00*vf.oldTime().oldTime().internalField()
);
}

return tfvm;
}

Seems to be the same... Not sure how this is possible for a higher order scheme.
Any clues? To my knowledge, mesh().moving() returns false, since the dynamic mesh class is not used.

Ah, I see now, there are three implicit functions defined in backwardbackwardD2dt2Scheme but only the first implemented backward scheme:
The code should give a not implemented error if either of the bottom two functions are called but the notImplemented functions were not added: I will let them know.

So, for now, you can try out the backward d2tdt2 scheme by changing:
to
in the solver code: this essentially assumes constant density for the calculation of inertia which is probably fine in most case.
It is probably straight-forward to implement the second two functions if you get an urge ;)

Best regards,
Philip

Regarding elasticNonLinTLSolidFoam (cantilever problem)
 

Hi Dr. Cardiff,

I was attempting to run the tutorial problem involving large deflections of a cantilever beam using the total lagrangian solver 'elasticNonLinTLSolidFoam' (in OpenFoam version 2.2.x). I received an error message stating that the extended least squares method was not a recognized solution method. I attempted to run the simulation using least squares, which resulted in small displacements, as can be seen in the enclosed image. I would like to know whether I might be implementing the test case incorrectly.

Thanks

Hi Danny,

OK, the "leastSquares" gradScheme is equivalent to "extendedLeastSquares 0" on orthogonal meshes, so that's fine. Actually, as a side point, I have yet to try out the new point/face least squares in official OpenFOAM, it would be interesting to see how good they are on unstructured meshes.

Also, you say small displacements but your mag(U) is almost 1 metre so that seems quite big to me; if you are referring to the mesh motion being small, then this is because the Total Lagrangian solver uses a non-moving mesh (it integrates over the reference configuration ie mesh); to see the actual deformation use the "Warp By Vector" filter in ParaView and uses the total displacement field "U" and a scale factor of 1.

Another point: your mesh has only one cell across its thickness (top to bottom): you will need at least 2 cells. This is because the computational nodes lie at the cell centres so you need at least 2 to approximate shear strains and the variation of tensile/compressive stress.

Best regards,
Philip

Hi Philip,

Thanks for all your help. I was wondering if it is possible to split up the time integration of DU into two steps. So first to solve for the velocity DV, and thereafter for the displacement DU.

I have tried the following:

But a fatal error is thrown when evaluating DVEqn. Am I doing something wrong or is this not possible with the current code setup.
The following error is thrown:

Thanks alot.

Best regards,
David Blom

Hi David,

I don't think it is possible to do it that way because there must be just one primary variable in the equation, either DV or DU but not both because the linear matrix will be A*DU = B OR A*DV = B; you can't mix and match.
It is probably possible to use implicit Euler and then iteratively correct the equation with explicit source term which is the explicit difference between Euler and backward (and it could be calculated using a mix of velocities or whatever).

However, my previous suggestion of rho*d2dt2(DU) should work fine and quite accurately, as long as you don't have large changes in rho: this is actually a good assumption here because the implemented constitutive law is St Venant Kirchhoff model (nonlinear version of Hooke's law) which doesn't actually match most materials for large elastic/elastoplastic strains.

If you really do want it, the derivation/implementation of d2dt2(rho, DU) is actually quite straight-forward.

Best regards,
Philip

using crackerFvMesh and another dynamic mesh
 

@ Phillip Cardiff -

First off thanks for your efforts packaging the elasticAcpSolidFoam and other solid mechanics solvers into foam-extend.

I experienced problems adapting that solver, for FSI, which I'm wondering if result from a basic conceptional mistake on my part, resulting from having read (oddly enough) the block coupled report.

Previously I wrote:
Edit #1
I have checked with a friend. There does not always seem to be a problem with two conflicting dynamic meshes.

Edit #2
Still puzzled by meshing differences between solid and fluid. However, suffice to say:

(1) If you are creating a FSI solver along the lines of icoFsiFoam
(2) Make sure to change the word in "createStressMesh.H" from "default" to "solid".

... problem solved.

broad elaticAcpQuestions
 

@ Dr. Cardiff

Three broad questions about the craking solver contributions from UCD, in my mind.

1. If you yourself were to attempt an FSI solver using the elaticAcp suite of solvers, how would you go about this?
(In broad-brush general terms, e.g.: "I would definitely use UCD's own FSI solvers, because ...")

2. Is there any fundamental, conceptual problem with using the icoFsiFoam methods for the above purpose?
(Looked at the solid models in foam-extend-3.0/src, but I couldn't resolve this in my mind...)

3. On defining patches:

a. What is the function of the 0/materials boundaryField, in crackingBiMatDcbLinear tut?

... I think they don't do anything, because just defines material heterogeneity on a per-cell basis.

b. How to visualize, for example with paraFoam, the 0/U crack patch?

... I am confused about the nature of the computed values for crack - e.g., as to what "1" or "2" means.


Thanks much for helping me understand, Eric

Hi Eric,

Finite volume crack-FSI solvers have been used successfully in previous analyses (e.g. see here) and they have used similar methods to elasticAcp, so I suppose they should be suitable.
As I see, the main difficulties will be how the fluid is represented inside the cracks; there are (at least) two options:

• Do not solve fluid equations inside cracks, just represent the fluid pressure as a boundary condition on the crack faces;
• Add fluid cells between the cracked faces to get more exact fluid pressure.

The second option may require quite a bit of coding.

Using icoFsiFoam as a basis should be fine.

The material properties are defined in constant/rheologyProperties; when multiMaterial is used (i.e. when there are more than one material present) then you must define/set the 0/materials field, where 0 corresponds to the first material in constant/rheologyProperties, 1 corresponds to the second material in constant/rheologyProperties, etc.

The crackingBiMatDcbLinear case has two materials defined: you can view the materials field in ParaView.

To visualise the crack, turn off the internalField and just turn on the crack patch in ParaView then you will see faces added to it in time.
As regards "1" and "2":

• a faces with "1" means that the face is damaging; this means the face is within the cohesive zone and is currently is the process of dissipating energy. Comparing with standard FE methods, this corresponds to when the cohesive element has initiated;
• a face with "2" means that the face is fully cracked i.e. the traction if zero, unless it is in contact.

Best regards,
Philip

links in Allrun for solidMechanics/depreciated/icoFsiFoam
 

Hi Dr. Cardiff -

I'm happy to put this on "Mantis", but I thought this might be the right place, cause it's not a core code issue. Just noticed that there's a minor duplication of work in:

... described by this post.

These lines can be deleted from the icoFsiFoam tutorial:

I myself have run the tutorial successfully after having deleted them.

Best, Eric

updates to cohesivePatch BC
 

@ Dr. Cardiff -

Thank you answering my questions, as well as the reference. I have contained in this and the following posst two questions about the code contained in elasticAcpSolidFoam.


QUESTION #1


How is it that I can represent the propagation of fluid, as a pressure boundary condition, inside a elasticAcpSolidFoam crack?


SOLUTION TO QUESTION #1


Basically, my trouble was not just to update the BC, but also to do so for each face of an interactively expanding patch.** Initially, I attempted to accomplish this by editing the patch source. But there's another way! I have worked out how to do this, from also another post by Dr. Cardiiff.

The first step involves getting a pointer to the expanding patch (documented in elasticAcpSolidFoam). Then, you go in and update the BC, by inserting a #include on the top of the solver loop... e.g., for me, the loop in which UEqn is defined.

My #include looks like:

I attach documentation of my misguided attempts at FvPatch editing as a warning to others :)


Best, Eric


**Also, with BC values specified by a non-uniform function, S.T. traction= traction[i].

guestion on globalCrackFaceAddressing
 

@ Dr. Cardiff -

The following question relates to an unexplained error I see during runtime of myIcoFsiElasticAcpSolidFoam, specifically relating to use of the crackerFvMesh library. Unlike Question #1 (above), I'm not exactly sure where to start on this - I believe because it seems to require a deeper understanding of meshing.

QUESTION #2

I have acted based upon your (tacit) approval of an icoFsiFoam modification. The (attached) "myIcoFsiElasticAcpSolidFoam" solver now complies and will execute over both fluid and solid case -- IFF the original crackingBiMatDcbLinear displacement boundary conditions have been set to zero (see attached images).

I was worried about runtime errors on stressedMesh.update() at the first crack event -- but no, that works. Instead, I'm seeing a runtime error at line in updateCrack.H:

The error corresponds to the first cracked face. In my log:

I have traced the error in file arbitraryCrack/crackerFvMesh/crackerFvMesh.C to (I think) a negative value of

... which gets checked for negativity at line 395:

I'm wondering what I'm doing wrong?

I though my code was fully segregated between (fluid) mesh and (cracking) stressedMesh - thus boundaryMesh() would produce identical results for "myIcoFsiElasticAcpSolidFoam" as does elasticAcpSolidFoam.

Thanks for any help, Eric

Hi Philip,

Small questions, I noticed that in writeFields.H of the elasticNonLinTLSolidFoam solver the density is updated with rho = rho / J.
I was wondering whether this is correct, since it is not executed at every time step, since runTime->outputTime() does not return true for every time step.
Any thoughts on this?

Thanks.

David

displacement limit
 

Hi, everybody,
I am trying to solve fsi problem with OpenFoam extended 3.0.
I'm using isoFsiElasticNonLinULSolidFoam solver. I have kind of beam, that is moving due to fluid pressure. But also there is support, that should limit displacement of beam.
Is there easy way to implement displacement limit in beam ?

Irek

Hmnn, there is always a way but I'm not sure if there is an easy way.
It sounds like you essentially want to activate a contact-style condition.

You would have to do a bit of coding in the solver and check the solid FSI interface if it goes outside the specified bounds; if it does then you would need to fix the normal displacement or apply a penalty force or something along those lines.

Philip

Hi David,

You are correct: this is a mistake.

rho should not be updated at all in the total Lagrangian solver, as it integrates over the initial mesh and hence always uses the initial density.

So this line should be removed:
This line was probably added during analysis of a quasi-static case as a quick way to check rho changes.

Thanks,
Philip

Roller BC in Simply Supported Beam
 

Hi,

I would like to know how a roller BC of a simply supported beam, the right BC in attachment, can be modeled.

Thanks.

If you apply a pressure to the top surface of the beam then it will bend downwards and the roller will not play any part, if I understand your drawing correctly.

If the pressure was on the bottom surface of the beam then there would be contact between the beam and the roller; in that case you would need to use a contact boundary condition: take a look at the elasticSolidFoam tutorials which use the solidContact boundary condition.

Philip

Simply Supported Beam
 

How can the right BC of attached beam be modeled, please? The right BC is like a pin so beam can rotate around the connection point.

Thanks.

I am not sure if I understand, but if you make a very small boundary patch where the pin is attached, with only one face, then you can set the displacement of this patch/face to be fixedValue ( 0 0 0 ): this will act like pin which the beam can rotate about.

Philip

Hi small question,

I am using the flappingConsole icoFsiFoaam Tutorial in Couette Channel.Therfore I used cyclicGgi for inlet and outlet. How do I have to define these boundaries in the motionU file, cause cyclicGgi is not supported.

I don't have much experience with cyclicGgi so I'm afraid I can't help. Maybe someone else has ideas.

Philip

The problem is that I get backflow at the outlet, I tried in motionU file for inlet and outlet
fixedValue and slip; both are resulting in backlow at the outlet