[Yoshiki]

Hello,

I am using OpenFOAM v2112 to calculate wall shear stress (WSS) in a straight circular pipe model and compare it with theoretical values. For the analysis, I am using the wallshearstress utility under postProcess.

When I compare the computational WSS values with the theoretical ones, the computed value is approximately 0.3 Pa, which is about 50% lower than the theoretical value of approximately 0.6 Pa. I am struggling to resolve this issue.

I have tried varying the mesh size (0.1 mm, 0.3 mm, 0.5 mm) and boundary layer settings, but the computational result remains around 0.3 Pa in all cases.

Model specifications:
Straight circular pipe:
Diameter: 0.005 m
Length: 0.05 m
Mesh:
Generated using cfMesh [cartesianMesh (hexahedral mesh)]
Mesh size: 0.3 mm
Boundary layer: 10 layers, each 0.01 mm thick
Boundary conditions:
Inlet: Average velocity of 0.21 m/s applied
Outlet: patch (pressure condition)
Walls: No-slip condition
Other settings:
Solver: simpleFoam (steady-state)
Laminar flow
Newtonian fluid with a kinematic viscosity of
𝜈=3.301×10−6m^2/s

For post-processing, I used ParaView to calculate WSS by multiplying the output of the wallshearstress utility with the fluid density.

I also tried running the same analysis with pisoFoam (unsteady simulation), but the WSS still remained around 0.3 Pa, which is lower than expected.

If anyone could provide suggestions on how to accurately calculate WSS, I would greatly appreciate it.

velocity.xlsx

[Severus]

Hello,
The length of your pipe isn't long enough for the flow to fully develop.
Since you compare it with the theoretical estimate, I assume the theoretical value is derived for a fully developed pipe flow.

Based on your inlet conditions, your Re is 318 and hence your development length (or entrance length) is 0.09m.
So, make your pipe longer, lets say > 0.1m and measure the wall shear stress at a distance after 0.09m.

To calculate the development length use:
https://www.engineeringtoolbox.com/e...low-d_615.html
https://en.wikipedia.org/wiki/Entran...ntrance_length


Good luck!

[Yoshiki]

Thank you, Severus.

At the pipe inlet, a parabolic velocity profile with a maximum velocity of 0.42 m/s and an average velocity of 0.21 m/s is specified.　U.txt
The pipe radius is 0.005 m. I apologize for the missing information and corrections in my earlier message.

This time, since the flow at the inlet is fully developed, the pipe length should not be an issue.

Additionally, when a uniform flow of 0.21 m/s is applied at the inlet to a pipe with a radius of 0.005 m and a length of 0.45 m, and the wall shear stress is calculated at a distance of 0.4 m downstream, the wall shear stress is about 50% smaller than the theoretical value.

I believe there might be another reason why the wall shear stress does not match the theoretical solution. If you have any ideas, please let me know.

[Severus]

Hello,
I understand, thanks for the clarification.

I have some suggestions:
1) If your flow is already fully developed, simply compute the wall shear stress from the pressure drop between the inlet and outlet. See for example here: https://resources.system-analysis.ca...through-a-pipe
This value should be the same as what you obtain from the functionObject. Also double check if the theoretical relations you are using are correct.

2) Make sure you are not missing a factor 2 when performing the theoretical calculation (also see the above link) or setting up your geometry for the simulation.
It is easy to mess it up when using one instead of the other, for e.g. radius, diameter, hydraulic diameter etc.

3) Check your solver residuals and make sure that your simulation has converged well. If not simulate it for a longer time until convergence is achieved.

4) Also check your mesh and make sure everything is fine with the checkMesh command.

Thanks

[Yoshiki]

Thank you for your suggestions.

1) When comparing the pressure loss with the theoretical values, the pressure loss closely matches the theoretical values.
2) I have verified the coefficients multiple times, and they are correct.
3) The number of iterations has been set to 1000, and based on the residuals, it seems to have converged.
4) Regarding the mesh, I ran checkMesh, and it passed with "mesh OK," indicating no abnormal mesh issues.

