[zhang-qiang]

I want to implement a simulation of blood flow in artery.


I copied configure file from TestCases/incomp_navierstokes/bend, and the modified part is:

INC_DENSITY_INIT= 998.2 % use the water density
MU_CONSTANT= 3E-3 % use the blood viscosity
MARKER_INLET= ( inlet1, 288.15, 0.1, -8.932349, 0.446989, -4.473627 )
MARKER_SYM= ( NONE )
ITER= 1000
MESH_FILENAME= VesselCFD.su2

However, the result is very weird. The configure & SU2 and result are uploaded in github: https://github.com/zhang-qiang-github/VesselSU2.

How should I modify the configure file to achieve a correct result?

Any suggestion is appreciated~~~

[bigfootedrockmidget]

Can you post some pictures of your mesh, especially at the inlet and outlet? It does not look good when I visualize it. Compare this to the mesh from the 90 degree bend.

[zhang-qiang]

The mesh is:



The left picture is my mesh, and the right picture is the result.

In the left picture, the yellow color indicates wall, the red color indicates inlet, and blue color indicates the outlet.

The mesh and configure file have been uploaded in github: https://github.com/zhang-qiang-github/VesselSU2.

The right picture is generated by flow.vtu with VTK:

import vtkmodules.all as vtk

reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName('flow.vtu')
reader.Update()

flow = reader.GetOutput()

mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(flow)

mapper.SetScalarModeToUsePointFieldData()
mapper.ScalarVisibilityOn()
mapper.SelectColorArray('Pressure')
mapper.Modified()

actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetOpacity(0.3)

renderer = vtk.vtkRenderer()
renderer.AddActor(actor)
renderer.SetBackground(1, 1, 1)
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(renderer)
iren = vtk.vtkRenderWindowInteractor()
iren.SetInteractorStyle(vtk.vtkInteractorStyleTrac kballCamera())
iren.SetRenderWindow(renWin)
iren.Initialize()
iren.Start()

[bigfootedrockmidget]

Can you show a picture that shows the surface edges of the individual cells that make up your model? Like in Figure 1 or Figure 2 of the 90 degree bend tutorial?
https://su2code.github.io/tutorials_...nshot_gmsh.png

[zhang-qiang]

The surface is:



The tetrahedron is:

[bigfootedrockmidget]

If you look at the surface mesh at the inlet and outlet, then you can see that the mesh quality is not good. Have a look at the mesh in the 90 degree Bend and try to reproduce something like that, with a layer of finer cells near the wall to capture the velocity gradients. From the center to the wall you need at least 10 cells to capture the solution with some decent accuracy if it's laminar, more for turbulence.

[zhang-qiang]

I have optimized the mesh of inlet and outlet. The complete mesh is:



The inlet/outlet is:





However, the result is different from the previous one, but it still is not correct:



The optimized SU2 file has been uploaded to github: https://github.com/zhang-qiang-github/VesselSU2

Is this inlet/outlet still not good enough to obtain a good SU2 result?

[bigfootedrockmidget]

In the second/third figure, the one of the inlet and outlet, you see that a lot of cells seem to converge to a point close to the top. This leads to a very bad mesh quality. How do you make this mesh? Is there a method for you to check the mesh quality, especially the cell skewness? Can't you gmsh for this geometry?

You have to make this surface mesh such that the cells are evenly distributed over the surface, and preferably that the cells are smaller close to the wall.
I made the 90 degree bend by defining a square region in the middle of the bend, and then I defined 4 additional squares where one edge was the edge of the middle square and the other edge was 1/4 of the edge of the pipe. This is easy to do in e.g. gmsh.

[zhang-qiang]

Many thanks for the kindly reply.

The initial inlet/outlet mesh is generated by VTK and it is optimized by gmsh.

The initial inlet has very poor quality, and it make the optimized mesh not good.

I will try to generate mesh from only mesh.

[zhang-qiang]

I have generated new mesh.

The inlet is:


The outlet is:


And the result is:


The new generated SU2 file has been updated in github: https://github.com/zhang-qiang-github/VesselSU2

[bigfootedrockmidget]

