[cfdopenfoam]

Hi,

Does the fixedshearstress B.C. set the patch normal component of velocity Un to 0? Could someone guide me the math of this B.C. and explain the formula in the source code?

regards,

[mAlletto]

From the description in the header file: Set a constant shear stress as tau0 = -nuEff dU/dn.


it is derived from th3 fixedValueFvPatchVectorField. That means it sets the value at the boundary.


rearrangin tau0 = -nuEff dU/dn gives:



operator==(tauHat*(tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc)));


vector tauHat = tau0_/(mag(tau0_) + ROOTVSMALL); -> unit vector in direction of the shear


tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc)) -> consider only shearstress (i.e. the shear in direction of tau0_ and exclude normal stress)



(tauHat*(tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc))); -> update velocity components only in direction of tau0_


i hope the explanation is understandeble.

[cfdopenfoam]

Originally Posted by mAlletto

From the description in the header file: Set a constant shear stress as tau0 = -nuEff dU/dn.


it is derived from th3 fixedValueFvPatchVectorField. That means it sets the value at the boundary.


rearrangin tau0 = -nuEff dU/dn gives:



operator==(tauHat*(tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc)));


vector tauHat = tau0_/(mag(tau0_) + ROOTVSMALL); -> unit vector in direction of the shear


tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc)) -> consider only shearstress (i.e. the shear in direction of tau0_ and exclude normal stress)



(tauHat*(tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc))); -> update velocity components only in direction of tau0_


i hope the explanation is understandeble.

Thank you for your explanation, Michael. A few more questions:

1) how does the negative "-" comes from in tau0 = -nuEff dU/dn? As I see in textbooks, its like tau0 = nuEff dU/dy.



2) look back to my original question, does this B.C. set Un to 0? I see a lot of code doing mathmatical operations like this:

const vectorField nHat(this->patch().nf()); // get patch normal vector (unit) and do the following
nHat*(nHat & U)

could you please explain what we can get from the above formula and what is the result? Is it giving the patch-normal component of U? and how? I cannot understand this formula mathmatically.




3) Also, we see the following at a lot of place:

transform(I - sqr(nHat), this->patchInternalField()) // see all mentioned code in fixedNormalSlipFvPatchField.C

Does this give us the component of patchInternalField() parallel to the patch?

So many questions. Thank you for your patience.

[mAlletto]

1) how does the negative "-" comes from in tau0 = -nuEff dU/dn? As I see in textbooks, its like tau0 = nuEff dU/dy.


The minus sign comes from the fact the the boundary conditions applied model the shear force exerted by the boundary on the fluid. If someone wants to calculate the force exerted by the fluid on the boundary the plus sign is the correct sign.


2) look back to my original question, does this B.C. set Un to 0? I see a lot of code doing mathmatical operations like this:

It does not set the velocity Un to zero. This boundary condition only imposes the velocty components parallel to the patch and extrapolates the components normal to the patch from the cell centre next to the boundary.


const vectorField nHat(this->patch().nf()); // get patch normal vector (unit) and do the following
nHat*(nHat & U)

could you please explain what we can get from the above formula and what is the result? Is it giving the patch-normal component of U? and how? I cannot understand this formula mathmatically.

this-> patch().nf() returns the list of face normal vectors at the boundary

nHat*(nHat & U) is the flow velocity component in the direction of the face normal.


nHat & U is the dot product between the velocity at the boundary and the face normals (it gives the projection of U in patch normal direction).

by multiplying the projection of U in patch normal direction with the patch normal vector you get the component of U in direction of the patch normal.



3) Also, we see the following at a lot of place:

transform(I - sqr(nHat), this->patchInternalField()) // see all mentioned code in fixedNormalSlipFvPatchField.C

Does this give us the component of patchInternalField() parallel to the patch?

So many questions. Thank you for your patience.


I cannot find the above code in the newer OF version (4.x or 5.x) which one are you using.

[cfdopenfoam]

Originally Posted by mAlletto

1) how does the negative "-" comes from in tau0 = -nuEff dU/dn? As I see in textbooks, its like tau0 = nuEff dU/dy.


The minus sign comes from the fact the the boundary conditions applied model the shear force exerted by the boundary on the fluid. If someone wants to calculate the force exerted by the fluid on the boundary the plus sign is the correct sign.


2) look back to my original question, does this B.C. set Un to 0? I see a lot of code doing mathmatical operations like this:

It does not set the velocity Un to zero. This boundary condition only imposes the velocty components parallel to the patch and extrapolates the components normal to the patch from the cell centre next to the boundary.


const vectorField nHat(this->patch().nf()); // get patch normal vector (unit) and do the following
nHat*(nHat & U)

could you please explain what we can get from the above formula and what is the result? Is it giving the patch-normal component of U? and how? I cannot understand this formula mathmatically.

this-> patch().nf() returns the list of face normal vectors at the boundary

nHat*(nHat & U) is the flow velocity component in the direction of the face normal.


nHat & U is the dot product between the velocity at the boundary and the face normals (it gives the projection of U in patch normal direction).

by multiplying the projection of U in patch normal direction with the patch normal vector you get the component of U in direction of the patch normal.