Since I am using the wallshearstress utility directly from OpenFOAM, I suspect that there may be an error in the wallshearstress utility itself (perhaps in the gradient calculation or other aspects).

Could you suggest a more accurate way to compute wall shear stress, or let me know if there is anything that needs to be added to the code of the wallshearstress utility?

Thank you.

[Tobermory]

If the pressure drop down the pipe matches the theoretical value, then so should the wall shear stress! You can do an integral momentum balance:




as a rough check on your theoretical values ... maybe there is a mistake there?

[Yoshiki]

The wall shear stress value calculated using the pressure drop equation is almost consistent with the theoretical value, so I believe there is no issue with the calculations for pressure and other wall shear stress-related calculations.

However, the value of wall shear stress calculated using the wallshearstress utility and then computed in ParaView's Calculator is about 50% smaller than the theoretical value.

I believe that changes to the source code are necessary when using the wallshearstress utility, but I'm not sure how to proceed.

If you have any other suggestions or potential solutions, I would appreciate it.

Thank you.

[Severus]

Hello,
I have two suggestions,
1) Now you use Paraview to simply multiply the wallShearStress output with the density. What you can additionally do is, you can also compute the wall shear stress directly in Paraview using the velocity field.
Lets say your flow is in the x-direction. Use the Probe filter, and measure the x-velocity value at the y location of the cell centre of the first cell from the wall (lets say at height h from wall). Then simply compute wall shear stress: rho*nu*du/dy = rho*nu*(u_probe_x - 0)/h
If you think the functionObject of OpenFOAM has a bug, compute it this way and compare it to the one you obtained using the functionObject.

2) I don’t expect any bug in the code. But you can check that as well. There are multiple tutorials that are used as validation cases, where you can employ the wallShearStress functionObject.
You can run one of them and verify if there is a bug.
See here: https://www.openfoam.com/documentati...alidation.html
Take a tutorial where Cf is compared and validated. Run the tutorial yourself and compute Cf using wallShearStress and see if it matches.

Let us know the results

[Yoshiki]

Thank you, Severus.

Apologies for the delayed response. I’d like to share the results of my verification:

Using the Paraview Probe filter, I calculated the wall shear stress from the velocity at the center of the closest cell to the wall. The results are as follows:

At the inlet cross-section, the value was 0.9 [Pa].
At the center cross-section (z = 0.025), it was approximately 0.6 [Pa] (which is very close to the theoretical value).
At the outlet cross-section, the value was 0.83 [Pa].
The values vary depending on the location. In contrast, the wall shear stress (WSS) calculated using the functionObject "wallshearstress" utility was about 0.3 [Pa], which is significantly different from the values I calculated above.

For the Turbulent Flow bump (2D) tutorial, I performed a validation of the Skin Friction Coefficient (Cf) using the following link:
https://www.openfoam.com/documentati...t-bump-2d.html

I used the WSS calculated from the "wallshearstress" utility and the inlet velocity of 69.44 [m/s] to compute Cf.
https://www.cfd-online.com/Wiki/Skin...on_coefficient
Comparing my calculated Cf with the results from the website, the shape of the distribution is the same, but the minimum value was around 0.002, which is about twice the minimum value of 0.001 on the website. The mesh and model size were the same as those used in the tutorial.

Cf_validation.xlsx

If you have any insights or ideas, I would greatly appreciate your help.

[Severus]

If your flow is fully developed, the value of the wallShearStress should not change with x-direction (along the pipe length after flow development)
Note that for a fully developed pipe flow, du/dx = 0 (i.e. the flow is homogenous in x).
Since, you mentioned that the value of the wallShearStress computed from Paraview changes so much, I doubt if the flow has converged and reached a fully developed state.
To check if the flow is fully developed, the inlet velocity profile that you feed in should not change with x. You can plot u (i.e. U_x) at the center of the pipe from x = 0 till x = L (inlet to outlet). It should remain a constant.

Also one question, does the value of wallShearStress computed using the OpenFOAM functionObject change along x?

