
[Sponsors] 
tractionDisplacement BC for solidDisplacementFoam with Body Forces 

LinkBack  Thread Tools  Search this Thread  Display Modes 
April 13, 2012, 21:39 
tractionDisplacement BC for solidDisplacementFoam with Body Forces

#1 
Senior Member
Hisham Elsafti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 257
Blog Entries: 10
Rep Power: 17 
Hi Foamers,
I need to add body forces to the solidDisplacementFoam solver. So it goes something similar to: Code:
fvVectorMatrix DEqn ( fvm::d2dt2(D) == fvm::laplacian(2*mu + lambda, D, "laplacian(DD,D)") + bodyForces + divSigmaExp ); Code:
gradient() = ( (traction_ + pressure_*n)/rho + twoMuLambda*fvPatchField<vector>::snGrad()  (n & sigmaD) )/twoMuLambda; if (thermalStress) { ... gradient() += n*threeKalpha*T/twoMuLambda; } Best regards, Hisham El Safti 

April 14, 2012, 02:40 

#2  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Quote:
Body forces act on the volume of each cell and are source terms so they do not have boundary conditions. So to add a body force term, just add 'rho*g' to the momentum equation, where g is a dimensionedVector 9.81 in the y direction. Philip 

April 14, 2012, 06:24 

#3  
Senior Member
Hisham Elsafti
Join Date: Apr 2011
Location: Braunschweig, Germany
Posts: 257
Blog Entries: 10
Rep Power: 17 
Hi Philip
Quote:
I also have some other explicit terms that come from another variable (pore fluid velocity). They are all explicit terms (as the case of the "T" term in the DEqn in solidDisplacementFoam.C) Code:
if (thermalStress) { const volScalarField& T = Tptr(); DEqn += fvc::grad(threeKalpha*T); } Best regards, Hisham 

April 14, 2012, 06:49 

#4 
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Hi Hisham,
The reason there is a Tterm in the boundary condition is because T appears in the constitutive equation i.e. Hooke's law for a thermal elastic solid is: sigma =2*mu*e +lambda*I*tr(e)  (2*mu + 3*lambda)*alpha*(T  Tref)*I Therefore the traction on a boundary face is Trac = n . sigma, and then you sub in Hooke's law and rearrange this equation to get the equivalent (implicit) gradient on the boundary, which is what you set in the traction BC. So the traction BC will contain a Tterm. But it will not contain any body force contributions. If you add rho*g to the momentum equation then it will work, I ran a quick test case and everything is fine on the traction boundaries. Philip 

April 29, 2013, 13:26 

#6 
Member
Samer
Join Date: Jan 2013
Posts: 31
Rep Power: 13 
Hi Philip
using the same methodology, if simulation fluidstructure interaction problem, and i want to transfer the temperature computed at the wall of a solid, as boundary condition for the fluid part (where there is a common patch between fluid and solid) ?? Best regards 

April 30, 2013, 16:27 

#7  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Quote:
yes you can follow the same DirichletNeumann coupling procedure to couple the temperature between the regions. Philip 

November 17, 2017, 09:27 

#8 
Member
Join Date: Sep 2016
Posts: 63
Rep Power: 9 
Hello Philip,
Thank you very much first of all for your solvers in solidMechanics incorporated in foamextend! I want to do a stress analysis of a fan with a certain rotation speed. The material remain elastic and I assume Newtonian properties. I wonder which solver can implement the centrifugal force please? I am new for your solvers. I checked through your solvers, but did't figure it out. Thank you very much in advance! Best regards, sitajeje 

January 9, 2018, 13:20 

#9  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Quote:
To include a body force, you just need to add this term to the momentum equation in the solver. For example, in elasticSolidFoam the momentum equation would become (on line 93 of elasticSolidFoam.C): Code:
fvVectorMatrix UEqn ( rho*fvm::d2dt2(U) == fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)") + divSigmaExp + rho*g + rho*magSqr(v)/r // centrifugal force per unit volume ); Philip 

January 9, 2018, 13:38 

#10 
Member
Join Date: Sep 2016
Posts: 63
Rep Power: 9 
Thank you very much Philip!
I get it done by myself, and my programming is improved a little. Thank you all the same! I have another question please, and I cannot figure it out until now. It would be very helpful if you could give me a suggestion. I wonder how to map the pressure field of the fan blades from results of pimpleDyMFoam to solidDisplacementFoam please? I did mapFields, but "p" is mapped to "p", and didn't effect. How can I map "p" to the "p" in the "tractionDisplacement" boundary condition in the "D" file please? I want to calculate the stress induced by the aerodynamic pressure. Thank you very much in advance! Best regards, sitajeje 

January 25, 2018, 07:41 

#11 
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Hi sitajeje,
mapFields is used to map from one domain to another domain that overlaps the original domain in some way e.g. see cavityClipped case in the documentation; however, it does not map (as far as I know) from a patch on one domain to a patch on another domain, where the patches overlap but the domains do not. So in this case, you need to use the functionality provided by classes like patchToPatchInterpolation or GGIInterpolation (foamextend) or AMIInterpolation (Foundation and ESI forks). If you send me a PM with your email, I have some code that may work for you. Philip 

February 15, 2018, 15:19 

#12  
New Member
Florida
Join Date: Feb 2018
Location: Florida
Posts: 12
Rep Power: 8 
Quote:
I was trying to introduce the body force as a centrifugal force per unit volume in a beam like the example you gave. However I am having trouble introducing the radius vector by reading it from the mesh(which is the z coordinate of the cell center of each cell in my mesh). Please help. Thanks. 

February 16, 2018, 04:44 