3) Also, we see the following at a lot of place:

transform(I - sqr(nHat), this->patchInternalField()) // see all mentioned code in fixedNormalSlipFvPatchField.C

Does this give us the component of patchInternalField() parallel to the patch?

So many questions. Thank you for your patience.


I cannot find the above code in the newer OF version (4.x or 5.x) which one are you using.

Thank you very much, Michael. I may have different views into these questions. Pls correct me if I am wrong anywhere.


For question 1), as said, tau0 = -nuEff dU/dn. I think the minus sign comes from the fact that n (vector) is the out-pointing normal vector from the patch face, and so is dn.

In many textbooks, we see tau0 = nuEff dU/dy because y is pointing into the fluid region. In this context, tau0 = nuEff dU/dy and dU = U_y - U_y0 and dy = y - y0.

When using tau0 = -nuEff dU/dn, there should be dU = U_cell - U_patch and dn = point_face - point_cell such that the minus sign is needed.



For question 2), I guess the fixedShearStress b.c. forces Un = 0.

Firstly, this b.c. is in group of grpWallBoundaryConditions (see fixedShearStressFvPatchVectorField.H). From the word "wall", I guess this b.c. is for patch of a wall type.

Secondly, following question 1), (tau0_*(1.0/(ry*nuEff)) + Uc) gives U_patch. Then tauHat*(tauHat & U_patch) gives U_patch component in the direction of tauHat (or parallel to the patch face if tau0_ is parallel to the patch face).

At last operator==(tauHat*(tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc))) only sets the velocity at patch to be in the direction of tau0_ and Un will be zero if tau0_ is parallel to the patch face.



For question 3), I am using OF-2.3.0. I've checked. You can see this piece of code in fixedNormalSlipFvPatchField.C (void Foam::fixedNormalSlipFvPatchField<Type>::evaluate( )), even in 5.x.

[mAlletto]

For question 1), as said, tau0 = -nuEff dU/dn. I think the minus sign comes from the fact that n (vector) is the out-pointing normal vector from the patch face, and so is dn.

In many textbooks, we see tau0 = nuEff dU/dy because y is pointing into the fluid region. In this context, tau0 = nuEff dU/dy and dU = U_y - U_y0 and dy = y - y0.

When using tau0 = -nuEff dU/dn, there should be dU = U_cell - U_patch and dn = point_face - point_cell such that the minus sign is needed.




Of course the sign of tau0 = -nuEff dU/dn. depends on where the normal vector is pointing. In classical textbook the expression tau0 = nuEff dU / dy is the force exerted by the fluid on the wall. If you think of a Couette flow which is flowing in x direction and the velocity is increasing in y direction (the bottom wall) tau0 points in positive x direction, hence it is the force exerted by the fluid on the wall.





For question 2), I guess the fixedShearStress b.c. forces Un = 0.

Firstly, this b.c. is in group of grpWallBoundaryConditions (see fixedShearStressFvPatchVectorField.H). From the word "wall", I guess this b.c. is for patch of a wall type.

Secondly, following question 1), (tau0_*(1.0/(ry*nuEff)) + Uc) gives U_patch. Then tauHat*(tauHat & U_patch) gives U_patch component in the direction of tauHat (or parallel to the patch face if tau0_ is parallel to the patch face).

At last operator==(tauHat*(tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc))) only sets the velocity at patch to be in the direction of tau0_ and Un will be zero if tau0_ is parallel to the patch face.

You're right: i oversee a bracket. the dot product tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc) leads to a zero component in the direction normal to tauHat. Hence the velocity in this direction is set to zero.

For question 3), I am using OF-2.3.0. I've checked. You can see this piece of code in fixedNormalSlipFvPatchField.C (void Foam::fixedNormalSlipFvPatchField<Type>::evaluate( )), even in 5.x.[/QUOTE]


I thought we where talking about fixedShearStressFvPatchVectorField.C

[cfdopenfoam]

Originally Posted by mAlletto

For question 1), as said, tau0 = -nuEff dU/dn. I think the minus sign comes from the fact that n (vector) is the out-pointing normal vector from the patch face, and so is dn.

In many textbooks, we see tau0 = nuEff dU/dy because y is pointing into the fluid region. In this context, tau0 = nuEff dU/dy and dU = U_y - U_y0 and dy = y - y0.

When using tau0 = -nuEff dU/dn, there should be dU = U_cell - U_patch and dn = point_face - point_cell such that the minus sign is needed.




Of course the sign of tau0 = -nuEff dU/dn. depends on where the normal vector is pointing. In classical textbook the expression tau0 = nuEff dU / dy is the force exerted by the fluid on the wall. If you think of a Couette flow which is flowing in x direction and the velocity is increasing in y direction (the bottom wall) tau0 points in positive x direction, hence it is the force exerted by the fluid on the wall.





For question 2), I guess the fixedShearStress b.c. forces Un = 0.

Firstly, this b.c. is in group of grpWallBoundaryConditions (see fixedShearStressFvPatchVectorField.H). From the word "wall", I guess this b.c. is for patch of a wall type.

