
[Sponsors] 
April 13, 2015, 02:06 
Mesh Motion is not consistent

#1 
Member
Pruthvi
Join Date: Feb 2014
Posts: 41
Rep Power: 12 
Hello guys I need to simulate an airfoil undergoing pitch + plunge as is done in the paper kinsey and dumas.
The pitch and plunge equations are as follows. Code:
H = h *sin(wt) theta = theta*cos(wt) Code:
dynamicFvMesh dynamicMotionSolverFvMesh; motionSolverLibs ("libsixDoFRigidBodyMotion.so"); solver sixDoFRigidBodyMotion; sixDoFRigidBodyMotionCoeffs { patches (wing); innerDistance 0.3; outerDistance 14; mass 10000; centreOfMass (0.3 0 0.125); momentOfInertia (1.958864357 3.920839234 10000); orientation ( 1 0 0 0 1 0 0 0 1 ); angularMomentum (0 0 1579.136704); //605.3357367); velocity (0 0 0); //0.1507964474 0); g (0 0 0); rhoName rhoInf; rhoInf 2.5; report on; constraints { yLine { sixDoFRigidBodyMotionConstraint line; centreOfRotation (0.3 0 0.125); direction (0 1 0); } zAxis { sixDoFRigidBodyMotionConstraint axis; axis (0 0 1); } } restraints { verticalSpring { sixDoFRigidBodyMotionRestraint linearSpring; anchor (0.3 1 0.125); refAttachmentPt (0.3 0 0.125); stiffness 227.3956854; damping 0; restLength 2; } axialSpring { sixDoFRigidBodyMotionRestraint linearAxialAngularSpring; axis (0.3 0 1); stiffness 227.3956854; damping 0; referenceOrientation $orientation; } } } Thanks, Pruthvi. Last edited by pruthvi1991; April 16, 2015 at 13:05. Reason: Added color to highlight mass, moment of inertia and stiffness 

April 14, 2015, 08:35 

#2 
Member
Thaw Tar
Join Date: Apr 2013
Location: Yangon, Myanmar
Posts: 35
Rep Power: 13 
How did you do this?
Does the mesh topology change? Did you use topology change dynamic mesh? 

April 16, 2015, 07:28 

#4 
Member
Alexander Bartel
Join Date: Feb 2015
Location: Germany
Posts: 97
Rep Power: 11 
Hi pruthvi,
Is w really the same for H an for theta? For me it seems that H has an higher frequency, and therefore the airfoil reaches highest/lowest point before being horizontal again. regards Alex 

April 16, 2015, 13:50 

#5 
Member
Pruthvi
Join Date: Feb 2014
Posts: 41
Rep Power: 12 
Hello Alex,
I think omega is the same for both. Check out the dynamicMeshdict file I posted. The mass, MOI and stiffness are highlighted in red. Thanks, Pruthvi. 

April 19, 2015, 12:00 

#6  
Senior Member
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 18 
I cannot understand what you are trying to model here.
Is it forced oscillation motion like T. Kinsey and G. Dumas paper 2008? OR you are trying to model a response problem? Because the equations of motion ,which are mentioned here, are for prescribed boundary motion, on the other hand you are using sixDoFRigidBodyMotion which calculate the displacement based on forces and moments every time step. Sorry, I think, I am missing something here. Bw, Hassan Quote:


April 20, 2015, 20:24 

#7 
Member
Pruthvi
Join Date: Feb 2014
Posts: 41
Rep Power: 12 
Hello Hassan,
I'm trying to do exactly what kinsey and dumas did. Forced oscillations. I'm not interested in the response of the system. I want to see the flow caused by the prescribed motion. Which other solver is apt for this kind of problem? I'm not very familiar with moving mesh techniques. Any help is appreciated in accurately describing the motion. Thanks, Pruthvi. 

April 21, 2015, 07:42 

#8 
Senior Member
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 18 
Hi Pruthvi,
There are two BC in OF for pointDisplacement which you can combine. I already used this approach before to test the dynamicMesh solvers in OF. I will describe it as a tutorial instead just uploading the files.They are oscillatingDisplacement and angularOscillatingDisplacement. The source code is located in Code:
$FOAM_SRC/fvMotionSolver/pointPatchFields/derived/ Code:
Field<vector>::operator=(amplitude_*sin(omega_*t.value())); Code:
scalar angle = angle0_ + amplitude_*sin(omega_*t.value()); vector axisHat = axis_/mag(axis_); vectorField p0Rel(p0_  origin_); vectorField::operator= ( p0Rel*(cos(angle)  1) + (axisHat ^ p0Rel*sin(angle)) + (axisHat & p0Rel)*(1  cos(angle))*axisHat ); (I will try to make it as simple as possible to be suitable for new OF users. Please excuse me if you found it boring but I am trying to make it as useful as possible.) 1  copy the "angularOscillatingDisplacement" to your run folder and rename it, lets call it for instance "twoDOscillatingDisplacement". Also replace "angularOscillatingDisplacement" with "twoDOscillatingDisplacement" inside the new twoDOscillatingDisplacementPointPatchVectorField.H and .C files. 2  Now you have to add two new private data to .H file Code:
// added for H vector amplitudeH_; //Heave motion amplitude (m) scalar omegaH_; // Heave motion frequency (rad/s) Code:
//first p0_(p.localPoints()), amplitudeH_(vector::zero), omegaH_(0.0) //second omega_(readScalar(dict.lookup("omega"))), amplitudeH_(dict.lookup("amplitudeH")), omegaH_(readScalar(dict.lookup("omegaH"))) //Third p0_(ptf.p0_, mapper), amplitudeH_(ptf.amplitudeH_), omegaH_(ptf.omegaH_) //fourth p0_(ptf.p0_), amplitudeH_(ptf.amplitudeH_), omegaH_(ptf.omegaH_) Code:
wing { type twoDOscillatingDisplacement; axis (0 0 1); // pitching axis origin (0.14944506531 0.01289078948 0.125); // elastic axis angle0 0; // mean theta amplitude 0.10; // theta amplitude (rad +ve C.C.W) omega 100; // frequency for theta (rad/s) value uniform (0 0 0); amplitudeH (0 0.1 0); // H amplitude omegaH 100; // H frequency } Code:
vectorField::operator= ( p0Rel*(cos(angle)  1) + (axisHat ^ p0Rel*sin(angle)) + (amplitudeH_*sin(omegaH_*t.value())) + (axisHat & p0Rel)*(1  cos(angle))*axisHat ); Code:
os.writeKeyword("amplitudeH") << amplitudeH_ << token::END_STATEMENT << nl; os.writeKeyword("omegaH") << omegaH_ << token::END_STATEMENT << nl; files Code:
twoDOscillatingDisplacementPointPatchVectorField.C LIB = $(FOAM_LIBBIN)/libmy2DOF options Code:
EXE_INC = \ I$(LIB_SRC)/triSurface/lnInclude \ I$(LIB_SRC)/meshTools/lnInclude \ I$(LIB_SRC)/dynamicMesh/lnInclude \ I$(LIB_SRC)/finiteVolume/lnInclude \ I$(LIB_SRC)/fileFormats/lnInclude LIB_LIBS = \ ltriSurface \ lmeshTools \ ldynamicMesh \ lfiniteVolume 8. finally try to complie it using wmake libso, then include the library in yourcase controlDict 9. modify dynamicMeshDict Code:
dynamicFvMesh dynamicMotionSolverFvMesh; motionSolverLibs ("libfvMotionSolvers.so"); solver displacementLaplacian; displacementLaplacianCoeffs { diffusivity inverseDistance (wing); // try different diffusivity } This is one approach to solve this problem but it isn't the only one. You can use NonConforming AMI Patches with solidBodyMotion like in the paper and in turbomachinery applications (relative reference frames). I am not familiar with this approach therefore I cannot give you any tips about it. Best wishes, Hassan 

May 1, 2015, 01:23 

#9 
Member
Thaw Tar
Join Date: Apr 2013
Location: Yangon, Myanmar
Posts: 35
Rep Power: 13 
Hi Pruthvi,
Which type of mesh motion solver do you use? I mean, how did you use topoChanging mesh and sixDofRigidBodyDisplacement together? Regards, Thaw Tar 

May 1, 2015, 01:26 

#10 
Member
Pruthvi
Join Date: Feb 2014
Posts: 41
Rep Power: 12 
@Hassan Many thanks for your solution. I will try to compile it and tell you about the results.
@Thaw I just followed the wing motion tutorial in pimpleDyMFoam. I changed the values of stiffness and mass to get the desired motion. 

May 1, 2015, 01:37 

#11 
Member
Thaw Tar
Join Date: Apr 2013
Location: Yangon, Myanmar
Posts: 35
Rep Power: 13 
Hi Pruthvi,
Thanks for your answer. Actually, I am doing 2D cylinder under vortex induced vibrations (viv) using sixDoFRigidBodyDisplacement. But if the motion amplitude is too high, the mesh distorts too much and the computation blows up. Thus, I wanted to use some kinds of topoChangerMesh together with sixDoFRigidBody motion solver. Regards thaw tar 

May 1, 2015, 01:47 

#12 
Member
Pruthvi
Join Date: Feb 2014
Posts: 41
Rep Power: 12 
Hello Thaw,
Whats the max amplitude that you may face? I could go upto 1 chord length without any computational issues. 

May 1, 2015, 04:09 

#13 
Senior Member
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 18 
@ThawTar, it is expected that the dynamic mesh solver cannot handle too high amplitude. But before trying topoChangMesh, try to use prescribed BC (oscillating displacement) with high amplitude. Solve the mesh only using moveDynamicMesh solver and try different diffusion methods. Also make sure that in your case, you are starting from convergent solution for fixed cylinder.
Bw Hassan 

May 1, 2015, 07:41 

#14  
Member
Pruthvi
Join Date: Feb 2014
Posts: 41
Rep Power: 12 
Hello Hassan,
I followed your instructions. They were pretty clear considering that I have little C++ knowledge. I couldn't compile it due to an error. Code:
wmakeLnInclude: linking include files to ./lnInclude make: *** No rule to make target `twoDOscillatingDisplacementPointPatchVectorField.dep', needed by `Make/linux64GccDPOpt/dependencies'. Stop. Quote:
Code:
#ifndef oscillatingDisplacementPointPatchVectorField_H #define oscillatingDisplacementPointPatchVectorField_H #include "fixedValuePointPatchField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { /**\ Class oscillatingDisplacementPointPatchVectorField Declaration \**/ class oscillatingDisplacementPointPatchVectorField : public fixedValuePointPatchField<vector> { // Private data vector amplitude_; scalar omega_; // added for H vector amplitudeH_; //Heave motion amplitude (m) scalar omegaH_; // Heave motion frequency (rad/s) public: // Runtime type information TypeName("oscillatingDisplacement"); // Constructors // Construct from patch and internal field oscillatingDisplacementPointPatchVectorField ( const pointPatch&, const DimensionedField<vector, pointMesh>& ); // Construct from patch, internal field and dictionary oscillatingDisplacementPointPatchVectorField ( const pointPatch&, const DimensionedField<vector, pointMesh>&, const dictionary& ); // Construct by mapping given patchField<vector> onto a new patch oscillatingDisplacementPointPatchVectorField ( const oscillatingDisplacementPointPatchVectorField&, const pointPatch&, const DimensionedField<vector, pointMesh>&, const pointPatchFieldMapper& ); // Construct and return a clone virtual autoPtr<pointPatchField<vector> > clone() const { return autoPtr<pointPatchField<vector> > ( new oscillatingDisplacementPointPatchVectorField ( *this ) ); } // Construct as copy setting internal field reference oscillatingDisplacementPointPatchVectorField ( const oscillatingDisplacementPointPatchVectorField&, const DimensionedField<vector, pointMesh>& ); // Construct and return a clone setting internal field reference virtual autoPtr<pointPatchField<vector> > clone ( const DimensionedField<vector, pointMesh>& iF ) const { return autoPtr<pointPatchField<vector> > ( new oscillatingDisplacementPointPatchVectorField ( *this, iF ) ); } // Member functions // Evaluation functions // Update the coefficients associated with the patch field virtual void updateCoeffs(); // Write virtual void write(Ostream&) const; }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif // ************************************************************************* // Thanks , Pruthvi. 

May 1, 2015, 08:10 

#15 
Senior Member
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 18 
You have to change each instance in .C and .H from oscillatingDisplacement to twoDoscillatingDisplacement. Check the Github, hpfully the changes will appear for you because I am not familiar with it
Computational methods for fluid dynamics has a small section about dynamic meshes. Also review Jasak's papers about dynamicmesh handling in OpenFOAM. About the diffusion coefficient effect there is a report from Charmers about it http://www.tfd.chalmers.se/~hani/kur...AMrapport.pdf Please note since OF2.3.0 there are major changes in the solid body motion library check http://www.openfoam.org/version2.3.0/meshmotion.php 

May 1, 2015, 09:12 

#16 
Senior Member
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 18 
I noticed extra few points needs to be modified in your code. I will search for my old code and upload it for you.


May 1, 2015, 10:17 

#17 
Senior Member
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 18 
Here is the link
https://github.com/HIKassem/twoDOscillatingDisplacement it is my first time to use github, so it could be not working or I messed up something. Bw, Hassan 

May 10, 2015, 07:45 

#18 
Member
Pruthvi
Join Date: Feb 2014
Posts: 41
Rep Power: 12 
Hi Hassan,
Your code worked fine. I was able to run a couple of simulations. Nice implementation. I ran into some SigFpe error. Now debugging. I will let you know if I can debug. Nice picking up github . Thanks, Pruthvi. 

June 19, 2015, 21:42 
Courant number suddenly blows up

#19 
Member
Pruthvi
Join Date: Feb 2014
Posts: 41
Rep Power: 12 
Hello Mr.Hassan,
I modified the twoDOscillaitngMotion.C file so that I can introduce a phase difference of pi/2 between the pitching and plunging motions. The simulation keeps blowing up. I used checkMesh and found that the quality was degrading at the cells generated by layer addition. To rectify I switched off layer addition and tried again. This has eleminated most of the mesh check errors but the courant number is still blowing up. All the files are here > https://github.com/pruthvi1991/solve...r/Kinsey_Dumas Here is a checkMesh log when the Courant number blows up but the mesh quality is OK Code:
Time = 9 Mesh stats points: 74594 internal points: 0 faces: 146827 internal faces: 73328 cells: 36510 faces per cell: 6.029991783 boundary patches: 6 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 35543 prisms: 11 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 0 polyhedra: 956 Breakdown of polyhedra by number of faces: faces number of cells 7 806 8 150 Checking topology... Boundary definition OK. Cell to face addressing OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces... Patch Faces Points Surface topology topAndBottom 160 324 ok (nonclosed singly connected) inlet 104 210 ok (nonclosed singly connected) outlet 104 210 ok (nonclosed singly connected) wing 111 222 ok (nonclosed singly connected) front 36510 37297 ok (nonclosed singly connected) back 36510 37297 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (35 35 0) (40 35 0.5) Mesh (nonempty, nonwedge) directions (1 1 0) Mesh (nonempty) directions (1 1 0) All edges aligned with or perpendicular to nonempty directions. Boundary openness (6.627104753e19 4.888447076e22 1.263269968e14) OK. Max cell openness = 2.199932693e16 OK. Max aspect ratio = 3.839901044 OK. Minimum face area = 4.396716366e05. Maximum face area = 1.945327376. Face area magnitudes OK. Min volume = 2.198358183e05. Max volume = 0.9726636881. Total volume = 2624.959185. Cell volumes OK. Mesh nonorthogonality Max: 64.3969137 average: 9.591723067 Nonorthogonality check OK. Face pyramids OK. Max skewness = 3.007630562 OK. Coupled point location match (average 0) OK. Mesh OK. End Code:
Create time Create polyMesh for time = 25 Time = 25 Mesh stats points: 74594 internal points: 0 faces: 146827 internal faces: 73328 cells: 36510 faces per cell: 6.029991783 boundary patches: 6 point zones: 0 face zones: 0 cell zones: 0 Overall number of cells of each type: hexahedra: 35543 prisms: 11 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 0 polyhedra: 956 Breakdown of polyhedra by number of faces: faces number of cells 7 806 8 150 Checking topology... Boundary definition OK. Cell to face addressing OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces... Patch Faces Points Surface topology topAndBottom 160 324 ok (nonclosed singly connected) inlet 104 210 ok (nonclosed singly connected) outlet 104 210 ok (nonclosed singly connected) wing 111 222 ok (nonclosed singly connected) front 36510 37297 ok (nonclosed singly connected) back 36510 37297 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (35 35 0) (40 35 0.5) Mesh (nonempty, nonwedge) directions (1 1 0) Mesh (nonempty) directions (1 1 0) All edges aligned with or perpendicular to nonempty directions. Boundary openness (6.61854997e19 9.776894152e22 6.381205446e15) OK. Max cell openness = 2.354408509e16 OK. Max aspect ratio = 3.293216054 OK. Minimum face area = 5.64186484e05. Maximum face area = 1.952868252. Face area magnitudes OK. Min volume = 2.82093242e05. Max volume = 0.9764341259. Total volume = 2624.959185. Cell volumes OK. Mesh nonorthogonality Max: 52.72379189 average: 8.614508389 Nonorthogonality check OK. Face pyramids OK. ***Max skewness = 4.529272567, 1 highly skew faces detected which may impair the quality of the results <<Writing 1 skew faces to set skewFaces Coupled point location match (average 0) OK. Failed 1 mesh checks. En Can this problem be solved by trying a different diffusivity? I tried quadratic already. Last edited by pruthvi1991; June 20, 2015 at 18:32. 

June 20, 2015, 06:42 

#20 
Senior Member
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 18 
I have few questions,
Could you please tell me which case on github you are using? Are you using physical solver of moveDynamicMesh? is the solver blowup even when you use small amplitude values? Bw, Hassan 

Tags 
flapping wing simulation, pimpledymfoam, sixdofrigidbodymotion 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
how to set periodic boundary conditions  Ganesh  FLUENT  15  November 18, 2020 06:09 
Update of the variables after dynamic mesh motion.  gtg258f  OpenFOAM Programming & Development  9  January 18, 2014 10:08 
Mesh motion with Translation & Rotation  Doginal  CFX  2  January 12, 2014 06:21 
3D Hybrid Mesh Errors  DarrenC  ANSYS Meshing & Geometry  11  August 5, 2013 06:42 
Dynamic moving mesh  PeiYing Hsieh (Hsieh)  OpenFOAM Running, Solving & CFD  64  June 7, 2012 10:04 