From: Ansys Fluent Rust Computational fluid dynamics

Code: https://github.com/A1i98/rust-cfd-lid-driven-cavity

Bind#

1
[[meta bind projects]]

Tasks#

Sessions `BUTTON[add-session]` `BUTTON[open-latest-session]`#

Simulation in ANSYS#

A two-dimensional square cavity with dimensions of 1 m × 1 m was created in ANSYS DesignModeler to represent the computational domain for the lid-driven cavity simulation.

Lid Driven Cavity Simulation in ANSYS Fluent-1783243704919.webp

After the sketch was converted into a surface body, the resulting geometry was assigned the Fluid body type. This step defines the computational domain as a fluid region, allowing ANSYS Fluent to generate the control volumes required for solving the governing flow equations.

Lid Driven Cavity Simulation in ANSYS Fluent-1783248425439.webp A structured quadrilateral mesh was generated with a global element size of 0.01 m. Local mesh refinement was then applied to improve spatial resolution in regions where higher velocity gradients are expected. Finally, appropriate Named Selections were created for all boundaries, including the moving wall (lid) and stationary walls, to facilitate the assignment of boundary conditions in Fluent.

Lid Driven Cavity Simulation in ANSYS Fluent-1783249314993.webp

For full access, please contact me.

Lid Driven Cavity Simulation in ANSYS Fluent-1783253543118.webp

Simulation in Rust#

Purpose

This note documents a complete 2D lid-driven cavity flow solver written in Rust. The solver uses a finite-difference staggered layout and an artificial-compressibility pseudo-time iteration to obtain the steady incompressible solution. It also writes CSV, Tecplot, and SVG plot files automatically.

Problem Definition#

The classical lid-driven cavity is a square 2D domain of side length $L = 1$ . The top wall moves horizontally with velocity $U = 1$ , while the other three walls are stationary. The flow is assumed incompressible, laminar, Newtonian, and two-dimensional.

The default case solved by the Rust program is:

Parameter	Symbol	Default value
Domain length	$L$	$1$
Lid velocity	$U$	$1$
Reynolds number	$Re$	$100$
Grid size	$N \times N$	$65 \times 65$
Pseudo-time step	$\Delta t$	$0.001$
Artificial-compressibility parameter	$\delta$	$4.5$

Physical meaning

The final result is the steady-state incompressible cavity flow. The artificial-compressibility time variable is a pseudo-time used for convergence, not necessarily the physical transient time.

Governing Equations#

For a 2D incompressible Newtonian flow, the dimensional governing equations are the continuity equation and the Navier-Stokes momentum equations:

\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0

\frac{\partial u}{\partial t} +u\frac{\partial u}{\partial x} +v\frac{\partial u}{\partial y} =-\frac{1}{\rho}\frac{\partial p}{\partial x} +\nu\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right)

\frac{\partial v}{\partial t} +u\frac{\partial v}{\partial x} +v\frac{\partial v}{\partial y} =-\frac{1}{\rho}\frac{\partial p}{\partial y} +\nu\left(\frac{\partial^2v}{\partial x^2}+\frac{\partial^2v}{\partial y^2}\right)

Non-Dimensional Form#

Using $L$ , $U$ , and $\rho U^2$ as reference scales, the equations become:

\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0

\frac{\partial u}{\partial t} +\frac{\partial uu}{\partial x} +\frac{\partial uv}{\partial y} =-\frac{\partial p}{\partial x} +\frac{1}{Re}\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right)

\frac{\partial v}{\partial t} +\frac{\partial uv}{\partial x} +\frac{\partial vv}{\partial y} =-\frac{\partial p}{\partial y} +\frac{1}{Re}\left(\frac{\partial^2v}{\partial x^2}+\frac{\partial^2v}{\partial y^2}\right)

where:

Re=\frac{UL}{\nu}

The Rust code solves the steady incompressible problem through pseudo-time marching. For this reason, the continuity equation is replaced during the iterative process by the artificial-compressibility equation:

\frac{1}{\delta}\frac{\partial p}{\partial t} +\frac{\partial u}{\partial x} +\frac{\partial v}{\partial y}=0

