
[Sponsors] 
[solidMechanics] Support thread for "Solid Mechanics Solvers added to OpenFOAM Extend" 

LinkBack  Thread Tools  Search this Thread  Display Modes 
October 3, 2013, 16:49 

#121  
Senior Member
Daniel
Join Date: Mar 2013
Location: Noshahr, Iran
Posts: 348
Rep Power: 21 
Quote:
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 Last edited by Daniel_Khazaei; October 3, 2013 at 18:14. 

October 4, 2013, 05:41 

#122 
New Member
Vincent
Join Date: Aug 2013
Location: Singapore
Posts: 7
Rep Power: 12 
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 

October 9, 2013, 02:38 

#123 
New Member
Vincent
Join Date: Aug 2013
Location: Singapore
Posts: 7
Rep Power: 12 
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/vin82/OpenFOAM/OpenFOAM1.6ext/ThirdParty/rpmBuild/BUILD/ParaView3.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 

October 29, 2013, 08:01 

#124 
New Member
Join Date: Oct 2013
Posts: 5
Rep Power: 12 
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. Code:
Time = 0.065, iteration: 2 Current fsi underrelaxation factor: 0.01 Maximal accumulated displacement of interface points: 0.0084761 Courant Number mean: 0.0443861 max: 0.522711 velocity magnitude: 1.31938 DILUPBiCG: Solving for Ux, Initial residual = 4.75368e05, Final residual = 3.94291e07, No Iterations 1 DILUPBiCG: Solving for Uy, Initial residual = 0.000110373, Final residual = 8.77263e07, No Iterations 1 GAMG: Solving for p, Initial residual = 0.000543828, Final residual = 8.25319e07, No Iterations 30 GAMG: Solving for p, Initial residual = 0.000102911, Final residual = 9.47413e07, No Iterations 10 time step continuity errors : sum local = 3.23041e10, global = 4.55244e11, cumulative = 5.6536e10 GAMG: Solving for p, Initial residual = 0.000131358, Final residual = 7.9384e07, No Iterations 24 GAMG: Solving for p, Initial residual = 1.99396e05, Final residual = 7.79756e07, No Iterations 6 time step continuity errors : sum local = 2.65822e10, global = 2.4541e11, cumulative = 5.89901e10 Setting traction on solid patch Total traction force = (0.00431445 7.4879e06 3.43629e22) Solving for DU, Initial residual = 0.00589588, Final residual = 9.88846e07, No outer iterations 41 Current fsi residual norm: 0.981515 Time = 0.065, iteration: 3 Current fsi underrelaxation factor (Aitken): 0.362946 Maximal accumulated displacement of interface points: 0.00876687 [0] [0] [0] > FOAM FATAL ERROR: [0] face 41 area does not match neighbour by 0.0103705%  possible face ordering problem. patch: procBoundary0to1 my area:0.000239443 neighbour area: 0.000239418 matching tolerance: 0.0001 Mesh face: 42167 vertices: 4((0.574801 0.214706 0.025334) (0.574854 0.20998 0.025334) (0.574854 0.20998 0.025334) (0.574801 0.214706 0.025334)) Rerun with processor debug flag set for more information. [0] [0] From function processorPolyPatch::calcGeometry() [0] in file [1] [1] [1] > FOAM FATAL ERROR: [1] face 41 area does not match neighbour by 0.0103705%  possible face ordering problem. patch: procBoundary1to0 my area:0.000239418 neighbour area: 0.000239443 matching tolerance: 0.0001 Mesh face: 42054 vertices: 4((0.574801 0.214706 0.025334) (0.574801 0.214706 0.025334) (0.574854 0.209981 0.025334) (0.574854 0.209981 0.025334)) Rerun with processor debug flag set for more information. [1] [1] From function processorPolyPatch::calcGeometry() [1] in file meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C at line 216. [1] FOAM parallel run exiting [1] meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C at line 216. [0] FOAM parallel run exiting [0]  MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD with errorcode 1. NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes. You may or may not see output from other processes, depending on exactly when Open MPI kills them.  Michael 

December 11, 2013, 07:49 

#125 
New Member
Join Date: Mar 2013
Posts: 17
Rep Power: 13 
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 

December 11, 2013, 07:53 