The surface at the inlet looks better, although it is still very coarse. But if I now look at the mesh at a cross-section through the pipe, the quality of the internal cells still looks bad. Is it possible for you to generate the mesh by extruding the mesh from the inlet to the outlet?

How did you generate this shape?

[zhang-qiang]

Which cross-section has poor quality? I have check the internal cell, and I don't find very poor cell. There are two cross section of my internal cell:





Actually, the mesh should be the surface of artery, and it should be generated from the real medical image.

The currently mesh is generated from a surface with gmsh. The surface is a constructed tube.

There are some printed information:


------------------- Geometry Preprocessing ( Zone 0 ) -------------------
Three dimensional problem.
11725 grid points.
48555 volume elements.
3 surface markers.
14210 boundary elements in index 0 (Marker = wall0).
440 boundary elements in index 1 (Marker = inlet1).
122 boundary elements in index 2 (Marker = outlet2).
48555 tetrahedra.
Setting point connectivity.
Renumbering points (Reverse Cuthill McKee Ordering).
Recomputing point connectivity.
Setting element connectivity.
Checking the numerical grid orientation.
All volume elements are correctly orientend.
There has been a re-orientation of 14210 TRIANGLE surface elements.
Identifying edges and vertices.
Setting the control volume structure.
Volume of the computational grid: 2.13539.
Searching for the closest normal neighbors to the surfaces.
Storing a mapping from global to local point index.
Compute the surface curvature.
Max K: 156.04. Mean K: 23.9922. Standard deviation K: 32.7136.
Checking for periodicity.
Computing mesh quality statistics for the dual control volumes.
+--------------------------------------------------------------+
| Mesh Quality Metric| Minimum| Maximum|
+--------------------------------------------------------------+
| Orthogonality Angle (deg.)| 45.1882| 86.3042|
| CV Face Area Aspect Ratio| 1.24601| 582.774|
| CV Sub-Volume Ratio| 1| 268.881|
+--------------------------------------------------------------+
Finding max control volume width.
Semi-span length = 1.3992 m.
Wetted area = 15.1148 m^2.
Area projection in the x-plane = 4.09077 m^2, y-plane = 3.60007 m^2, z-plane = 4.29287 m^2.
Max. coordinate in the x-direction = 1.18879 m, y-direction = 1.3992 m, z-direction = 4.77619 m.
Min. coordinate in the x-direction = -1.39968 m, y-direction = -1.18908 m, z-direction = -0.357651 m.
Computing wall distances.

-------------------- Solver Preprocessing ( Zone 0 ) --------------------
Incompressible flow: rho_ref, vel_ref, and temp_ref
are based on the initial values, p_ref = rho_ref*vel_ref^2.
Force coefficients computed using initial values.
The reference area for force coeffs. is 1 m^2.
The reference length for force coeffs. is 1 m.
The pressure is decomposed into thermodynamic and dynamic components.
The initial value of the dynamic pressure is 0.
Mach number: 0.00838426, computed using the Bulk modulus.
For external flows, the initial state is imposed at the far-field.
Angle of attack (deg): 0, computed using the initial velocity.
Side slip angle (deg): 0, computed using the initial velocity.
Reynolds number per meter: 33273.3, computed using initial values.
Reynolds number is a byproduct of inputs only (not used internally).
SI units only. The grid should be dimensional (meters).
No energy equation.

