Mesh Generation and Solver Setup

A mesh (grid) divides the flow domain into small cells on which the equations are solved. Its type affects accuracy, memory and effort.

 Structured        Unstructured       Hybrid
 +--+--+--+        /\  /\ /\          ||||||||  prism
 +--+--+--+       /__\/__\/__\        ||||||||  layers
 +--+--+--+       \  /\  /\  /        /\/\/\/\  tetra
 +--+--+--+        \/__\/__\/         core

• Structured: cells are arranged in a regular pattern; each interior node has the same number of neighbours and is addressed by indices $(i,j,k)$. Cells are quadrilaterals (2D) or hexahedra (3D).
• Unstructured: cells (triangles/tetrahedra, or any polyhedra) have arbitrary connectivity, stored in a table; the number of neighbours varies.
• Hybrid: a combination, e.g. thin prism or hexahedral layers next to walls and tetrahedral cells in the core, or structured blocks joined by unstructured regions.

2D elements

• Triangle: three nodes; fits any boundary, fully automatic generation, but needs more cells and gives more numerical diffusion if not aligned with flow.
• Quadrilateral: four nodes; good accuracy when aligned with the flow, fewer cells, preferred in boundary layers and simple shapes.

3D elements

 tetra     hexa        prism       pyramid
   /\      +----+      /|--|\       /\
  /__\     |    |     / |  | \     /  \
  \  /     +----+     \_|__|_/     +---+

Prism (inflation) layers

Prism layers are several thin layers of prismatic (or hexahedral) cells stacked from a wall, with the first-layer height small and each layer thicker than the previous by a growth ratio (usually 1.1 to 1.3).

Reasons for use:

1. Velocity and temperature gradients are steepest normal to the wall, so many cells are needed in that direction only; stretched cells save cell count.
2. The first cell height can be set to meet the required $y^+$ for the turbulence model (about 1 for wall-resolved, 30 to 300 for wall functions).
3. Cells aligned with the wall reduce numerical diffusion and improve wall shear stress and heat flux prediction.
4. Total layer thickness should be about the boundary layer thickness, with smooth transition to the core cells.

Mesh quality strongly affects accuracy, stability and convergence.

Skewness

It measures how far a cell is from an ideal shape (equilateral triangle, rectangle).

• Equiangular skewness:

$$S=\max\left[\frac{\theta_{max}-\theta_e}{180^\circ-\theta_e},\ \frac{\theta_e-\theta_{min}}{\theta_e}\right]$$

where $\theta_e=60^\circ$ for triangles/tetrahedra and $90^\circ$ for quads/hexahedra.

• Range 0 (perfect) to 1 (degenerate). Typical guide: 0 to 0.25 excellent, 0.25 to 0.5 good, 0.5 to 0.8 acceptable, 0.8 to 0.95 poor, above 0.95 unacceptable. Aim for maximum below 0.85 (below 0.9 in 3D tetra).

Aspect ratio

Ratio of the longest to the shortest edge (or dimension) of a cell:

$$AR=\frac{\Delta_{max}}{\Delta_{min}}$$

• Ideal value is 1. High values (10 to 1000) are acceptable only inside boundary layers where gradients are in one direction and the flow is aligned with the long side.
• Elsewhere keep $AR<5$; avoid sudden changes in cell size (size change ratio below about 1.2 to 1.3).

Orthogonality

Measures the angle between the face normal vector and the vector joining the centres of the two adjacent cells (and between the face normal and the vector to the face centre).

• Orthogonal quality $=\min\left(\cos\theta\right)$ where $\theta$ is that angle; range 0 (bad) to 1 (best). Acceptable above about 0.1 (best above 0.2), and non-orthogonality (the angle itself) should be below about 70 to 75°.
• A perfect Cartesian cell has orthogonality 1.

Effects of poor quality

• Large truncation error and numerical diffusion, as the central gradient assumption fails.
• Slow convergence, divergence or negative cell volumes.
• Wrong pressure gradients and forces, even when residuals look small.
• Need for extra non-orthogonal correction loops.

Quality should be checked before solving, and bad cells fixed by re-meshing, smoothing or local refinement.

Equiangular skewness: $S=\max\left[\dfrac{\theta_{max}-\theta_e}{180-\theta_e},\ \dfrac{\theta_e-\theta_{min}}{\theta_e}\right]$, with $\theta_e=60^\circ$ (triangle) and $90^\circ$ (quad). Aspect ratio = longest edge / shortest edge.

Triangle

Edge lengths:

• A(0,0)-B(2,0): $2.000$ mm
• B(2,0)-C(0.5,1.2): $\sqrt{1.5^2+1.2^2}=1.921$ mm
• C-A: $\sqrt{0.5^2+1.2^2}=1.300$ mm

Angles from the dot product:

• At A: vectors AB=(2,0), AC=(0.5,1.2): $\cos\theta=\dfrac{1.0}{2\times1.3}=0.3846$, $\theta_A=67.4^\circ$
• At B: BA=(-2,0), BC=(-1.5,1.2): $\cos\theta=\dfrac{3.0}{2\times1.921}=0.7809$, $\theta_B=38.7^\circ$
• At C: $180-67.4-38.7=74.0^\circ$

$$S=\max\left[\frac{74.0-60}{120},\ \frac{60-38.7}{60}\right]=\max[0.116,\ 0.356]=0.356$$ $$AR=\frac{2.000}{1.300}=1.54$$