#126  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Quote:
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: Code:
if (runTime.timeIndex() == 1) { // Add gravity increment on first time step DUEqn = rho*dimensionedVector("gravity", dimLength/(dimTime/dimTime), vector(0,9.81,0)); } Last edited by bigphil; December 11, 2013 at 08:00. Reason: Adding example code 

December 11, 2013, 08:29 

#127  
New Member
Join Date: Mar 2013
Posts: 17
Rep Power: 13 
Quote:
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 fluidstructure 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 

December 11, 2013, 08:36 

#128  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Quote:
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. Quote:
Philip 

December 11, 2013, 09:22 

#129  
New Member
Join Date: Mar 2013
Posts: 17
Rep Power: 13 
Quote:
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. flappingConsole_velocity2.jpg 

December 11, 2013, 09:30 

#130  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Quote:
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:
Philip 

December 12, 2013, 08:27 

#131  
New Member
Join Date: Mar 2013
Posts: 17
Rep Power: 13 
Quote:
I'm currently looking into the structure solver used by the fluidstructure 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. 

December 12, 2013, 09:46 

#132  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
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: OpenFOAM1.6ext/src/finiteVolume/finiteVolume/d2dt2Schemes/backwardD2dt2Scheme Alternatively a smaller timestep should help. Quote:
Philip 

December 14, 2013, 06:09 

#133 
Member
Ben
Join Date: Sep 2012
Posts: 40
Rep Power: 13 
I have solved my previous problem but have a new one. I used the github for the official OpenFOAM 2.2x (I am using OpenFOAM2.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:
Code:
/**\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 2.2.2   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ Build : 2.2.29240f8b967db Exec : elasticSolidFoam Date : Dec 15 2013 Time : 12:37:57 Host : "skVirtualBox" PID : 4790 Case : /home/sk/OpenFOAM/sk2.2.2/run/solidMechanics/tutorials/elasticSolidFoam/plateHole nProcs : 1 sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE). fileModificationChecking : Monitoring runtime modified files using timeStampMaster allowSystemOperations : Disallowing usersupplied system call operations // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Create time Create mesh for time = 0 Reading field U Patch hole Traction boundary field: U > FOAM FATAL IO ERROR: Unknown grad scheme extendedLeastSquares Valid grad schemes are : 9 ( Gauss cellLimited cellMDLimited edgeCellsLeastSquares faceLimited faceMDLimited fourth leastSquares pointCellsLeastSquares ) file: /home/sk/OpenFOAM/sk2.2.2/run/solidMechanics/tutorials/elasticSolidFoam/plateHole/system/fvSchemes.gradSchemes.grad(U) at line 26. From function gradScheme<Type>::New(const fvMesh& mesh, Istream& schemeData) in file /opt/openfoam222/src/finiteVolume/lnInclude/gradScheme.C at line 72. FOAM exiting Thank you. Last edited by ben1793; December 15, 2013 at 12:15. 

December 21, 2013, 17:12 

#134 
Senior Member
Daniel
Join Date: Mar 2013
Location: Noshahr, Iran
Posts: 348
Rep Power: 21 
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 (atleast 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: Code:
do { outerCorr++; # include "setInterfaceDisplacement.H" # include "moveFluidMesh.H" # include "solveFluid.H" if (runTime.timeIndex() >= 4) { #include "setInterfaceForce.H" } # include "solveSolid.H" # include "calcFsiResidual.H" } Is there any side effects (skipping a few time steps) that I need to know? regards 

January 15, 2014, 20:48 

#135  
Member
Eric Bryant
Join Date: Sep 2013
Location: Texas
Posts: 44
Rep Power: 12 
@ Daniel_Khazaei
Quote:
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. 

January 28, 2014, 05:59 

#136  
New Member
Join Date: Mar 2013
Posts: 17
Rep Power: 13 
Quote:
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 fluidstructure 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. Zeroth order behaviour is shown for the total displacement U, whereas the first order scheme shows the correct order. Any ideas? Many thanks. Best, David 

January 28, 2014, 06:15 

#137  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Quote:
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 builtin backward/Euler schemes like elasticNonLinULSolidFoam. Best regards, Philip 

January 28, 2014, 08:03 

#138  
New Member
Join Date: Mar 2013
Posts: 17
Rep Power: 13 
Quote:
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. 

January 28, 2014, 08:47 

#139 
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
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>& ); So, for now, you can try out the backward d2tdt2 scheme by changing: Code:
fvm::d2dt2(rho,DU) Code:
rho*fvm::d2dt2(DU) It is probably straightforward to implement the second two functions if you get an urge Best regards, Philip 

January 28, 2014, 16:45 
Regarding elasticNonLinTLSolidFoam (cantilever problem)

#140 
New Member
Join Date: Oct 2011
Posts: 28
Rep Power: 14 
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 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
GPU Linear Solvers for OpenFOAM  gocarts  OpenFOAM Announcements from Other Sources  37  August 17, 2022 14:22 
[Virtualization] OpenFOAM oriented tutorial on using VMware Player  support thread  wyldckat  OpenFOAM Installation  2  July 11, 2012 16:01 
New OpenFOAM Forum Structure  jola  OpenFOAM  2  October 19, 2011 06:55 
Crosscompiling OpenFOAM 1.7.0 on Linux for Windows 32 and 64bits with Mingww64  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 