-- Models:
+------------------------------------------------------------------------------+
| Viscosity Model| Conductivity Model| Fluid Model|
+------------------------------------------------------------------------------+
| CONSTANT_VISCOSITY| CONSTANT_PRANDTL| CONSTANT_DENSITY|
+------------------------------------------------------------------------------+
-- Fluid properties:
+------------------------------------------------------------------------------+
| Name| Dim. value| Ref. value| Unit|Non-dim. value|
+------------------------------------------------------------------------------+
| Viscosity| 0.003| 99.82| N.s/m^2| 3.00541e-05|
+------------------------------------------------------------------------------+
| Prandtl (Lam.)| -| -| -| 0.72|
| Prandtl (Turb.)| -| -| -| 0.9|
+------------------------------------------------------------------------------+
| Bulk Modulus| 142000| 1| Pa| 142000|
+------------------------------------------------------------------------------+
-- Initial and free-stream conditions:
+------------------------------------------------------------------------------+
| Name| Dim. value| Ref. value| Unit|Non-dim. value|
+------------------------------------------------------------------------------+
| Dynamic Pressure| 0| 9.982| Pa| 0|
| Total Pressure| 4.991| 9.982| Pa| 0.5|
| Density| 998.2| 998.2| kg/m^3| 1|
| Velocity-X| 0.1| 0.1| m/s| 1|
| Velocity-Y| 0| 0.1| m/s| 0|
| Velocity-Z| 0| 0.1| m/s| 0|
| Velocity Magnitude| 0.1| 0.1| m/s| 1|
+------------------------------------------------------------------------------+
| Viscosity| 0.003| 99.82| N.s/m^2| 3.00541e-05|
| Conductivity| -| 0.00346417| W/m^2.K| -|
+------------------------------------------------------------------------------+
| Mach Number| -| -| -| 0.00838426|
| Reynolds Number| -| -| -| 33273.3|
+------------------------------------------------------------------------------+
Initialize Jacobian structure (Navier-Stokes). MG level: 0.

------------------- Numerics Preprocessing ( Zone 0 ) -------------------

----------------- Integration Preprocessing ( Zone 0 ) ------------------

------------------- Iteration Preprocessing ( Zone 0 ) ------------------
Euler/Navier-Stokes/RANS fluid iteration.

------------------------------ Begin Solver -----------------------------

Simulation Run using the Single-zone Driver
+----------------------------------------------------------------+
| Inner_Iter| rms[P]| rms[U]| CL| CD|
+----------------------------------------------------------------+
| 0| -2.943554| -3.294951| -0.627246| 35.153973|
| 1| -3.194206| -2.665360| 0.034259| 19.225236|
| 2| -3.284977| -3.024991| 0.197407| 4.931505|
| 3| -3.412528| -3.359970| 0.033803| -0.194493|
| 4| -3.780346| -3.550017| -0.101823| -1.377322|
| 5| -4.046881| -3.863619| -0.121923| -1.404524|
| 6| -4.090989| -3.998806| -0.088121| -1.221971|
| 7| -4.121224| -4.033508| -0.029000| -1.064993|
| 8| -4.140167| -4.058693| 0.022296| -1.003519|
| 9| -4.153915| -4.096703| 0.048060| -0.991977|
| 10| -4.160987| -4.139719| 0.063729| -0.974949|
| 11| -4.169387| -4.205363| 0.080532| -0.943486|
| 12| -4.169842| -4.229793| 0.106870| -0.906462|
| 13| -4.170700| -4.221808| 0.139401| -0.883823|
| 14| -4.171313| -4.205070| 0.173042| -0.887454|
| 15| -4.170217| -4.176839| 0.204751| -0.918347|
| 16| -4.166617| -4.126737| 0.234515| -0.971228|
| 17| -4.164631| -4.091859| 0.262378| -1.038613|
| 18| -4.163657| -4.065459| 0.290341| -1.111412|
| 19| -4.163199| -4.042905| 0.320506| -1.182465|
| 20| -4.162446| -4.017410| 0.354236| -1.248463|
| 21| -4.161615| -3.992102| 0.391467| -1.309480|
| 22| -4.160992| -3.969860| 0.431208| -1.367421|
| 23| -4.160463| -3.949180| 0.472442| -1.424667|
| 24| -4.159998| -3.929486| 0.514464| -1.483277|
| 25| -4.159588| -3.910588| 0.556993| -1.544606|
| 26| -4.159226| -3.892413| 0.600089| -1.609229|
| 27| -4.158908| -3.874907| 0.644002| -1.677063|
| 28| -4.158632| -3.858026| 0.689027| -1.747607|
| 29| -4.158399| -3.841744| 0.735405| -1.820196|
| 30| -4.158210| -3.826051| 0.783279| -1.894216|
| 31| -4.158068| -3.810945| 0.832694| -1.969237|
| 32| -4.157972| -3.796415| 0.883625| -2.045048|
| 33| -4.157876| -3.781966| 0.936036| -2.121672|
| 34| -4.157791| -3.767736| 0.989858| -2.199300|
| 35| -4.157757| -3.754148| 1.045013| -2.278154|
| 36| -4.157763| -3.741057| 1.101476| -2.358395|
| 37| -4.157808| -3.728416| 1.159241| -2.440092|
| 38| -4.157886| -3.716199| 1.218313| -2.523235|
| 39| -4.157996| -3.704380| 1.278693| -2.607766|
| 40| -4.158133| -3.692937| 1.340377| -2.693596|
| 41| -4.158297| -3.681846| 1.403360| -2.780632|
| 42| -4.158485| -3.671088| 1.467635| -2.868793|
| 43| -4.158696| -3.660645| 1.533198| -2.958023|
| 44| -4.158928| -3.650500| 1.600045| -3.048286|
| 45| -4.159181| -3.640639| 1.668171| -3.139570|
| 46| -4.159454| -3.631052| 1.737570| -3.231872|
| 47| -4.159744| -3.621724| 1.808235| -3.325192|
| 48| -4.160053| -3.612646| 1.880157| -3.419529|
| 49| -4.160378| -3.603807| 1.953325| -3.514875|
| 50| -4.160718| -3.595197| 2.027729| -3.611215|