At convergence, $\partial p/\partial t \rightarrow 0$ , so the incompressibility constraint is recovered:

\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\rightarrow 0

Boundary Conditions#

The cavity walls use no-penetration and no-slip boundary conditions.

Boundary	Velocity condition
Top wall / moving lid	$u=U$ , $v=0$
Bottom wall	$u=0$ , $v=0$
Left wall	$u=0$ , $v=0$
Right wall	$u=0$ , $v=0$

The pressure uses a zero-normal-gradient condition at the walls:

\frac{\partial p}{\partial n}=0

Because incompressible pressure is defined up to an arbitrary constant, the code subtracts the mean pressure after every pressure update.

Numerical Method#

The code uses a staggered finite-difference arrangement similar to a MAC grid:

Quantity	Storage idea	Purpose
$u$	horizontal velocity stored on vertical control-volume faces	reduces pressure-velocity decoupling
$v$	vertical velocity stored on horizontal control-volume faces	improves discrete continuity coupling
$p$	pressure-like field at pressure nodes	enforces incompressibility through pseudo-time correction
$u_c$ , $v_c$ , $p_c$	cell-centered post-processed values	exported for plotting and comparison

Why staggered storage is useful

On a collocated grid, pressure and velocity can decouple if the interpolation is not handled carefully. The staggered layout naturally couples pressure gradients and mass fluxes.

Staggered Grid Arrangement#

For a grid size n, the Rust solver stores the fields as follows:

Field	Rust vector size	Logical index range
`u`	`n * (n + 1)`	$i=0..n-1$ , $j=0..n$
`v`	`(n + 1) * n`	$i=0..n$ , $j=0..n-1$
`p`	`(n + 1) * (n + 1)`	$i=0..n$ , $j=0..n$

The grid spacing is:

\Delta x=\Delta y=\frac{1}{n-1}

The exported cell-centered values are computed as:

u_c(i,j)=\frac{1}{2}\left[u(i,j)+u(i,j+1)\right]

v_c(i,j)=\frac{1}{2}\left[v(i,j)+v(i+1,j)\right]

p_c(i,j)=\frac{1}{4}\left[p(i,j)+p(i+1,j)+p(i,j+1)+p(i+1,j+1)\right]

Discrete Equations Used in the Code#

X-Momentum Equation#

The continuous non-dimensional x-momentum equation is:

\frac{\partial u}{\partial t} +\frac{\partial uu}{\partial x} +\frac{\partial uv}{\partial y} =-\frac{\partial p}{\partial x} +\frac{1}{Re}\nabla^2u

The pseudo-time finite-difference update used in the code is:

u^{new}_{i,j}=u_{i,j} -\Delta t\left[ \frac{u^2_{i+1,j}-u^2_{i-1,j}}{2\Delta x} +\frac{(uv)_n-(uv)_s}{\Delta y} +\frac{p_{i+1,j}-p_{i,j}}{\Delta x} \right] +\frac{\Delta t}{Re}\left[ \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2} +\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{\Delta y^2} \right]

where:

(uv)_n=\frac{1}{4}(u_{i,j}+u_{i,j+1})(v_{i,j}+v_{i+1,j})

(uv)_s=\frac{1}{4}(u_{i,j}+u_{i,j-1})(v_{i,j-1}+v_{i+1,j-1})

Y-Momentum Equation#

The y-momentum equation is:

\frac{\partial v}{\partial t} +\frac{\partial uv}{\partial x} +\frac{\partial vv}{\partial y} =-\frac{\partial p}{\partial y} +\frac{1}{Re}\nabla^2v

The pseudo-time update used in the code is:

v^{new}_{i,j}=v_{i,j} -\Delta t\left[ \frac{(uv)_e-(uv)_w}{\Delta x} +\frac{v^2_{i,j+1}-v^2_{i,j-1}}{2\Delta y} +\frac{p_{i,j+1}-p_{i,j}}{\Delta y} \right] +\frac{\Delta t}{Re}\left[ \frac{v_{i+1,j}-2v_{i,j}+v_{i-1,j}}{\Delta x^2} +\frac{v_{i,j+1}-2v_{i,j}+v_{i,j-1}}{\Delta y^2} \right]

