[bigphil]

Originally Posted by danny261083

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

[davidsblom]

Originally Posted by bigphil

Ah, I see now, there are three implicit functions defined in backwardbackwardD2dt2Scheme but only the first implemented backward scheme:

        tmp<fvMatrix<Type> > fvmD2dt2 // This is correct
        (
            GeometricField<Type, fvPatchField, volMesh>&
    );

    tmp<fvMatrix<Type> > fvmD2dt2 // this is just Euler
        (
            const dimensionedScalar&,
            GeometricField<Type, fvPatchField, volMesh>&
    );

    tmp<fvMatrix<Type> > fvmD2dt2 // this is just Euler
        (
            const volScalarField&,
            GeometricField<Type, fvPatchField, volMesh>&
        );

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:

                  fvm::d2dt2(rho,DU)

to

                  rho*fvm::d2dt2(DU)

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

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:

fvVectorMatrix DVEqn
(
    fvm::ddt ( rho, DV )
    ==
    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 = DVEqn.solve ( );

fvVectorMatrix DUEqn ( fvm::ddt ( DU ) == DV );

DUEqn.solve();

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:

--> FOAM FATAL ERROR:

    request for volTensorField grad(DV) from objectRegistry region0 failed
    available objects of type volTensorField are

3
(
((1)+grad(DU).T())
grad(DU)
(((1)+grad(DU).T())-(1))
)

Thanks alot.

Best regards,
David Blom

[bigphil]

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

[codder]

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

Two questions:

(1) Does the use of the crackerFvMesh library for the solid, then preclude use of a second dynamic meshing library for an additional fluid domain?

(2) If so, are their steps I can take to segregate the crackerFvMesh from the fluid's dynamicFvMesh?

Not sure if I'm totally making sense; especially, as I initially thought that using icoFsiFoam's tet decomposition would protect me from this issue.

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.

[codder]

@ 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

[bigphil]

Hi Eric,

Originally Posted by codder

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 …")

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.

Originally Posted by codder

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

Using icoFsiFoam as a basis should be fine.

Originally Posted by codder

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.

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.

Originally Posted by codder

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.

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

[codder]

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:

~/foam/foam-extend-3.0/tutorials/solidMechanics/deprecatedTutorials/icoFsiFoam/flappingConsoleSmall/Allrun

... described by this post.

These lines can be deleted from the icoFsiFoam tutorial:

cd constant
ln -s ../../solid/constant solid
cd ..
cd 0
ln -s ../../solid/0 solid
cd ..

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

Best, Eric

[codder]

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

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;

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:

 {
    int cohesivePatchPresSize(cohesivePatchUPtr ? cohesivePatchUPtr->size() : cohesivePatchUFixedModePtr->size());
    for (label i=0; i<cohesivePatchPresSize; i++)
        {   
              vector pressure = 10e9 * mesh.Sf().boundaryField()[cohesivePatchID][i]; 
              //Info << "\nLabel "<< i <<" ."  << nl << endl;
              //Info << "\nPressure vector "<< pressure <<" ."  << nl << endl;

          if(cohesivePatchUPtr)
            {   
              cohesivePatchUPtr->traction()[i] -= pressure;
            }   
          else
            {   
              cohesivePatchUFixedModePtr->traction()[i] -= pressure;
              cohesivePatchUFixedModePtr->initiationTraction()[i] -= pressure;
            }   
    
        }   
}

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

[codder]

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

        const labelList& gcfa = stressMesh.globalCrackFaceAddressing();

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

Internal face to break: 1(1391)
Coupled face to break: 0()
--> FOAM FATAL ERROR: 
problem with defining global crack face addressing

    From function crackerFvMesh::makeGlobalCrackFaceAddressing() const
    in file arbitraryCrack/crackerFvMesh/crackerFvMesh.C at line 400.

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

    const vectorField::subField crackCf =
        boundaryMesh()[crackPatchID_.index()].faceCentres();

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

        if (gcfa[faceI] < 0)
        {   
            FatalErrorIn

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

[davidsblom]

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

[iczajka]

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

[bigphil]

Originally Posted by iczajka

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

[bigphil]

Originally Posted by davidsblom

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

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:

rho = rho / J

This line was probably added during analysis of a quasi-static case as a quick way to check rho changes.

Thanks,
Philip

[Forum]

Hi,

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

Thanks.

[bigphil]

Originally Posted by Forum

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

[Forum]

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.

[bigphil]

Originally Posted by Forum

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

[hxaxtma]

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.

[bigphil]

Originally Posted by hxaxtma

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

[hxaxtma]

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