V&V » ChoCG: Poiseuille flow

This example uses ChoCG in Inciter to compute the stationary viscous constant-density (incompressible) laminar flow field in a channel due to an imposed pressure gradient. Since the analytic solution is known the objective is to establish and verify the accuracy of the numerical method and the correctness of the software implementation.

Problem setup

The problem domain consists of a rectangular channel with an imposed streamwise pressure gradient. While the flow field is nominally 2D, we compute it in 3D. The length of the channel in the X coordinate direction is 20 units, its height in Y is $L_y=1$ , and the width in Z is adjusted to contain a single element in the Z direction depending on the mesh resolution.

An example mesh is displayed below. The initial and boundary conditions prescribe a pressure gradient of $\text{d}p/\text{d}x=-0.12$ , thus the analytic solution for the streamwise velocity is

\[ v_x = \frac{1}{2\mu}\left(-\frac{\text{d}p}{\text{d}x}\right)y\left(L_y-y\right) \]

where $p$ is the pressure and $\mu=0.01$ is the dynamic viscosity which yields the laminar Reynolds number as $ \text{Re} = \rho \overline{v}_x L_y / \mu = 100 $ , defined based on the fluid density $\rho=1$ and the average streamwise velocity $ \overline{v}_x = 1 $ . The flow starts with quiescent initial conditions: $\mbox{\boldmath$v$}(\mbox{\boldmath$x$},t=0)=(0,0,0)$ . At the Z extremes of the computational domain Dirichlet boundary conditions are set to $v_y=v_z=0$ and no-slip/no-penetration boundary conditions, $v_x=v_y=v_z=0$ , are applied at the Y extremes representing the channel walls.

Image Coarse surface mesh for computing Poiseuille flow, nelem=3K, npoin=1.2K.

Code revision to reproduce

To reproduce the results below, use code revision 291b6f1 and the control file below.

Control file

-- vim: filetype=lua:

print "Poiseuille flow"

-- mesh: poiseuille?tetz.exo

term = 300.0
ttyi = 100

cfl = 0.7

solver = "chocg"
flux = "damp2"

fct = false

part = "rcb"

mat = { dyn_viscosity = 0.01 }

momentum = {
  iter = 50,
  tol = 1.0e-3,
  pc = "jacobi"
}

pressure = {
  iter = 500,
  tol = 1.0e-2,
  pc = "jacobi",
  bc_dir = { { 1, 2 },
             { 2, 2 } },
  bc_dirval = { { 1, 2.4 },
                { 2, 0.0 } }
}

ic = {
  velocity = { 0.0, 0.0, 0.0 }
}

bc_noslip = {
  sideset = { 3, 4 }
}

bc_dir = {
  { 5, 0, 1, 1 }
}

fieldout = {
  iter = 1000,
  sideset = { 1, 2, 3 }
}

diag = {
  iter = 1,
  format = "scientific",
  precision = 12
}

Run on a single CPU

Main/inciter -i poiseuille3tetz.exo -c poiseuille_chocg.q -r 10000

Numerical solutions and errors

To estimate the order of convergence, the numerical solution is computed using three different meshes, whose properties of are tabulated below.

MeshPointsTetrahedrah
11,2123,0000.25
24,42212,0000.16
316,84248,0000.08

Here h is the average edge length.

Plotted below are the computed streamwise velocity profiles sampled in the middle of the domain along a line in the Y direction using three different meshes, obtained after the solution has converged to a stationary state.

Image

The L1 and L2 errors are computed as

\[ \begin{split} \left| \left| \varepsilon \right| \right|_1 & = \frac{\sum_{i=1}^n V^i \left| \hat{v}_x^i - v_x^i \right|} {\sum_{i=1}^n V^i} \\ \left| \left| \varepsilon \right| \right|_2 & = \left[\frac{\sum_{i=1}^n V^i \left( \hat{v}_x^i - v_x^i \right)^2} {\sum_{i=1}^n V^i}\right]^{1/2} \end{split} \]

where $n$ is the number of points sampled across the channel, $\hat{v}_x^i$ and $v_x^i$ are the exact and computed solutions at mesh point $i$ , and $V^i$ denotes the volume associated to mesh point $i$ .

The tables below display the numerical errors for two different types of damping for the advection term. The order of convergence is estimated using the two finer meshes as

\[ p = \frac{ \log \left| \left| \varepsilon \right| \right|^{m+1} - \log \left| \left| \varepsilon \right| \right|^m}{ \log h^{m+1} - \log h^m} \]

where $m$ is a mesh index.

Errors and convergence rates using second-order damping for advection.

MeshhL1(err)L2(err)
10.250.05981390.0636658
20.160.01333450.0144788
30.080.00237680.0027099
p2.4882.4176

Errors and convergence rates using fourth-order damping for advection.

MeshhL1(err)L2(err)
10.250.0555480.0592916
20.160.0132520.0143954
30.080.0017350.0018838
p2.9332.934