[vitorspadetoventurin]

Hello Friends, My present results of validation blood flow (non-Newtonian) using the model BirdCarreau viscosity (selected in transportProperties) are:

The scientific paper profiles compared with my profiles are:
vertical profile:


horizontal profile:



(the graph above is inverted)

It can be seen in greater rounding my results. The profiles are very different. What's going on? I am missing some parameter or the OpenFOAM is not able to simulate the precision of scientific paper? I found it very far these results: /

The adopted parameters are of a paper where it adopts the following properties:
- Laminar;
- Reynolds < 500;
- Non-newtonian;
- tetraedrical mesh (CHECKMESH gives "Mesh OK");

my mesh:


paper viscosity model:

Note: I just modified the BirdCarreau viscosity model. I changed the value of "a" viscosity (above). OpenFOAM inserts a= 2 as standard and I needed a=0.644 as described above in the paper). The procedure I did was just change the value of "a" constant in BirdCarreau.C, I saved the changes and directory of BirdCarreau.C typed "wmake libso" which gave the following result (which I think is ok):

a@a-Aspire-V3-571:/opt/openfoam230/src/transportModels/incompressible/viscosityModels/BirdCarreau$ wmake libso
wmake: 'Make' directory does not exist in /opt/openfoam230/src/transportModels/incompressible/viscosityModels/BirdCarreau
    Searching up directories tree for Make directory
    Found target directory  ./../..
Making dependency list for source file viscosityModels/BirdCarreau/BirdCarreau.C
SOURCE=viscosityModels/BirdCarreau/BirdCarreau.C ;  g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3  -DNoRepository -ftemplate-depth-100 -I.. -I../twoPhaseMixture/lnInclude -I/opt/openfoam230/src/finiteVolume/lnInclude -IlnInclude -I. -I/opt/openfoam230/src/OpenFOAM/lnInclude -I/opt/openfoam230/src/OSspecific/POSIX/lnInclude   -fPIC -c $SOURCE -o Make/linux64GccDPOpt/BirdCarreau.o
'/opt/openfoam230/platforms/linux64GccDPOpt/lib/libincompressibleTransportModels.so' is up to date.
a@a-Aspire-V3-571:/opt/openfoam230/src/transportModels/incompressible/viscosityModels/BirdCarreau$

Then, my transportProperties file will be:

transportModel  BirdCarreau; 

nu              nu [ 0 2 -1 0 0 0 0 ] 3e-06;


BirdCarreauCoeffs
{
    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 1.560283688e-05; //μ0 = 0.022 Pa.s >> divided by density of 1410 kg/s = 1.56e-05
    nuInf           nuInf [ 0 2 -1 0 0 0 0 ] 1.560283688e-06;//μ∞ = 0.0022 Pa.s >> divided by density of 1410 kg/s = 1.56e-06 
    k               k [ 0 0 1 0 0 0 0 ] 0.11;// k = 0.11 Seconds
    n               n [ 0 0 0 0 0 0 0 ] 0.392; //n= 0.392
}

Here are my other configuration files of the simulated case (USING SIMPLEFOAM SOLVER):
U file:

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet_host
    {
        type            fixedValue;
        value           uniform (0.04284 0 0);
    }

    inlet_graft
    {
        type            fixedValue;
        value           uniform (0.0905 0 -0.0905);
    }

    outlet
    {
        type            zeroGradient;
    }

    wall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

}

p file:

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{

    inlet_host
    {
        type            zeroGradient;
    }

    inlet_graft
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    wall
    {
        type            zeroGradient;
    }

    
}

RASproperties:

RASModel        laminar;

turbulence      off;

printCoeffs     on;

fvschemes:

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss upwind;
    div(phi,k)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div(phi,R)      bounded Gauss upwind;
    div(R)          Gauss linear;
    div(phi,nuTilda) bounded Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p               ;
}

fvsolution file:

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        cacheAgglomeration on;
        agglomerator    faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels     1;
    }

    "(U|k|epsilon|R|nuTilda)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-07;
        relTol          0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;

    residualControl
    {
        p               1e-5;
        U               1e-6;
        "(k|epsilon|omega)" 1e-6;
    }
}

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
        k               0.7;
        epsilon         0.7;
        R               0.7;
        nuTilda         0.7;
    }
}