To answer your 2nd result, where you computed Cf.
I suspect that some changes with the case setting or the choice of the turbulence model is affecting your result.
I simply ran the same tutorial with OpenFOAM2206 and the result looks fine (see the plot attached).
Did you setup the case from scratch? I simply ran the ./Allrun command and I obtained the nice plots.
Simply try it out first. Copy the case from the tutorials/ folder and then run it using the ./Allrun command.

[Yoshiki]

1) In the straight pipe, the wall shear stress varies depending on the location. （The theoretical value is 0.6 Pa）
wss.png

Additionally, the velocity distribution in the flow direction (z-direction) at the y-coordinate of the cell center closest to the wall was not constant. (Y=0.00498462)
U.png wss_magnitude.png

I believe that the variation in the wall shear stress distribution across locations is a problem.

Model specifications:

Straight circular pipe:
Radius: 0.005 m
Length: 0.05 m
Mesh:
cfMesh [Cartesian mesh (hexahedral mesh)]
Mesh size: 0.3 mm
Boundary layer: 5 layers, each with thickness of 0.03 mm
Boundary conditions:
Inlet: Parabolic velocity profile with an average velocity of 0.21 m/s
Outlet: Patch (pressure condition)
Wall: No-slip condition
Other settings:
Solver: simpleFoam (steady-state)
Kinematic viscosity: ν = 3.301 × 10⁻⁶ m²/s

2) I changed from 8 cores to 6 cores for the OpenFOAM 2112 tutorial on bump 2D, and ran ./Allrun without changing other parts, but an error occurred. I then analyzed the turbulence models one by one.
The results for k-epsilon almost matched your Cf distribution.
The results for the Spalart model also nearly matched.
However, for k-omega SST, the solution converged quickly, but the distribution shape was different.
For the k-epsilon and Spalart models, I think the Cf values were successfully derived from the wall shear stress.

Could you please suggest any improvements?
Thank you.

[Severus]

Hello,
Now it is clear that the flow has either not yet converged or developed and that is why the velocity along the pipe length doesn’t stay constant.

I think there could be some problems with the case setup or the mesh.

I suggest you to do the following:
1) Check if your case converged
- plot the residuals
- probe the value of velocity at the centre at a specific location. By doing this you can monitor centre velocity data (https://doc.openfoam.com/2306/tools/...mpling/probes/)
2) Perform a grid insensitivity study and refine your mesh. May be you need a finer mesh.
3) Do not use a fixed time step but switch on adjustTimeStep with maxCo = 0.3 or 0.5
4) Increase outerCorrectors to 3 or even higher in fvSolutions
5) Share with us your fvSchemes and fvSolutions file. We can try to see if there are any issues with that.

Other option - trying out a new case
Simulate a long pipe with uniform velocity at inlet and then measure wallShearStress when it is fully developed.
Check out these videos for your reference: https://www.youtube.com/watch?v=6aIK...hoUnBVf5x6egPp
When doing this, you could maybe realize if you missed something in the case setup.

[Yoshiki]

Thank you.
I will try your suggestions.

5) I will post the current code of fvSchemes and fvSolution below.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2112                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;

    div(phi,U)      bounded Gauss linearUpwindV grad(U);

    turbulence      bounded Gauss upwind;
    div(phi,k)      $turbulence;
    div(phi,omega)  $turbulence;

    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method          meshWave;
}


// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2112                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-7;
        relTol          0.01;
    }

    Phi
    {
        $p;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
        nSweeps         1;
    }

    k
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
        nSweeps         1;
    }

    omega
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
        nSweeps         1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent yes;
}

potentialFlow
{
    nNonOrthogonalCorrectors 10;
}

relaxationFactors
{
    equations
    {
        U               0.8;
        k               0.6;
        omega           0.6;
    }
}

cache
{
    grad(U);
}


// ************************************************************************* //

[Severus]

Hello,
I forgot that you were using simpleFoam. You can ignore my 3rd suggestion.
Also, I am not how much the 4th suggestion might help, you can try in any case.

After taking a look at your fvSolutions, the final values are missing (pFinal UFinal etc). If you don't specify, I am not sure what values it will take by default.
Follow here: https://develop.openfoam.com/Develop...?ref_type=tags