#13 
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Hello shinjanghosh,
What is your problem? You will need to give more information if you want someone to help you. Also, looking back on the body force term I suggested above, the dimensions are wrong: it should be a vector not a scalar; the following makes more sense: Code:
+ rho*(magSqr(v)/mag(R))*(R/mag(R)) Code:
// Points on the axis const vector pointOnAxis1(0, 0, 0); const vector pointOnAxis2(1, 0, 0); // Axis unit direction vector const vector axisDir = (pointOnAxis1  pointOnAxis2)/mag(pointOnAxis1  pointOnAxis2); // Radial vector field: it only has a component in the radial direction const volVectorField R ( "R", (I  sqr(axisDir)) & (mesh.C()  pointOnAxis1) ) 

February 16, 2018, 13:51 

#14 
New Member
Florida
Join Date: Feb 2018
Location: Florida
Posts: 12
Rep Power: 8 
Hello Philip,
Thanks for the reply. Sorry for the inconvenience. My problem : Calculation of solid stresses in a hollow beam with a circular cross section(basically a pipe with a finite wall thickness) which rotates at a rpm of 3600. the length of the pipe is along the 'zaxis' and it's axis of rotation is the x axis. I figured out the boundary conditions (fixed end cantilever beam) and did a simple end load condition simulation solution (without any rotation) using solidDisplacementFoam. Then I started having problems adding the body forces inside the vector equation when I decided to include the rotation. I tried hardcoding the body force in the solver itself like you suggested in one of your previous replies. Here's where I had issues : 1.) Reading the cellcenter values for calculation of 'R' vector. In my case it is the z component of the cell center. 2.) Specifying it as a vector. Thanks, Shinjan 

February 19, 2018, 07:23 

#15 
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Hi Shinjan,
As the xaxis is your axis of rotation, the updated code I suggested in my previous post should work for you. If you have issues with compilation of running, you can post them here. Philip 

February 19, 2018, 13:09 

#16 
New Member
Florida
Join Date: Feb 2018
Location: Florida
Posts: 12
Rep Power: 8 
Hey Philip,
The compilation didn't give any errrors, but I encountered the following errors while running the solver :  > FOAM FATAL ERROR: LHS and RHS of  have different dimensions dimensions : [0 1 0 0 0 0 0]  [0 0 0 0 0 0 0] From function Foam::dimensionSet Foam:perator(const Foam::dimensionSet&, const Foam::dimensionSet&) in file dimensionSet/dimensionSet.C at line 521. FOAM aborting Generating stack trace...  From what I understood, It has to be a dimension mismatch regarding subtracting a vector which is dimensionless, from one which has length as it's dimension. I think the error came from here : const volVectorField R ( "R", (I  sqr(axisDir)) & (mesh.C()  pointOnAxis1) ); 

February 20, 2018, 06:56 

#17 
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Hi Shinjan,
The pointOnAxis1, pointOnAxis2 and aixsDir should be dimensionedVector instead of vector: try make the changes below: Code:
// Points on the axis const dimensionedVector pointOnAxis1("pointOnAxis1", dimLength, vector(0, 0, 0)); const dimensionedVector pointOnAxis2("pointOnAxis2", dimLength, vector(1, 0, 0)); // Axis unit direction vector const dimensionedVector axisDir = (pointOnAxis1  pointOnAxis2)/mag(pointOnAxis1  pointOnAxis2); // Radial vector field: it only has a component in the radial direction const volVectorField R ( "R", (I  sqr(axisDir)) & (mesh.C()  pointOnAxis1) ) 

March 27, 2018, 13:28 

#18  
New Member
Dmitry
Join Date: Dec 2016
Location: Leoben, Austria
Posts: 1
Rep Power: 0 
Quote:
I'm also as sitajeje interested in mapping "p" from one region (calculated in fluid) to another (solid). I'm doing a steadystate heat transfer simulation using one of multiregion solvers (I've modified one from chtMultiRegionSimpleFoam for my purpose) and now I'm interested in stresses in the solid region. I've copied the last time step of the solid region from my multiregion simulation, added a "D" file and I wanted to simply run a oneregion simulation on it using the "solidDisplacementFoam" solver. My question is how I correctly set pressure in the "tractionDisplacement" boundary condition which is, from my understanding, the pressure field on the interface between the fluid and solid regions in my previous simulation. However, the boundary condition for "p" in fluid was set as "type zeroGradient" (initially I expected to see a nonuniform list which I wanted to simply copy in the "D file" for the pressure field instead of "uniform 0"). I'm wondering whether Philip's answer (using patchToPatchInterpolation or GGIInterpolation or AMIInterpolation) is the right direction for me or there is a simpler solution. I'm using OpenFOAM 2.4.0. Best regards, Dmitry 

March 29, 2018, 04:32 

#19  
Super Moderator
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,089
Rep Power: 34 
Quote:
From my understanding of your case, it should be possible to explicitly interpolate fields using patchToPatchInterpolation/GGIInterpolation/AMIInterpolation. However, be careful with "pressure" in the tractionDisplacement specification: here pressure is the surface normal component of the wall traction vector, this is not the same (in general) as the hydrostatic pressure (i.e. "p" in the fluid solver). So you should reconstruct the total stress/traction vector on the fluid side and interpolate this across to the solid side (tractionDisplacement). Philip 

Tags 
bodyforces, soliddisplacementfoam, tractiondisplacement 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
DPM body force to express electric forces  adam14qin  FLUENT  24  August 1, 2013 13:34 
Body forces on hovering airfoil  sanmysterio  Main CFD Forum  1  July 13, 2010 03:14 
Fan modelling using body forces  Fred  Siemens  2  June 13, 2007 08:56 
strong body forces  Ankan  Main CFD Forum  0  August 2, 2005 02:13 
Body Forces  Thomas P. Abraham  Main CFD Forum  1  April 21, 1999 21:24 