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