
[Sponsors] 
October 17, 2011, 08:34 
Instability in transport solution on tet mesh

#1 
New Member
Peter Maday
Join Date: Aug 2011
Posts: 14
Rep Power: 7 
Dear forum members,
I have created a solver based on icoFoam to solve for the time dependent concentration distribution of a contrast agent dissolved in the fluid. The aim of the simulation is to recover the concentration values of the agent for applications with blood flow simulation. The considered domain consists of a Y shaped junction, time varying boundary conditions have been specified on the inlet for the velocity and the concentration values describing a periodic change in the concentration and velocity magnitude values. The problem is that when using a tetrahedral mesh generated with a library (vmtk based on TetGen) some abrupt values are obtained for the concentration values, that are clearly wrong (like 50000 when the values should be between 0 and 1) and contrary to the expected smoothness of the distribution the concentration values have discontinuities (see images). The pressure and velocity fields seem to be fine, its only the transported scalar field that seem the be problematic. Based on visual inspection the mesh seems fine (contains uniformly distributed tetrahedral elements), and checkMesh sais its ok. The used timestep is 0.01s. When running the same simulation on the same geometry with a hex mesh (generated with snappyHexMesh) with comparable mesh resolution the simulation seems to run fine. The question is if anyone should have an idea about what could possibly go wrong in this (very simple case)? The boundary condition settings (/0/) and the solver settings are based on the elbow tutorial example (2D pipe flow with junction) + modifications for T (the scalar field that is transported) The mesh: The results: ================================================== ============ Settings: The following boundary conditions have been used: pressure (p) walls zeroGradient inlet zeroGradient outlets fixedValue uniform 0 velocity (U) walls fixedValue uniform (0,0,0) inlet timeVaryingUniformFixedValue outlets zeroGradient transported scalar field (contrast concentration) (T) walls fixedValue; uniform 0; inlet type timeVaryingUniformFixedValue outlets zeroGradient; the schemes for T: div(phi,T) Gauss linear; laplacian(nu,T) Gauss linear corrected; the solver settings for T: solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; the solver for T: solve ( fvm::ddt(T) + fvm::div(phi, T)  fvm::laplacian(nu, T) // nu ); 

October 17, 2011, 10:02 

#2  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 368
Rep Power: 12 
Quote:
div(phi,T) Gauss Gamma 1; If the problem persist you can try also other limited schemes (see the User/Programmers guides or take a look here in the forum about the usage of limited schemes in general). Hope this helps V. PSBTW, what scheme are you using for solving the velocity field? And what about the maxCo value related to your DeltaT? 

October 18, 2011, 07:36 

#3 
New Member
Peter Maday
Join Date: Aug 2011
Posts: 14
Rep Power: 7 
Dear Vesselin,
Thank you very much for the idea, using the Gamma differencing scheme you have suggested seemed to resolve the stability problem for the transported scalar field. To answer your questions, for the velocity the following settings are used: div(phi,U) Gauss limitedLinearV 1; laplacian(nu,U) Gauss linear corrected; laplacian((1A(U)),p) Gauss linear corrected; The settings were taken from the elbow (icoFoam) example. In case the maxCo you referred to means the Courant number, its maximum is larger than 1, around 2.5 most of the time. The mean CN is around 0.2. As far as I know the Courant number value should be less then 1 to gain meaningful results. However further reducing the time step should make the simulation time higher, and my question is if it really is worth the additional effort to warrant that the condition is not violated? Is there any generic principle to say that for example in case the residual stays low enough and the results are believable than the settings should be fine even thought the CN is higher than 1 in some cases? Peter P.S.: The purpose of the simulation is to compute the concentration of contrast agent in the blood flow from which simulated xray projection images (angiograms) can be generated similar to those acquired for diagnostic purposes. The simulation is needed to provide ground truth information about velocity values for image analysis algorithms working with the angiographic data, and thus extreme accuracy in the simulation is not required. 

October 18, 2011, 10:08 

#4  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 368
Rep Power: 12 
Quote:
I Quote:
1) when possible, try to keep the max Co number around 1 also with implicit unsteady solvers like icoFoam (a good way to do it could be introducing adaptive time stepping: you can very easily modify your code looking at the pimpleFoam solver, which has this capability built in). 2) if you want to run yor simulation with large timesteps, you should base your code on the pimpleFoam solver procedure: here, an outer loop correction procedure is introduced which, cycling many times on a single time step, allows you to retain good accuracy and stability with a max Co far larger then 1. Of course, there should be a trade off between the number of outer cycles and the max Co (instead you will not have any time saving), but I can tell you that I was able to run a highRe turbulent unsteady simulation with a max Co of about 3.4 and only one outer iteration, which is approximately the same as running with a max Co of about 1.7 and a single time step integration. Hope this (also) helps V. 

October 18, 2011, 11:05 

#5 
New Member
Peter Maday
Join Date: Aug 2011
Posts: 14
Rep Power: 7 
If you assume your flow to be laminar (based on the expected behavior of the system) than should the ability to use increased time steps in pimpleFoam make the computations faster than when using the simpler icoFoam solver? (pimpleFoam does contain turbulence modeling right?)
Thanks a lot, Peter 

October 18, 2011, 11:25 

#6  
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 368
Rep Power: 12 
Quote:
Best V. 

Tags 
tet mesh instability 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
how to set periodic boundary conditions  Ganesh  FLUENT  13  January 22, 2014 05:11 
3D Hybrid Mesh Errors  DarrenC  ANSYS Meshing & Geometry  11  August 5, 2013 06:42 
Convergence moving mesh  lr103476  OpenFOAM Running, Solving & CFD  30  November 19, 2007 15:09 
Mesh for 3 dim Geometry  Phil  FLUENT  9  July 12, 2000 04:39 