Skewness 0.36 falls in the "good" range (0.25 to 0.5), and AR is small, so the triangle is acceptable.

Quadrilateral (a leaning parallelogram)

Edges: bottom $4.0$ mm, right side $\sqrt{0.6^2+0.5^2}=0.781$ mm, top $4.0$ mm, left side $0.781$ mm.

Angles: at (0,0) between (4,0) and (0.6,0.5): $\cos\theta=\dfrac{2.4}{4\times0.781}=0.768$, $\theta=39.8^\circ$. The adjacent angle is $180-39.8=140.2^\circ$.

$$S=\max\left[\frac{140.2-90}{90},\ \frac{90-39.8}{90}\right]=0.558$$ $$AR=\frac{4.0}{0.781}=5.12$$

Comment

The quadrilateral is both strongly skewed and stretched. A stretched cell is fine in a boundary layer if aligned with the wall, but the sheared shape (not rectangular) reduces orthogonality and accuracy; it should be improved by reducing the shear (making the angles closer to 90°).

Answer: Triangle $S=0.36$, $AR=1.54$ (good); quad $S=0.56$, $AR=5.1$ (fair).

Definition: $y^+=\dfrac{\rho u_\tau y}{\mu}$ with $u_\tau=\sqrt{\tau_w/\rho}$. The estimate is done before meshing using flat-plate correlations.

(a) Reynolds number

$$Re_L=\frac{\rho UL}{\mu}=\frac{1.2\times30\times1.5}{1.8\times10^{-5}}=3.0\times10^{6}$$

(turbulent over most of the plate).

(b) Skin friction coefficient

$$C_f=0.026\,(3\times10^6)^{-1/7}=0.026\times0.1188=3.09\times10^{-3}$$

(c) Wall shear stress and friction velocity

$$\tau_w=\tfrac12C_f\rho U^2=0.5\times3.09\times10^{-3}\times1.2\times30^2=1.67\ \text{Pa}$$ $$u_\tau=\sqrt{\frac{\tau_w}{\rho}}=\sqrt{\frac{1.667}{1.2}}=1.179\ \text{m/s}$$

(d) First cell height

Kinematic viscosity $\nu=\mu/\rho=1.5\times10^{-5}$ m²/s. For $y^+=1$ the cell centre is at

$$y=\frac{y^+\nu}{u_\tau}=\frac{1\times1.5\times10^{-5}}{1.179}=1.27\times10^{-5}\ \text{m}$$

The first cell has its centre at half its height, so the first cell height is

$$\Delta y_1=2y=2.55\times10^{-5}\ \text{m}\approx 0.025\ \text{mm}$$

(e) Number of layers

Boundary layer thickness:

$$\delta=0.37\,L\,Re_L^{-1/5}=0.37\times1.5\times(3\times10^6)^{-0.2}=0.0281\ \text{m}$$

For $n$ layers with growth ratio $g=1.2$, total thickness is $\Delta y_1\dfrac{g^n-1}{g-1}\ge\delta$:

$$g^n\ge1+\frac{\delta(g-1)}{\Delta y_1}=1+\frac{0.0281\times0.2}{2.55\times10^{-5}}=221.4$$ $$n\ge\frac{\ln221.4}{\ln1.2}=29.6\ \Rightarrow\ n=30\ \text{layers}$$

Check: with 30 layers the total is $2.55\times10^{-5}\times\dfrac{1.2^{30}-1}{0.2}=0.0301$ m, which is above $\delta$.

Answer: $Re_L=3\times10^6$; $C_f=3.09\times10^{-3}$; $\tau_w=1.67$ Pa; $u_\tau=1.18$ m/s; first cell height $\approx2.5\times10^{-5}$ m; about 30 prism layers.

Solver setup

1. Select the solver type: pressure-based (incompressible and low Mach) or density-based (high-speed compressible); steady or transient.
2. Choose the physics: laminar or turbulence model (e.g. k-$\varepsilon$, k-$\omega$ SST), energy equation on/off, gravity, other models.
3. Materials: density, viscosity, specific heat, conductivity.
4. Boundary conditions: velocity inlet (with turbulence intensity and length scale), pressure outlet, no-slip wall, symmetry, periodic.
5. Discretisation schemes: first-order upwind to start, then second-order upwind or QUICK for accuracy; gradients by least squares or Green-Gauss; pressure scheme (PRESTO or second order).
6. Pressure-velocity coupling: SIMPLE, SIMPLEC or PISO (transient); set under-relaxation factors (e.g. 0.3 for pressure and 0.7 for momentum).
7. Initialise the field (from inlet values or hybrid initialisation) and run.

Convergence criteria

• Residuals (imbalance of each discretised equation) normalised; typically fall by three to four orders of magnitude ($10^{-3}$ continuity/momentum, $10^{-6}$ energy).
• Monitors of quantities of interest (drag, pressure drop, outlet temperature) become steady.
• Global balances: mass flow in equals mass flow out (error below 0.1 to 1%); energy balance closes. Low residuals alone are not proof, so the monitored values must be checked.

Grid independence study

1. Create at least three meshes (coarse, medium, fine) with a uniform refinement ratio $r\ge1.3$.
2. Solve each with identical settings and compare a key output (e.g. pressure drop or $C_D$).
3. When the change between the last two meshes is small (say below 1 to 2%), the solution is considered grid independent; use the coarser of these two.
4. Optionally apply Richardson extrapolation and the Grid Convergence Index (GCI) to estimate the discretisation error.

Finally, validate the results against experimental or analytical data.