Measuring wall shear stress in bend pipe
My geometry is a bifurcation and on of the vessel is bent. I want to measure the wall shear stress along the bent vessel. The wall shear stress is already calculated. I only find a way to measure the wall shear stress along the straight pipe but no idea how to do it in bent pipe.
Any hind would be helpful. Cheers! |
You can use the utility wallShearStress.
|
Quote:
I think wallShearStress utility is for turbulent flow. My case is laminar flow and I want to plot the wall shear stress along the pipe. My code is to find laminar flow but don't know how to measure it. Thanks for the hind anyway. |
Of course you can still use that utility. I suppose it reads turbulenceProperties or RASProperties. You set it to laminar there. Otherwise you make your own version of the wallShearStress utility and update it for laminar flow. I think the first method will work, since the RASModels have a dummy laminar model included.
|
Quote:
I think this utility needs some modifications because it computes wall traction instead of wall shear stress, right? |
Quote:
What modification have you done to measure the wall shear stress ? Thanks |
Quote:
I will try it later today. However, i still don't understand how to plot wall shear stress along a bent pipe. Have you got any idea for that? Thanks |
Quote:
In many solvers such as FLUENT and OpenFOAM wall traction is considered as wall shear stress and I don't know why! If we assume T as stress tensor and n as unit normal vector of desired face we have: t= T.n ; t= traction (force vector exerted to desired face per unit area); the utility compute this, but we want: shear stress= t - (t.n)n ;which shear stress is tangential vector here. Bests, |
Can any one tell me the difference between wallShearStress and the definition?
For newtonian flow the wall shear stress is defined as mu*du/dy which is proportional to the normal velocity gradient to the wall But in wallShearStress it is defined: wallShearStress.boundaryField()[patchi] = ( -mesh.Sf().boundaryField()[patchi] /mesh.magSf().boundaryField()[patchi] ) & Reff.boundaryField()[patchi]; -- What does this code mean? why it calculated wall shear stress ? In wallGradU: The velocity gradient is defined as wallGradU.boundaryField()[patchi] = -U.boundaryField()[patchi].snGrad(); -- This makes sense to me. And mu*wallGradU is the definition one. Can anyone explain? |
Quote:
Thanks for your reply. To my knowledge, the wall shear stress is mu*velocity gradient. In the code of wallShearStress: -mesh.Sf().boundaryField()[patchi] /mesh.magSf().boundaryField()[patchi] Do you mean this is the T(stress tensor)? Why the say shear stress = t-(t.n)n ? If T is the shear tensor, isn't T's normal component T.n already the shear stress? I am confused here. |
Quote:
T is stress tensor and n is unit normal vector, according to openFOAM definitions: Code:
n=-mesh.Sf().boundaryField()[patchi] / mesh.magSf().boundaryField()[patchi] Code:
T=Reff.boundaryField()[patchi] Code:
wallShearStress.boundaryField()[patchi] = Code:
shear stress= t - (t.n)n |
Quote:
Thanks for the reply. I found that the tangetient component of traction tensor is t-(t.n)n as you decribed. For the code stuff: Is this the right thing: forAll(real_wall_shear_stress.boundaryField(),patc hi) { wallShearStress.boundaryField()[patchi] = ( -mesh.Sf().boundaryField()[patchi] /mesh.magSf().boundaryField()[patchi] ) & Reff.boundaryField()[patchi];[/CODE] real_wall_shear_stress.boundaryField()[patchi]= wallShearStress-(wallShearStress&(mesh.Sf().boundaryField()[patchi] /mesh.magSf().boundaryField()[patchi]))&(mesh.Sf().boundaryField()[patchi]/mesh.magSf().boundaryField()[patchi]) } the right results since mesh.Sf().boundaryField()[patchi] /mesh.magSf().boundaryField()[patchi] is the "n" Please let me know if there is something wrong. Btw, why the tangential component of t is not t&n but t-(t&n)&n ? Thanks again! |
Quote:
First of all: I'm not sure that inner product is overloaded for 2 variables which one of them is scalar and another is vector, so it's better to use: Code:
t-(t&n)n Code:
t-(t&n)n Code:
volTensorField gradU=fvc::grad(U); Obviously, you can change some part to use other stress tensor objects. Bests, |
Quote:
Thanks for the correction. This code works however the results are a little bit different from it should be in my case. After a little bit research. I found that the definition of Shear tensor is T=2*mu*(gamma_hat)*D where D=1/2(gradU+gradU.T()). gamma_hat is shear rate. For newtonian case , mu is independent of gamma_hat, for non-Newtonian case, mu is a function of gamma_hat. gamma_hat=sqrt(2.0)*mag(symm(gradU)) for Newtonian case. In your code, I found "volTensorField T=mu*(gradU+gradU.T())" and there is a difference in gamma_hat(shear rate) term. Did you miss it or I did something wrong( please correct me if I am wrong). Cheers! |
Quote:
Thanks for the correction. I misunderstood the definition. I try to finer my mesh and see what happened next. Regards |
Dear Amir,
I tried and your code to get some wall shear stress plots. In published paper, they normally define wall shear stress as the way you did. After some comparison, I found that my wall shear stress is quite low than what they did on paper. There are still more than 20% difference between my plots and their plots after I multiple a factor to my plots. Have you got any experience with these kind of problem? I am stuck with it for quite a while. From the other method:The plots are totally different. Do you have any idea with this? Any idea would be appreciated:) Thanks |
thx u very much. it helps a lot.
|
Quote:
First of all, note that this relation for shear stress we discussed about is valid just for incompressible flow; for compressible cases, another term should be added. As you know, in many papers, wall traction is reported instead of wall shear stress, so I suggest you check OpenFOAM utility without changes and see what will happen. Ensure you have reached grid independence solution by examining wall shear stress of boundaries. Another suggestion; there is another utility for evaluating velocity gradient @ boundaries (wallGradU), obviously it can compute grad(U) more precisely which you can also use that in your code. Also try high order schemes for gradient and others. Bests, |
Quote:
Can you tell me how to measure the wall shear stress? I'm confused why my wall shear stress is zero and only a tiny part in inlet is non-zero I am modelling a 3D laminar straight pipe flow by using star ccm |
Can you post your bifurcation geometry as well as how you measure the WSS?
Quote:
|
5 Attachment(s)
Quote:
I'm not sure how to measure WSS I just created a scalar plot and then selected wall in parts and WSS in function Appreciate any useful suggestion Thanks |
I think there is tool under filter is plot over intersection curve.
1) In property inspection field, choose the wall object (whoever the name you named it). 2) Use "Clip" to get rid of the part you do not want to measure. 3) Use "Plot over intersection curve" to plot the wall shear stress over the line which is the intersection of a plane and the wall. Quote:
|
Hi guys,
Is there any way that you can output the wall shear stress at every calculated time step (to use this as a measure of convergence to a steady state)? I would like to incorporate it into the solver but not quite sure what the script would look like. Any help would be greatly appreciated. Cheers, Tian |
Hi,
Try this: solve->monitor->surface .... wall fluxes (wall shear stress) Bests, |
Quote:
// Print out average wall shear stress over topSurface // // find the identification number (e.g. label) for our boundary of interest. label topSurfacePatch = mesh.boundaryMesh().findPatchID("topSurface"); // if we don't have such a patch, warn the user // if (topSurfacePatch==-1) { Info << "Failure to find patch named topSurface for average wall shear stress calc." << endl; } else // calculate the result and do output { // the sum operator implicity loops over the boundary faces and stores the result in avgWSS scalar avgWSS = 0.0; avgWSS = sum(wallShearStress.boundaryField()[topSurfacePatch]); reduce(avgWSS,sumOp<scalar>()); Info << "Monitor: at Time = " << runTime.timeName() << " -- Average wall shear stress over top surface= " << avgWSS <<" [kg/(m*s^2)] " ; } When compiling with wmake, the following message appears: Making dependency list for source file mypisoFoam1.C SOURCE=mypisoFoam1.C ; g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3 -DNoRepository -ftemplate-depth-100 -I/hpcwarwick/openfoam/2.1.0/OpenFOAM-2.1.0/src/turbulenceModels/incompressible/turbulenceModel -I/hpcwarwick/openfoam/2.1.0/OpenFOAM-2.1.0/src/transportModels -I/hpcwarwick/openfoam/2.1.0/OpenFOAM-2.1.0/src/transportModels/incompressible/singlePhaseTransportModel -I/hpcwarwick/openfoam/2.1.0/OpenFOAM-2.1.0/src/finiteVolume/lnInclude -IlnInclude -I. -I/hpcwarwick/openfoam/2.1.0/OpenFOAM-2.1.0/src/OpenFOAM/lnInclude -I/hpcwarwick/openfoam/2.1.0/OpenFOAM-2.1.0/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/mypisoFoam1.o mypisoFoam1.C: In function ‘int main(int, char**)’: mypisoFoam1.C:175: error: ‘wallShearStress’ was not declared in this scope /hpcwarwick/openfoam/2.1.0/OpenFOAM-2.1.0/src/finiteVolume/lnInclude/readPISOControls.H:3: warning: unused variable ‘nOuterCorr’ /hpcwarwick/openfoam/2.1.0/OpenFOAM-2.1.0/src/finiteVolume/lnInclude/readPISOControls.H:15: warning: unused variable ‘transonic’ make: *** [Make/linux64GccDPOpt/mypisoFoam1.o] Error 1 The problem appears to be that it can't find wallShearStress as a variable. How do I make it calculate/read this? Cheers Tian |
1 Attachment(s)
Oh, sorry about previous post, that was a forum conflict! :o
But in openFoam it seems that the original utility compute wall traction instead of wall shear stress and I posted the revised form you can easily embed it in your solver and assess its changing. (note that it's a revised utility now which can be used in this form as a utility or you can put in your desired code with few manipulation) see the attachment ... (the error of your code is that you haven't defined wall shear stress yet) |
Thanks a lot for your reply Amir. Does the old wallShearStress function not actually give wall shear stress then? In your utility it requires a tau file in the 0 directory and therefore gives the resulting error message:
--> FOAM FATAL IO ERROR: cannot find file file: /gpfs/home/eng/mauiie/OpenFOAM/mauiie-2.1.0/run/project/cavity2D/solverTest/0/tau at line 0. From function regIOobject::readStream() in file db/regIOobject/regIOobjectRead.C at line 73. FOAM exiting Is this a file used by the RAS solvers? It appears that this utility is designed for RAS not LES. With the old utility it was possible to make some alterations to the script (mainly changing RAS to LES wherever it occurred) and the utility would just as well with an LES case. Is this possible also with this solver? Also, how would you embed the function into the solver? Would you simply copy and paste it in? Or would you do something more a long the lines I started with? Sorry for all the questions, I am a bit new to OpenFOAM. Please bear with me :) |
1 Attachment(s)
Quote:
Quote:
Quote:
Quote:
|
Thank you again for your help so far.
I've managed to incorporate the wall shear stress into my solver now (after some modifications to createFields.H). Ideally what I want is to output an average of wall shear stress over the top surface of my simulation into a log file (for every deltaT). I have put this script in at the end of my script: // Print out average wall shear stress over topSurface // // find the identification number (e.g. label) for our boundary of interest. label topSurfacePatch = mesh.boundaryMesh().findPatchID("topSurface"); // if we don't have such a patch, warn the user // if (topSurfacePatch==-1) { Info << "Failure to find patch named topSurface for average wall shear stress calc." << endl; } else // calculate the result and do output { // the sum operator implicity loops over the boundary faces and stores the result in avgWSS scalar avgWSS = 0.0; avgWSS = sum(mag(shear.boundaryField()[topSurfacePatch])); reduce(avgWSS,sumOp<scalar>()); Info << "Monitor: at Time = " << runTime.timeName() << " -- Average wall shear stress over top surface= " << avgWSS <<" [kg/(m*s^2)] " ; } It is now outputting the sum of over all of that patch. The only thing I need to know is how to say divide by the total number of points on that face??(this is because I need the average not the sum) Cheers Tian |
Also, the script that gives a line break would be very useful :)
|
Quote:
As you can see; you'll need to compute patch area there (not the numbers!): scalar area = gSum(mesh.magSf().boundaryField()[patchI]); Quote:
Bests, |
Quote:
Why should it be over the patch area rather than the number of points on that patch? Since the shear utility outputs a vector at every point on the patch, surely the more points you have (i.e. a finer resolution of mesh) the greater the sum of these vectors will be. So it is by dividing by the number of points that you would achieve the average. Cheers, Tian |
Quote:
- your formulation: (1/N)*sum(shear) - The correct formulation: (1/area)*sum(shear*area) Bests, |
Ah, I see... I did not realise the calculated values were weighted by the size of the cells!
Looking over the script I have an understanding of all aspects except, why do you write this: volTensorField T=mu*(gradU+gradU.T()); T you say is the stress tensor which is the Jacobian of the velocity (i.e. gradU) and I assume you multiply by mu because it's as good a stage as any to do so. However, why the "+gradU.T()"? Cheers Tian |
Quote:
(T=2*mu*D) for incompressible flow D : symmetric part of grad(U) Bests, |
Hi again!
I seem to be getting far too high values out for wall shear stress. I have implemented a monitor in my solver, so that it outputs average wall shear stress over the top surface of my simulation. The code looks as follows: volTensorField gradU=fvc::grad(U); volTensorField T=mu*(gradU+gradU.T()); volVectorField nn ( IOobject ( "nn", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, vector::zero ); forAll(nn.boundaryField(), patchI) { nn.boundaryField()[patchI] = ( -mesh.Sf().boundaryField()[patchI] /mesh.magSf().boundaryField()[patchI] ); } forAll(wallShearStress.boundaryField(), patchI) { wallShearStress.boundaryField()[patchI] = nn.boundaryField()[patchI] & T.boundaryField()[patchI]; } forAll(normal.boundaryField(), patchI) { normal.boundaryField()[patchI] = (nn.boundaryField()[patchI] & wallShearStress.boundaryField()[patchI]) *nn.boundaryField()[patchI]; } forAll(shear.boundaryField(), patchI) { shear.boundaryField()[patchI] = wallShearStress.boundaryField()[patchI]-normal.boundaryField()[patchI]; } // Print out average wall shear stress over topSurface // find the identification number (e.g. label) for our boundary of interest. label topSurfacePatch = mesh.boundaryMesh().findPatchID("topSurface"); // if we don't have such a patch, warn the user if (topSurfacePatch==-1) { Info << "Failure to find patch named topSurface for average wall shear stress calc." << endl; } else // calculate the result and do output { // the sum operator implicity loops over the boundary faces and stores the result in avgWSS scalar area = gSum(mesh.magSf().boundaryField()[topSurfacePatch]); scalar avgWSS = gSum(mag(shear.boundaryField()[topSurfacePatch]))/area; // reduce(avgWSS,sumOp<scalar>()); Info << "Monitor: at Time = " << runTime.timeName() << " -- Average wall shear stress over top surface= " << avgWSS <<" [kg/(m*s^2)] " << "Area = " << area << " [m^2] "; } The values I am getting out are far too high for my simulation, for example average wall shear stress of 40 where mean velocity is magnitude 1. I have searched through the code trying to work out where the issue is occurring. Do you have any suggestions please? Cheers, Tian |
Hi Tian,
I think (am certain) that the error is in the following line: Code:
scalar avgWSS = gSum(mag(shear.boundaryField()[topSurfacePatch]))/area; Code:
scalar avgWSS = gSum(mag(shear.boundaryField()[topSurfacePatch]) * mesh.magSf().boundaryField()[topSurfacePatch])/area; Niels |
Ahhh... That is so much better.
Thank you Niels :) you're a life saver! Tian |
No problem.
BTW: In stead of making the stress tensor yourself, you could benefit from the turbulence model by calling devReff(). The only thing to remember is that this tensor is also non-tangential with the wall, so you have to perform the same projection (as you already do) onto the boundary face. The upside of using the turbulence model is that you do not have to implement turbulence model specific shear stress methods, because the turbulence model itself tells, how the shear stress is defined. / Niels |
All times are GMT -4. The time now is 16:05. |