Thanks

[Yoshiki]

Originally Posted by Yoshiki

1)
The velocity distribution in the flow direction (z-direction) at the y-coordinate of the cell center closest to the wall was not constant. (Y=0.00498462)

The velocity distribution in the flow direction (z-direction) at the y-coordinate of the cell center closest to the wall showed that, except for regions near the inlet and outlet, the velocity along the z-axis had an error of less than 1%.

Additionally, the velocity at the center (Y=0) also showed a small error of around 0.05%.

Also, when pFinal and UFinal were set to 1e-2 and 1e-3 respectively, the analysis converged in 114 iterations, and the results were almost the same as those obtained from a 1000-iteration analysis without setting pFinal and UFinal.

[Severus]

The tolerance for the UFinal and pFinal should not be that high. Note that these are tolerances of the linear solver and not the one you see when you monitor the residuals.
Since you set them so high, it converges faster. Setting such high values for the tolerance is not generally recommended.
Check the tutorial in the link from my previous post and follow that.

[Yoshiki]

Thank you for your feedback.

I have posted the corrected fvSolution here.

solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-7;
        relTol          0.01;
        outerCorrectors 3;

    }

    pFinal
    {
        $p;
        tolerance       1e-06;
        relTol          0;
    }


    Phi
    {
        $p;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
        nSweeps         1;
        outerCorrectors 3;
    }

    UFinal
    {
        $U;
        tolerance       1e-05;
        relTol          0;
    }

    k
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
        nSweeps         1;
    }

    omega
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
        nSweeps         1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent yes;
}

potentialFlow
{
    nNonOrthogonalCorrectors 10;
}

relaxationFactors
{
    equations
    {
        U               0.8;
        k               0.6;
        omega           0.6;
    }
}

cache
{
    grad(U);
}

Also, when I used the probe function to check the velocity at the center of the outlet cross-section,
the center of the center plane in the length direction,
and the velocity at the cell center closest to the wall, all of them stabilized at a constant value by iteration 120.
Also, the results for wall shear stress (wallshearstress) did not change much, even after modifying the fvSolution.

[Severus]

Your fvSolutions looks fine.
As you use simpleFoam, outerCorrectors shouldn't really matter.
One issue that I see is where you have placed the outerCorrectors.
It has nothing to do with your current issue with this problem, but letting you know so that you don't make this mistake in the future when simulating transient problems using pimpleFoam.
The outerCorrectors should be inside the PIMPLE section. Not in the solver section.
And the syntax is It is nOuterCorrectors not outerCorrectors. See here: https://develop.openfoam.com/Develop...?ref_type=tags


Now coming back to the original problem, you mentioned that you probed the centreline velocity and the velocity of the first cell centre from the wall. How does these values (after stabilizing) compare to the values used at the inlet (using the parabolic velocity profile relation)? They should be the same. If they are not, then I would more strongly suspect the mesh as the main reason for this issue. You should perform a grid insensitivity and make your mesh finer.

One more thing, your wallShearStress at inlet should be the same as theoretical value (because it is the boundary condition). If it doesn't match, then there is some issue with your inlet parabolic velocity profile boundary condition.

[Alczem]

Hey!


Quick digression, but UFinal, pFinal etc should not be read in a steady-state run, right? It should only affect the final outer iteration of the pimple loop in a transient simulation if I am not mistaken. How could it change the number of iterations required for your case? To support this, if OpenFoam needed these settings, you should not have been able to run the simulation when you were missing them.

[Severus]

Thanks a lot for pointing it out.
Since I have always used either PISO or PIMPLE I was not aware of this, good to know.
I did some reading and indeed you dont need the UFinal and pFinal values when running SIMPLE.

Just to check if OpenFOAM chooses Final values over the regular U/p values, I took a tutorial: https://develop.openfoam.com/Develop...?ref_type=tags
And ran two cases - one case with and one case without UFinal/pFinal. The log file was exactly the same.
The solver simple ignores the Final values and I get the same final residual values.
So I can confirm that the solver simply neglects it.

I am also puzzled as to how and why it converged faster for Yoshiki.