[hsieh]

Hi,

I completed a case using interFoam. How can I compute shear stress of the 3D domain based on the transient data?

Thanks!

Pei

[hartinger]

// you get the stress tensor with:

volTensorField stressTensor = - mu * (fvc::grad(U) + (fvc::grad(U)).T());
// there might be a better way describing that, if haven't used the stressTensor directly

//if you're only interested in the shearRate, for shear-thinnning for example:

volScalarField shearRate = mag(fvc::grad(U));

for the wall shear stress check the tool 'wallShearStress'

regards
markus

[hsieh]

Thanks a lot Markus!

Pei

[santos]

Hello Pei

The tool 'wallShearStress' didnt work for me (it kept asking me about turbulence model, but I was using DNS), so what I did was:

- Calculate U gradient at wall (using wallGradU utility);

- Then calculate the magnitude of wall shear stress by using the following:
volScalarField shear
(
IOobject
(
"shear",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
8.9e-4*sqrt(sqr(wallGradU.component(vector::X)) + sqr(wallGradU.component(vector::Z)))
);
shear.write();

X and Z components were relevant in my case, your situation may be different. 8.9e-4 is the dynamic viscosity for water at room temperature, this way the shear stress comes in Pa.

Hope that this can help you!

Regards,
José Santos

[hsieh]

Hi, José Santos,

Thanks for the reply!

Do you have any idea if I would like to compute the tangential direction of wall shear stress? Because my wall is curved, so, the X and Z components may not necessarily tangent.

Pei

[liu]

dot product with n (normal direction)

[hsieh]

Hi, Xiaofeng:

Thanks for the answer. I did try it, but, I got a message said that n is not defined. What will I have to include to get n in post-processing?

Pei

[liu]

I remember there are some discussion on this issure. Search "wall shear stress" may do the work.

[hsieh]

Hi, Xiaofeng,

I did do a search on "wall shear stress". I have also looked at the source code of wallShearStress under postprocessing/wall. It is not obvious to me how to get n.

Pei

[philippose]

Hi,
I just happened to browse through this dialogue regarding calculating wall shear stresses.

In the simple case, if you take laminar flow, the wall shear stress is calculated using the equation:

Tau (shear stress) = mu * (dU/dy)

In this equation, "mu" is the dynamic viscosity of the medium, and "(dU/dy)" is the change in velocity "U" with change in "y", which is the distance from the wall in a direction normal to the wall.

Now, in the case of incompressible solvers in openFOAM, the whole system is normalised with respect to the fluid density "rho".... so in this case, the dynamic viscosity (mu) is given by:

mu = nu * rho

Where "nu" is the kinematic viscosity, defined in the "transportProperties" dictionary, and "rho" would depend on the fluid you are working with.

Now, to implement this in openFOAM, there exists a class member called "snGrad". This directly gives you the gradient of velocity (dU/dy) normal to each face of the wall patches.

So, you would use:

wallShear = nu * rho * mesh.boundaryField()[patchID].snGrad()

This directly gives you the wall shear stress (assuming the "patchID" refers to a wall patch ofcourse). Here, you can read "nu" from the transportProperties dictionary, and you will have to define "rho" somewhere.

When the case is not laminar, but turbulent, then things look a bit different, which is why you may have gotten confused when you looked through the wallShearStress source-code.

In this case, the Reynolds Stress Tensor gives you the shear stress, but not necessarily normal to the face, and you need to resolve the tensor along the normal to the face in order to obtain the wall shear stress.

This is because in the case of turbulent simulations, the viscosity is not a constant, and is an effective value which is a result of the turbulence model, and the Reynolds Stress Tensor.

So, to calculate the wall shear stress in the case of a turbulent simulation, you need to do something like:

wallShear = (-mesh.Sf().boundaryField()[patchID]/(mag(mesh.Sf().boundaryField()[patchID])) ) & (turbulence->R()().boundaryField()[patchID])

The above equation will give you the wall shear stress with the velocity gradient resolved in the direction normal to each face in the patch.

Have a nice day!

Philippose

[hsieh]

Hi, Philipose,

Thanks for the detailed post.

The solver I used was standard interFoam. So, this is a laminar case.

The question I am having is:

Does U.boundaryField()[patchi].snGrad() contains normal component of shear (shear stress normal to wall)?

In my case, I have a curve wall with strong centrifugal force. I am only interested in the shear stress tangent to the wall.

Do a n dot U.boundaryField()[patchi].snGrad() should remove the normal component, correct?

What do I have to do to get n in post-processing?

Pei

[philippose]

Hello Pei,
A Good day to you!

Looks like you have a small confusion regarding the issue of shear stress at a wall.

Lets go from the basics.... we know that for a Newtonian fluid, shear stress (Tau) is given by the equation:

Tau = mu * (dU/dy)

This can be seen in the image below:



In your case, the "boundary plate (stationary)" would be the wall (i.e., the fluid is moving relative to the wall).

In the image, the term "gradient, dU/dy" is what is given by the openFOAM term:

U.boundaryField()[patchi].snGrad()

It is the "change in velocity U, with distance y, normal to the wall", which is basically, the gradient "dU/dy".

When you multiply this term with the dynamic viscosity of the medium, you get the shear stress on the wall due to the motion of the fluid, and ofcourse, this shear stress is by definition, tangential to the wall.

So you dont need to take a dot product of the snGrad with a vector tangential to the surface, because by definition, tangential wall shear stress is proportional to the gradient of velocity normal to the wall.

[santos]

Hello Philippose

I have a question regarding your good explanation of wall shear stress.

When we calculate:

U.boundaryField()[patchi].snGrad()

we get the velocity gradient at the wall, and in OpenFOAM definition, U is a vector with components x, y and z, so the previous calculation would give us (assuming y is normal to the wall):

(dUx/dy, dUy/dy, dUz/dy)

Isnt dUy/dy a normal stress? And the shear stress wouldnt be calculated using the other two components?

Regards,
José Santos

[philippose]

Hi José,
Interesting question! In order to answer that, I guess we will have to look deeper into the implementation of the "snGrad" function in OpenFOAM.

I was trying to locate the definition of the function, and found something in the "fvPatchField.H" file in the "finiteVolume/lnInclude" directory of the openFOAM source.

As far as I understood, the snGrad function takes the difference between the value (in this case U) at the patch face, and in the internal field next to the patch, and divides it by the distance from the patch face to the centre of the corresponding cell.

So technically, if there exists a velocity component in a direction normal to the wall, that would come up as a term in the stress calculation (I hope someone corrects me if I am wrong :-)!)