Is there some abnormal thing for the printed information?

The CL is grow to 180 for 999 iteration. Is it normal?

[bigfootedrockmidget]

Your mesh quality needs to be much better. Your othogonality range is OK, but you have a range of very small and very large cells, and it looks like sometimes small and large cells are next to each other.

Try to make the cells have more or less the same size, and make them at least twice as small as what you have now. If possible, try to create cells that are smaller close to the wall.



Your Reynolds number is also very high, so you need to use one of the turbulence models.

[zhang-qiang]

Thanks for warning of Reynolds number.

Re=ρvd/μ. In my cofigure file:

ρ=998.2 kg/m^3, INC_DENSITY_INIT= 998.2 Kg/m^3
v =0.1m/s, INC_VELOCITY_INIT= ( 0.1, 0.0, 0.0 )
d is about 1
μ=3e-3. MU_CONSTANT= 3E-3
Thus, the Re=998.2*0.1*1/(3e-3)=33000, which is close to the printed Reynolds number: "Reynolds number per meter: 33273.3".

Is my calculation correct?

If the above calculation correct, I need to adjust the input.
The d=1 is wrong for artery, and it should be about 5 cm = 5*10^-2 m. In addition, the velocity should be 30~40 cm/s for blood flow. If d=0.05, v=0.3 and the other parameters keep the same, the Re should be about 5000. Is it a suitable number for laminar model? If not, which parameter should be adjusted?

v=30 cm/s comes from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4058969/:
The usual normal velocity of the common carotid artery is 30-40 cm/sec [19]

MU_CONSTANT= 3E-3 comes from https://en.wikipedia.org/wiki/Hemorheology:
In pascal-seconds (Pa·s), the viscosity of blood at 37 °C is normally 3 × 10−3 to 4 × 10−3,[8]

On the other hand, I would try to generate good mesh with small cell for inlet/outlet.

[bigfootedrockmidget]

Hi,


Reynolds nr in a pipe is usually turbulent when Re>1800, so you will need a turbulence model. But a diameter of 5 cm sounds too large for a human artery, or should it be 0.005 m = 5mm?

[bigfootedrockmidget]

Do you know about this paper:
https://gmsh.info/doc/preprints/gmsh_bio2_preprint.pdf

[zhang-qiang]

Actually, the diameter of aorta may be 3cm, the carotid artery diameter may be 1 cm.

Any I find a reference: https://www.sciencedirect.com/topics...turbulent-flow.

It said:

In the circulatory system, the Reynolds number is 3,000 (mean value) and 7,500 (peak value) for the aorta, 500 for a typical artery, 0.001 in a capillary, and 400 for a typical vein.

Does it means the aorta should be a turbulence model?

I will learn about the paper: https://gmsh.info/doc/preprints/gmsh_bio2_preprint.pdf. And I will carefully generate my mesh.

[zhang-qiang]

Is Re=ρvd/μ? Does the d meaning the diameter of inlet? But, when I change the diameter of tube, the printed Reynolds number do not change. Does the diameter of inlet impact the Reynolds number?