where:

(uv)_e=\frac{1}{4}(u_{i,j}+u_{i,j+1})(v_{i,j}+v_{i+1,j})

(uv)_w=\frac{1}{4}(u_{i-1,j}+u_{i-1,j+1})(v_{i,j}+v_{i-1,j})

Artificial-Compressibility Pressure Update#

The discrete pressure update is:

p^{new}_{i,j}=p_{i,j} -\Delta t\,\delta\left[ \frac{u^{new}_{i,j}-u^{new}_{i-1,j}}{\Delta x} +\frac{v^{new}_{i,j}-v^{new}_{i,j-1}}{\Delta y} \right]

The convergence criterion is based mainly on the mean absolute divergence:

\epsilon_{div}=\frac{1}{N_{int}}\sum_{i,j} \left| \frac{u_{i,j}-u_{i-1,j}}{\Delta x} +\frac{v_{i,j}-v_{i,j-1}}{\Delta y} \right|

The default stopping criterion is:

\epsilon_{div}<10^{-8}

Rust Code#

Configuration and CLI Arguments#

The solver can be controlled from the terminal without editing the source file.

1
#[derive(Debug, Clone)]
2
struct Config {
3
    n: usize,
4
    re: f64,
5
    lid_velocity: f64,
6
    dt: f64,
7
    delta: f64,
8
    max_steps: usize,
9
    min_steps: usize,
10
    tolerance: f64,
11
    report_every: usize,
12
    out_dir: PathBuf,
13
    stream_sor_omega: f64,
14
    stream_max_iter: usize,
15
    stream_tol: f64,
16
}

The most important parameters are n, re, dt, delta, and tolerance.

Staggered Indexing#

The arrays are one-dimensional Vec<f64> buffers, but the code maps them to 2D staggered indices.

1
fn iu(&self, i: usize, j: usize) -> usize {
2
    // u: i = 0..n-1, j = 0..n
3
    j * self.cfg.n + i
4
}
5

6
fn iv(&self, i: usize, j: usize) -> usize {
7
    // v: i = 0..n, j = 0..n-1
8
    j * (self.cfg.n + 1) + i
9
}
10

11
fn ip(&self, i: usize, j: usize) -> usize {
12
    // p: i = 0..n, j = 0..n
13
    j * (self.cfg.n + 1) + i
14
}

X-Momentum Update#

This section applies convection, pressure gradient, and viscous diffusion for the horizontal velocity.

1
let u_ij = self.u[self.iu(i, j)];
2

3
let duu_dx = (self.u[self.iu(i + 1, j)].powi(2)
4
    - self.u[self.iu(i - 1, j)].powi(2))
5
    / (2.0 * dx);
6

7
let uv_n = 0.25
8
    * (self.u[self.iu(i, j)] + self.u[self.iu(i, j + 1)])
9
    * (self.v[self.iv(i, j)] + self.v[self.iv(i + 1, j)]);
10
let uv_s = 0.25
11
    * (self.u[self.iu(i, j)] + self.u[self.iu(i, j - 1)])
12
    * (self.v[self.iv(i, j - 1)] + self.v[self.iv(i + 1, j - 1)]);
13
let duv_dy = (uv_n - uv_s) / dy;
14

15
let dp_dx = (self.p[self.ip(i + 1, j)] - self.p[self.ip(i, j)]) / dx;
16

17
let lap_u = (self.u[self.iu(i + 1, j)] - 2.0 * u_ij + self.u[self.iu(i - 1, j)]) / dx2
18
    + (self.u[self.iu(i, j + 1)] - 2.0 * u_ij + self.u[self.iu(i, j - 1)]) / dy2;
19

20
let new_u = u_ij - dt * (duu_dx + duv_dy + dp_dx) + dt * inv_re * lap_u;

Y-Momentum Update#

This section is the analogous update for vertical velocity.

1
let uv_e = 0.25
2
    * (self.u[self.iu(i, j)] + self.u[self.iu(i, j + 1)])
