# tractionDisplacement BC for solidDisplacementFoam with Body Forces

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

 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: 16 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 );``` Where the bodyForces vectorField represents accelerations (e.g. gravity). Nevertheless, I have not found a reference on how to add a term to gradient() in tractionDisplacement boundary condition. I mean that for example a term for strain resulting from thermal stresses is added to the gradient() in tractionDisplacementFvPatchVectorField.C like: Code: ```gradient() = ( (traction_ + pressure_*n)/rho + twoMuLambda*fvPatchField::snGrad() - (n & sigmaD) )/twoMuLambda; if (thermalStress) { ... gradient() += n*threeKalpha*T/twoMuLambda; }``` I really appreciate any ideas or directions to references (or even how to do it ) Best regards, Hisham El Safti

April 14, 2012, 02:40
#2
Super Moderator

Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,078
Rep Power: 33
Quote:
 Originally Posted by Hisham 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 );``` Where the bodyForces vectorField represents accelerations (e.g. gravity). Nevertheless, I have not found a reference on how to add a term to gradient() in tractionDisplacement boundary condition. I mean that for example a term for strain resulting from thermal stresses is added to the gradient() in tractionDisplacementFvPatchVectorField.C like: Code: ```gradient() = ( (traction_ + pressure_*n)/rho + twoMuLambda*fvPatchField::snGrad() - (n & sigmaD) )/twoMuLambda; if (thermalStress) { ... gradient() += n*threeKalpha*T/twoMuLambda; }``` I really appreciate any ideas or directions to references (or even how to do it ) Best regards, Hisham El Safti
Hi Hisham,

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: 16
Hi Philip

Quote:
 Originally Posted by bigphil 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
I agree with you that body forces do not need a BC for themselves but don't they contribute to displacement gradient at the boundary? I know that the boundary values can be retrieved from cells next to it.

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();
}```
So is the existence of a "T"-term added to the gradient() [at D's tractionDisplacement BC] due to the existence of a "T" BC or because "T" contributes to the displacement gradient?

Best regards,
Hisham

 April 14, 2012, 06:49 #4 Super Moderator     Philip Cardiff Join Date: Mar 2009 Location: Dublin, Ireland Posts: 1,078 Rep Power: 33 Hi Hisham, The reason there is a T-term 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 T-term. 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 14, 2012, 07:03 #5 Senior Member     Hisham Elsafti Join Date: Apr 2011 Location: Braunschweig, Germany Posts: 257 Blog Entries: 10 Rep Power: 16 Hi Philip Thanks a lot for your help! Now it is clear to me Best regards Hisham

 April 29, 2013, 13:26 #6 Member   Samer Join Date: Jan 2013 Posts: 31 Rep Power: 12 Hi Philip using the same methodology, if simulation fluid-structure 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,078
Rep Power: 33
Quote:
 Originally Posted by SamerAli Hi Philip using the same methodology, if simulation fluid-structure 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
Hi Samer,

yes you can follow the same Dirichlet-Neumann 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 foam-extend! 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,078
Rep Power: 33
Quote:
 Originally Posted by sitajeje Hello Philip, Thank you very much first of all for your solvers in solidMechanics incorporated in foam-extend! 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
Hi sitajeje,

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
);```
where you would need to define your velocity volVectorField "v" and the radial distance field "r". You could read these fields from the case.

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,078 Rep Power: 33 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 (foam-extend) 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: 7
Quote:
 Originally Posted by bigphil Hi sitajeje, 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 );``` where you would need to define your velocity volVectorField "v" and the radial distance field "r". You could read these fields from the case. Philip
Hello Philip,
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,078 Rep Power: 33 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))` where R is the radial volVectorField, which could be calculated as: 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) )``` Philip

 February 16, 2018, 13:51 #14 New Member   Florida Join Date: Feb 2018 Location: Florida Posts: 12 Rep Power: 7 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 'z-axis' 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 hard-coding 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 cell-center 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,078 Rep Power: 33 Hi Shinjan, As the x-axis 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: 7 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,078 Rep Power: 33 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) )``` Philip

March 27, 2018, 13:28
#18
New Member

Dmitry
Join Date: Dec 2016
Location: Leoben, Austria
Posts: 1
Rep Power: 0
Quote:
 Originally Posted by bigphil 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 (foam-extend) 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
Hello Philip, hi sitajeje and Foamers,

I'm also as sitajeje interested in mapping "p" from one region (calculated in fluid) to another (solid). I'm doing a steady-state 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 one-region 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 non-uniform 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,078
Rep Power: 33
Quote:
 Originally Posted by parnumeric Hello Philip, hi sitajeje and Foamers, I'm also as sitajeje interested in mapping "p" from one region (calculated in fluid) to another (solid). I'm doing a steady-state 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 one-region 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 non-uniform 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
Hi Dmitry,

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