On the other hand, when one talks about laminar flow, one of the basic assumptions is that the fluid flows in the form of shear layers in regions close to walls, and at a distance of one cell height from the wall, you will not (must not / should not) have flow in a direction normal to the wall... if a normal component exists so close to a wall, it would be more or less violating the laminar flow assumptions.

If I am wrong, it would be great if someone could correct me :-)!

Have a nice day!

Philippose

[hsieh]

Thanks Philippose!

I get it now.

Pei

[giovanni]

Hi all,

I'm trying to use Markus's

volTensorField stressTensor = - mu * (fvc::grad(U) + (fvc::grad(U)).T());

in a rhoSimpleFoam case for OpenFOAM 1.4, but appears this error message:

error: 'class Foam::tmp<foam::geometricfield<foam::tensor<double >, Foam::fvPatchField, Foam::volMesh> >' has no member named 'T'

Can anyone help me to solve this problem?

Thanks a lot!


Giovanni

[hartinger]

Hi Giovanni,

the stuff i posted wasn't correct, do it like that:

volTensorField gradU = fvc::grad(U);

// stress tensor
volTensorField tau =
- mu * (gradU + gradU.T());

[giovanni]

Thank you Markus,

I'll try as soon as possible (this afternoon I think)!

Have a nice day!

Giovanni

[ngj]

Hi and good evening

I have been reading the above discussion several times, and I have been thinking what to do about the bed shear stresses. I like the procedure of simply calculating

dU_dn = U.boundaryField().snGrad()

though in my case, I will have vertical velocity components in the cell next to the wall, as I consider non-stationary flows such as wave boundary layers. Thus I will suggest a slightly modified version of Phillipposes (please consider if I am off track).

As I get a non-zero surface-normal component in dU_dn the bed shear stress is not aligned with the bed, which is wrong. Thus I suggest to make a projection of dU_dn onto the boundary patch. Then the bed shear stress will definitely be aligned with the bed, and therefor being shear stress components only and not a combination of shear and normal stresses.

Based on this: http://www.euclideanspace.com/maths/geometry/elements/plane/lineOnPlane/index.ht m l
the projection reads:

tau / mu = sf x (dU_dn x sf / magSf) / magSf

where sf is the surface normal gradient, magSf is the magnitude of sf and x takes the role of the cross product. I have projected dU_dn, since snGrad() takes care of the non-orthogonality of the mesh.

I would like to hear if you agree and especially if there are any numerical traps which I have walked straight into.

Best regards,

Niels

[sven]

Hello,

in this thread you it is often mentioned to use something like "U.boundaryField()[patchi].snGrad()" to calculate a gradient. How do I do this? is this a command for bash or what do I have to do to use this function for calculating a gradient??