3
    * (self.v[self.iv(i, j)] + self.v[self.iv(i + 1, j)]);
4
let uv_w = 0.25
5
    * (self.u[self.iu(i - 1, j)] + self.u[self.iu(i - 1, j + 1)])
6
    * (self.v[self.iv(i, j)] + self.v[self.iv(i - 1, j)]);
7
let duv_dx = (uv_e - uv_w) / dx;
8

9
let dvv_dy = (self.v[self.iv(i, j + 1)].powi(2)
10
    - self.v[self.iv(i, j - 1)].powi(2))
11
    / (2.0 * dy);
12

13
let dp_dy = (self.p[self.ip(i, j + 1)] - self.p[self.ip(i, j)]) / dy;
14

15
let lap_v = (self.v[self.iv(i + 1, j)] - 2.0 * v_ij + self.v[self.iv(i - 1, j)]) / dx2
16
    + (self.v[self.iv(i, j + 1)] - 2.0 * v_ij + self.v[self.iv(i, j - 1)]) / dy2;
17

18
let new_v = v_ij - dt * (duv_dx + dvv_dy + dp_dy) + dt * inv_re * lap_v;

Pressure and Continuity Update#

The pressure equation drives the divergence toward zero.

1
for i in 1..=(n - 1) {
2
    for j in 1..=(n - 1) {
3
        let divergence = (self.un[self.iu(i, j)] - self.un[self.iu(i - 1, j)]) / dx
4
            + (self.vn[self.iv(i, j)] - self.vn[self.iv(i, j - 1)]) / dy;
5
        let idx = self.ip(i, j);
6
        self.pn[idx] = self.p[idx] - dt * delta * divergence;
7
    }
8
}

Boundary Conditions#

The top lid is imposed through a ghost-cell relation. For the moving lid:

\frac{u_{top\,ghost}+u_{top\,inside}}{2}=U

therefore:

u_{top\,ghost}=2U-u_{top\,inside}

1
for i in 0..n {
2
    let bottom = self.iu(i, 0);
3
    let bottom_inside = self.iu(i, 1);
4
    let top = self.iu(i, n);
5
    let top_inside = self.iu(i, n - 1);
6
    self.un[bottom] = -self.un[bottom_inside];
7
    self.un[top] = 2.0 * self.cfg.lid_velocity - self.un[top_inside];
8
}

Post-Processing Fields#

The solver computes vorticity as:

\omega=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}

It also computes a streamfunction from:

\nabla^2\psi=-\omega

The post-processing step writes the final fields, centerline profiles, and plots.

How to Run#

1. Build in Release Mode#

1
cargo build --release

2. Run the Default Case#

1
cargo run --release -- --n 65 --re 100 --dt 0.001 --delta 4.5 --out results

Lid Driven Cavity Simulation-1783534125409.webp

3. Quick Test Case#

For a fast check before the final run:

1
cargo run --release -- --n 33 --re 100 --max-steps 30000 --out quick_results

4. Useful Parameter Changes#

Goal	Suggested command option
Finer grid	`--n 129`
Lower Reynolds number	`--re 10`
Higher Reynolds number	`--re 400`
More strict convergence	`--tol 1e-9`
Faster rough run	`--tol 1e-6 --n 33`
Change output folder	`--out my_results`

Stability warning

For larger Reynolds numbers or finer grids, reduce --dt. A safe first adjustment is to halve dt when the residual oscillates or grows.

Output Files#

After running, the solver creates an output folder such as results/.

File	Description
`run_config.txt`	Simulation parameters and final step
`convergence.csv`	Iteration history: divergence and velocity change
`field.csv`	Full field data: `x,y,u,v,p,speed,vorticity,streamfunction`
`field_tecplot.dat`	Tecplot-compatible field file
`centerline_u_vertical.csv`	$u$ profile along $x=0.5$
`centerline_v_horizontal.csv`	$v$ profile along $y=0.5$
`images/u_velocity.svg`	SVG heatmap for $u$
`images/v_velocity.svg`	SVG heatmap for $v$
`images/speed.svg`	SVG heatmap for velocity magnitude
`images/pressure.svg`	SVG heatmap for pressure
`images/vorticity.svg`	SVG heatmap for vorticity
`images/streamfunction.svg`	SVG heatmap for streamfunction
`images/centerline_u_vertical.svg`	Centerline $u$ profile plot
`images/centerline_v_horizontal.svg`	Centerline $v$ profile plot
`images/convergence.svg`	Convergence history plot