Secondly, following question 1), (tau0_*(1.0/(ry*nuEff)) + Uc) gives U_patch. Then tauHat*(tauHat & U_patch) gives U_patch component in the direction of tauHat (or parallel to the patch face if tau0_ is parallel to the patch face).

At last operator==(tauHat*(tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc))) only sets the velocity at patch to be in the direction of tau0_ and Un will be zero if tau0_ is parallel to the patch face.

You're right: i oversee a bracket. the dot product tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc) leads to a zero component in the direction normal to tauHat. Hence the velocity in this direction is set to zero.

For question 3), I am using OF-2.3.0. I've checked. You can see this piece of code in fixedNormalSlipFvPatchField.C (void Foam::fixedNormalSlipFvPatchField<Type>::evaluate( )), even in 5.x.

I thought we where talking about fixedShearStressFvPatchVectorField.C[/QUOTE]

Thank you for your help, Michael. Now clear. Hope discussion here would help someone else.

[cfdopenfoam]

Originally Posted by mAlletto

For question 1), as said, tau0 = -nuEff dU/dn. I think the minus sign comes from the fact that n (vector) is the out-pointing normal vector from the patch face, and so is dn.

In many textbooks, we see tau0 = nuEff dU/dy because y is pointing into the fluid region. In this context, tau0 = nuEff dU/dy and dU = U_y - U_y0 and dy = y - y0.

When using tau0 = -nuEff dU/dn, there should be dU = U_cell - U_patch and dn = point_face - point_cell such that the minus sign is needed.




Of course the sign of tau0 = -nuEff dU/dn. depends on where the normal vector is pointing. In classical textbook the expression tau0 = nuEff dU / dy is the force exerted by the fluid on the wall. If you think of a Couette flow which is flowing in x direction and the velocity is increasing in y direction (the bottom wall) tau0 points in positive x direction, hence it is the force exerted by the fluid on the wall.





For question 2), I guess the fixedShearStress b.c. forces Un = 0.

Firstly, this b.c. is in group of grpWallBoundaryConditions (see fixedShearStressFvPatchVectorField.H). From the word "wall", I guess this b.c. is for patch of a wall type.

Secondly, following question 1), (tau0_*(1.0/(ry*nuEff)) + Uc) gives U_patch. Then tauHat*(tauHat & U_patch) gives U_patch component in the direction of tauHat (or parallel to the patch face if tau0_ is parallel to the patch face).

At last operator==(tauHat*(tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc))) only sets the velocity at patch to be in the direction of tau0_ and Un will be zero if tau0_ is parallel to the patch face.

You're right: i oversee a bracket. the dot product tauHat & (tau0_*(1.0/(ry*nuEff)) + Uc) leads to a zero component in the direction normal to tauHat. Hence the velocity in this direction is set to zero.

For question 3), I am using OF-2.3.0. I've checked. You can see this piece of code in fixedNormalSlipFvPatchField.C (void Foam::fixedNormalSlipFvPatchField<Type>::evaluate( )), even in 5.x.

I thought we where talking about fixedShearStressFvPatchVectorField.C[/QUOTE]


Hi，I got one more question for the fixedShearStress boundary condition. Hope you could help. Its about nuEff in tau0 = -nuEff dU/dn. I think we should use nuEff at said patch. However from the source it does not look like in that way. See line 112.

scalarField nuEff(turbModel.nuEff(patch().index())); // what does patch().index() give us? Does this mean we use nuEff at patch connected cell?

Should not it be something like:

const fvPatchField<scalar>& nuEff =
patch().lookupPatchField<volScalarField, scalar>("nuEff");

Any hints?

[mAlletto]

patch().index() gives back the index of this patch in the boundary mesh.



https://www.openfoam.com/documentati...732edc13c119c7

[mAlletto]

the call you described probably works too, since it it used also in https://www.openfoam.com/documentati...ce.html#l00104


so probably bath ways to assess the values of a given field at a given patch is works. I actually never tried it out. The best would be to write both commands in an on boundary condition and try it out if the list of values they give back are actually the same

[cfdopenfoam]

Originally Posted by mAlletto

patch().index() gives back the index of this patch in the boundary mesh.



https://www.openfoam.com/documentati...732edc13c119c7

Thank you so much for your discussion here. One more question for fixedShearStress B.C. in fixedShearStressFvPatchVectorField.C.

In Eqn tau0 = -nuEff dU/dn, from the source code (of-2.3.0) I see it uses the Kinematic viscosity. Should not it be the dynamic viscosity?

Correct me if I am wrong anywhere.

[mAlletto]

this is because for the incompressible solvers the equations solved are divided by the density. Hence you get the kinematic viscosity.

[cfdopenfoam]

Originally Posted by mAlletto

this is because for the incompressible solvers the equations solved are divided by the density. Hence you get the kinematic viscosity.

Do you mean the shear stress, tau0_ to be specified, has been divided by the density? If it is, we should be careful to specify the value of tau0_ because physically its tau/rho.

[mAlletto]

Exactly. Since for the incompressible solvers the momentum equation is divided by the density, it is tau_0/rho which has to be specified.