[Hisham]

Hi Foamers,

I need to add body forces to the solidDisplacementFoam solver. So it goes something similar to:

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:

gradient() =
    (
        (traction_ + pressure_*n)/rho
      + twoMuLambda*fvPatchField<vector>::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

[bigphil]

Originally Posted by Hisham

Hi Foamers,

I need to add body forces to the solidDisplacementFoam solver. So it goes something similar to:

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:

gradient() =
    (
        (traction_ + pressure_*n)/rho
      + twoMuLambda*fvPatchField<vector>::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

[Hisham]

Hi Philip

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)

 if (thermalStress)
                {
                    const volScalarField& T = Tptr();
                    DEqn += fvc::grad(threeKalpha*T);
                }

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

[bigphil]

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

[Hisham]

Hi Philip

Thanks a lot for your help! Now it is clear to me

Best regards
Hisham

[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

[bigphil]

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

[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

[bigphil]

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

            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

[sitajeje]

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

[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

[shinjanghosh]

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

            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.

[bigphil]

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:

        + rho*(magSqr(v)/mag(R))*(R/mag(R))

where R is the radial volVectorField, which could be calculated as:

// 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

[shinjanghosh]

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

[bigphil]

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

[shinjanghosh]

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)
);

[bigphil]

Hi Shinjan,

The pointOnAxis1, pointOnAxis2 and aixsDir should be dimensionedVector instead of vector: try make the changes below:

// 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

[parnumeric]

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

[bigphil]

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