Plot#

Velocity Magnitude#

Lid Driven Cavity Simulation-1783534250339.webp

Figure 1. Velocity magnitude field.

Horizontal Velocity Contour#

Lid Driven Cavity Simulation-1783534310294.webp

Figure 2. Horizontal velocity field.

Vertical Velocity Contour#

Lid Driven Cavity Simulation-1783534316065.webp

Figure 3. Vertical velocity field.

Pressure Field#

Lid Driven Cavity Simulation-1783534330264.webp

Figure 4. Pressure field from artificial-compressibility iteration.

Vorticity Field#

Lid Driven Cavity Simulation-1783534341066.webp

Figure 5. Vorticity field.

Streamfunction Field#

Lid Driven Cavity Simulation-1783534349156.webp

Figure 6. Streamfunction field.

Vertical Centerline U Profile#

Lid Driven Cavity Simulation-1783534376199.webp

Figure 7. $u$ velocity profile along $x=0.5$ .

Horizontal Centerline V Profile#

Lid Driven Cavity Simulation-1783534383976.webp

Figure 8. $v$ velocity profile along $y=0.5$ .

Convergence History#

Lid Driven Cavity Simulation-1783534392017.webp

Figure 9. Convergence history based on mean absolute divergence.

Verification and Comparison with ANSYS Fluent#

For full access, please contact me.

Numerical Stability Notes#

The explicit pseudo-time update is simple and transparent, but stability depends on grid size, Reynolds number, and pseudo-time step.

A conservative starting point is:

\Delta t \lesssim \min\left(0.5\frac{\Delta x}{U},\,0.25Re\Delta x^2\right)

For the default case:

\Delta x=\frac{1}{64}=0.015625

0.25Re\Delta x^2 \approx 0.0061

The default $\Delta t=0.001$ is therefore conservative for $Re=100$ .

Higher-Reynolds-number warning

At higher Reynolds numbers, the cavity develops thinner boundary layers and stronger gradients. Use a finer grid and reduce dt. For example, try --n 129 --dt 0.0005 as a starting point for $Re=400$ .

References#

Ghia, U., Ghia, K. N., & Shin, C. T. (1982). High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics, 48(3), 387-411.
Peyret, R., & Taylor, T. D. (1983). Computational Methods for Fluid Flow. Springer.
Ferziger, J. H., Perić, M., & Street, R. L. (2020). Computational Methods for Fluid Dynamics. Springer.
Patankar, S. V. (1980). Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing.

Bind#

Tasks#

Sessions BUTTON[add-session] BUTTON[open-latest-session]#

Simulation in ANSYS#

Simulation in Rust#

Problem Definition#

Governing Equations#

Non-Dimensional Form#

Boundary Conditions#

Numerical Method#

Staggered Grid Arrangement#

Discrete Equations Used in the Code#

X-Momentum Equation#

Y-Momentum Equation#

Artificial-Compressibility Pressure Update#

Rust Code#

Configuration and CLI Arguments#

Staggered Indexing#

X-Momentum Update#

Y-Momentum Update#

Pressure and Continuity Update#

Boundary Conditions#

Post-Processing Fields#

How to Run#

1. Build in Release Mode#

2. Run the Default Case#

3. Quick Test Case#

4. Useful Parameter Changes#

Output Files#

Plot#

Velocity Magnitude#

Horizontal Velocity Contour#

Vertical Velocity Contour#

Pressure Field#

Vorticity Field#

Streamfunction Field#

Vertical Centerline U Profile#

Horizontal Centerline V Profile#

Convergence History#

Verification and Comparison with ANSYS Fluent#

Numerical Stability Notes#

References#

Pages mentioning this page

Sessions `BUTTON[add-session]` `BUTTON[open-latest-session]`#