CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Support thread for "Solid Mechanics Solvers added to OpenFOAM Extend"

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree18Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   January 29, 2014, 08:17
Default
  #141
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 599
Rep Power: 20
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by danny261083 View Post
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
bigphil is offline   Reply With Quote

Old   January 29, 2014, 08:47
Default
  #142
New Member
 
Join Date: Mar 2013
Posts: 17
Rep Power: 5
davidsblom is on a distinguished road
Quote:
Originally Posted by bigphil View Post
Ah, I see now, there are three implicit functions defined in backwardbackwardD2dt2Scheme but only the first implemented backward scheme:
Code:
        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:
Code:
                  fvm::d2dt2(rho,DU)
to
Code:
                  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:

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

Code:
--> 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
davidsblom is offline   Reply With Quote

Old   January 29, 2014, 10:42
Default
  #143
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 599
Rep Power: 20
bigphil will become famous soon enoughbigphil will become famous soon enough
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
bigphil is offline   Reply With Quote

Old   January 29, 2014, 13:56
Default using crackerFvMesh and another dynamic mesh
  #144
Member
 
Eric Bryant
Join Date: Sep 2013
Location: Texas
Posts: 44
Rep Power: 4
codder is on a distinguished road
@ 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:
Quote:
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.
bigphil likes this.

Last edited by codder; February 4, 2014 at 21:06.
codder is offline   Reply With Quote

Old   February 4, 2014, 21:46
Default broad elaticAcpQuestions
  #145
Member
 
Eric Bryant
Join Date: Sep 2013
Location: Texas
Posts: 44
Rep Power: 4
codder is on a distinguished road
@ 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
codder is offline   Reply With Quote

Old   February 5, 2014, 06:10
Default
  #146
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 599
Rep Power: 20
bigphil will become famous soon enoughbigphil will become famous soon enough
Hi Eric,

Quote:
Originally Posted by codder View Post
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.

Quote:
Originally Posted by codder View Post
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.

Quote:
Originally Posted by codder View Post
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.

Quote:
Originally Posted by codder View Post
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
bigphil is offline   Reply With Quote

Old   February 8, 2014, 01:26
Default links in Allrun for solidMechanics/depreciated/icoFsiFoam
  #147
Member
 
Eric Bryant
Join Date: Sep 2013
Location: Texas
Posts: 44
Rep Power: 4
codder is on a distinguished road
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:

Code:
~/foam/foam-extend-3.0/tutorials/solidMechanics/deprecatedTutorials/icoFsiFoam/flappingConsoleSmall/Allrun
... described by this post.

These lines can be deleted from the icoFsiFoam tutorial:

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

Last edited by codder; February 9, 2014 at 19:47.
codder is offline   Reply With Quote

Old   February 13, 2014, 01:05
Default updates to cohesivePatch BC
  #148
Member
 
Eric Bryant
Join Date: Sep 2013
Location: Texas
Posts: 44
Rep Power: 4
codder is on a distinguished road
@ 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?

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

Code:
 {
    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].
Attached Files
File Type: txt warning.txt (4.0 KB, 1 views)

Last edited by codder; March 8, 2014 at 15:01. Reason: solved question, hopefully
codder is offline   Reply With Quote

Old   February 23, 2014, 15:37
Default guestion on globalCrackFaceAddressing
  #149
Member
 
Eric Bryant
Join Date: Sep 2013
Location: Texas
Posts: 44
Rep Power: 4
codder is on a distinguished road
@ 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:

Code:
        const labelList& gcfa = stressMesh.globalCrackFaceAddressing();
The error corresponds to the first cracked face. In my log:

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

Code:
    const vectorField::subField crackCf =
        boundaryMesh()[crackPatchID_.index()].faceCentres();
... which gets checked for negativity at line 395:

Code:
        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
Attached Images
File Type: jpg t1.jpg (57.2 KB, 48 views)
File Type: jpg t2.jpg (68.0 KB, 57 views)
Attached Files
File Type: gz myIcoFsiElasticAcpSolidFoam.tar.gz (19.7 KB, 13 views)
File Type: gz tutorial.tar.gz (7.4 KB, 13 views)
pizicai likes this.
codder is offline   Reply With Quote

Old   February 27, 2014, 10:55
Default
  #150
New Member
 
Join Date: Mar 2013
Posts: 17
Rep Power: 5
davidsblom is on a distinguished road
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
davidsblom is offline   Reply With Quote

Old   March 5, 2014, 09:44
Default displacement limit
  #151
New Member
 
Ireneusz Czajka
Join Date: Nov 2013
Posts: 6
Rep Power: 4
iczajka is on a distinguished road
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
iczajka is offline   Reply With Quote

Old   March 6, 2014, 06:08
Default
  #152
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 599
Rep Power: 20
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by iczajka View Post
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 is offline   Reply With Quote

Old   March 6, 2014, 06:21
Default
  #153
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 599
Rep Power: 20
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by davidsblom View Post
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:
Code:
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
bigphil is offline   Reply With Quote

Old   March 6, 2014, 06:52
Default Roller BC in Simply Supported Beam
  #154
New Member
 
Join Date: Mar 2014
Posts: 8
Rep Power: 4
Forum is on a distinguished road
Hi,

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

Thanks.
Attached Files
File Type: docx Beam Shape.docx (30.1 KB, 13 views)
Forum is offline   Reply With Quote

Old   March 6, 2014, 08:01
Default
  #155
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 599
Rep Power: 20
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by Forum View Post
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
bigphil is offline   Reply With Quote

Old   March 6, 2014, 09:30
Default Simply Supported Beam
  #156
New Member
 
Join Date: Mar 2014
Posts: 8
Rep Power: 4
Forum is on a distinguished road
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.
Forum is offline   Reply With Quote

Old   March 6, 2014, 09:41
Default
  #157
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 599
Rep Power: 20
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by Forum View Post
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
bigphil is offline   Reply With Quote

Old   March 6, 2014, 10:12
Default
  #158
Member
 
Join Date: Jan 2014
Posts: 93
Rep Power: 4
hxaxtma is on a distinguished road
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.
hxaxtma is offline   Reply With Quote

Old   March 6, 2014, 10:14
Default
  #159
Senior Member
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 599
Rep Power: 20
bigphil will become famous soon enoughbigphil will become famous soon enough
Quote:
Originally Posted by hxaxtma View Post
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
bigphil is offline   Reply With Quote

Old   March 6, 2014, 10:31
Default
  #160
Member
 
Join Date: Jan 2014
Posts: 93
Rep Power: 4
hxaxtma is on a distinguished road
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
Attached Images
File Type: jpg couette_fsi.jpg (32.6 KB, 40 views)
hxaxtma is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
ESI-OpenCFD Releases OpenFOAM v2.2.0 opencfd OpenFOAM Announcements from ESI-OpenCFD 13 March 30, 2013 17:52
[Virtualization] OpenFOAM oriented tutorial on using VMware Player - support thread wyldckat OpenFOAM Installation 2 July 11, 2012 16:01
GPU Linear Solvers for OpenFOAM gocarts OpenFOAM Announcements from Other Sources 35 March 1, 2012 21:41
Cross-compiling OpenFOAM 1.7.0 on Linux for Windows 32 and 64bits with Mingw-w64 wyldckat OpenFOAM Announcements from Other Sources 3 September 8, 2010 06:25
OpenFOAM Debian packaging current status problems and TODOs oseen OpenFOAM Installation 9 August 26, 2007 13:50


All times are GMT -4. The time now is 21:33.