ControDict file:

application     simpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         2000;

deltaT          0.4;

writeControl    timeStep;

writeInterval   800;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

the simulation converged in 704 iterations.





Where is the error?


Best Regards,
Vitor Spadeto Venturin

[vitorspadetoventurin]

Someone help me? I desire so much to validate that simulation... is a blood flow

[vitorspadetoventurin]

I share with you my case simulation..

http://www.filedropper.com/1p5dq34

[vitorspadetoventurin]

checkMesh gives OK:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           142741
    faces:            1484913
    internal faces:   1412843
    cells:            724439
    faces per cell:   4
    boundary patches: 4
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     0
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    724439
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                  
    inlet_host          430      241      ok (non-closed singly connected)  
    inlet_graft         380      216      ok (non-closed singly connected)  
    wall                70844    35497    ok (non-closed singly connected)  
    outlet              416      235      ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-0.03 -0.0015 -0.0015) (0.03 0.0015 0.0260563)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-4.95962e-17 -4.0039e-17 -3.44229e-16) OK.
    Max cell openness = 2.55702e-16 OK.
    Max aspect ratio = 4.99948 OK.
    Minimum face area = 2.70255e-09. Maximum face area = 5.70468e-08.  Face area magnitudes OK.
    Min volume = 8.05461e-14. Max volume = 4.19327e-12.  Total volume = 6.48978e-07.  Cell volumes OK.
    Mesh non-orthogonality Max: 65.888 average: 14.4021
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.650632 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

[vitorspadetoventurin]

Friends I runned a simulation in order to validate it. I will appreciate if you gives your opinion to improve accuracy. My results already are good, but can improve more through fvschemes and fvsolution.

http://www.filedropper.com/jbiomech2006b
(this is an important paper)

And This is my description and
goals:my_simulation.txt


This is part of my college work.

Thank you for attention

[vitorspadetoventurin]

I forgot the case.. Here is it:
http://www.filedropper.com/showdownload.php/sendfalko

[vonboett]

Hi,
in "Numerical investigation of inertia and shear-thinning effects in axisymmetric flows of carreau fluids by a galerkin least-squares method", Lat. Am. appl. res. from October 2008 they quote Neofytou (2005) to have used the SIMPLE pressure corrector step together with QUICK difference scheme. Maybe that is necessary. I have only one thing to add: As soon as the viscosity is shear-rate dependend, things become strongly grid dependent- so you should make shure that your grid resolution vertical to the flow direction is the same as in your comparison paper. And if your flow is not steady state, it becomes time-step dependent, too, because since your viscosity changes with velocity, and your velocity profile depends on viscosity, you can amplify errors.
One thing popular with turbulence modelling, where the turbulent part of the viscosity makes your viscosity flow dependent, is to introduce an outer loop, "nOuterCorrectors" in the PIMPLE loop in fvSolution dictionary. Maybe that exists for your solver, too. And keeping Courant number below 0.5 may help, too. (controlDict)
I implemented an iterative step, where the viscosity is calculated at the beginning of the time step, and again at the end, and then the time step is performed again with an averaged viscosity, but still dam-break simulations with interFOAM are heavily timestep and grid dependent.
So, I think, you need to get the mesh and timestep the same as in the paper.

[vitorspadetoventurin]

thank for attention friend... I have other case (from the same validation paper)....residuals are +/- (e-06), but in this case... did not converged...... Can you take a look in order to improve results(if possible for you)?

http://www.filedropper.com/forumcase2dq12

really thanks!

(Don't need simulate again.. just take a look in fvschemes and fvsolution.)

[olivierG]

Hello,

Just to say that, if i remeber corectly, you don't have to make a custom carreau viscosity model. In OF carreau model, a=2 by default, but just add a=0.6 in your transportProperties file, and the new "a" value will be taken into account.
PS: use 2nd order scheme also instead of Gauss upwind.

Regards,
olivier