263 lines
9.6 KiB
Text
263 lines
9.6 KiB
Text
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "f61d8fbc",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"# AirfRANS finite-volume hand simulation\n",
|
||
|
|
"\n",
|
||
|
|
"Purpose: learn the computation, not the Python stepper or OpenFOAM object graph.\n",
|
||
|
|
"\n",
|
||
|
|
"This notebook reduces the AirfRANS/OpenFOAM momentum loop to the smallest useful stencil:\n",
|
||
|
|
"\n",
|
||
|
|
"1. one internal face shared by two cells;\n",
|
||
|
|
"2. one boundary face on one of those cells.\n",
|
||
|
|
"\n",
|
||
|
|
"Repo grounding:\n",
|
||
|
|
"\n",
|
||
|
|
"- AirfRANS case scale: $U_\\infty=(93.009,6.161,0)$ m/s, $\\nu=1.56\\times10^{-5}$ m²/s, fields `U`, `p`, `phi`, `nut`, `k`, `omega`.\n",
|
||
|
|
"- OpenFOAM v14 momentum predictor source: `OpenFOAM-14/applications/modules/incompressibleFluid/momentumPredictor.C`.\n",
|
||
|
|
"- The real predictor assembles\n",
|
||
|
|
"\n",
|
||
|
|
"$$\n",
|
||
|
|
"\\texttt{fvm::ddt(U)} + \\texttt{fvm::div(phi,U)} + \\texttt{momentumTransport->divDevSigma(U)}\n",
|
||
|
|
"= \\texttt{fvModels().source(U)}\n",
|
||
|
|
"$$\n",
|
||
|
|
"\n",
|
||
|
|
"then solves\n",
|
||
|
|
"\n",
|
||
|
|
"$$\n",
|
||
|
|
"\\texttt{UEqn == -fvc::grad(p)}.\n",
|
||
|
|
"$$\n",
|
||
|
|
"\n",
|
||
|
|
"For the hand calculation, keep one velocity component, ignore time/MRF/source/non-orthogonal details, and replace the viscous/turbulent stress term by scalar diffusion with\n",
|
||
|
|
"\n",
|
||
|
|
"$$\n",
|
||
|
|
"\\Gamma = \\nu + \\nu_t.\n",
|
||
|
|
"$$\n",
|
||
|
|
"\n",
|
||
|
|
"For a cell $P$, use the residual convention\n",
|
||
|
|
"\n",
|
||
|
|
"$$\n",
|
||
|
|
"R_P =\n",
|
||
|
|
"\\sum_f \\underbrace{\\phi_f U_f}_{\\text{convection leaving cell}}\n",
|
||
|
|
"-\n",
|
||
|
|
"\\sum_f \\underbrace{\\Gamma_f \\nabla U_f\\cdot S_f}_{\\text{diffusion leaving cell}}\n",
|
||
|
|
"+\n",
|
||
|
|
"\\sum_f \\underbrace{p_f S_{f,x}}_{\\text{pressure-gradient contribution}}.\n",
|
||
|
|
"$$\n",
|
||
|
|
"\n",
|
||
|
|
"The equation wants $R_P=0$. A face loop is bookkeeping for these three numbers.\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 1,
|
||
|
|
"id": "176670ce",
|
||
|
|
"metadata": {
|
||
|
|
"execution": {
|
||
|
|
"iopub.execute_input": "2026-07-24T06:22:08.071527Z",
|
||
|
|
"iopub.status.busy": "2026-07-24T06:22:08.071431Z",
|
||
|
|
"iopub.status.idle": "2026-07-24T06:22:08.076755Z",
|
||
|
|
"shell.execute_reply": "2026-07-24T06:22:08.076272Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"internal face: P -> N\n",
|
||
|
|
"AirfRANS anchor: U_inf=(93.009, 6.161, 0.0) m/s, nu=1.560e-05 m^2/s\n",
|
||
|
|
"face flux phi_f = 91\n",
|
||
|
|
"upwind U_f = 92\n",
|
||
|
|
"\n",
|
||
|
|
"convection owner P 8372 | neighbour N -8372\n",
|
||
|
|
"diffusion owner P 0.0400312 | neighbour N -0.0400312\n",
|
||
|
|
"pressure owner P 0.15 | neighbour N -0.15\n",
|
||
|
|
"residual owner P 8372.19 | neighbour N -8372.19\n",
|
||
|
|
"\n",
|
||
|
|
"conservation check: R_P + R_N = 0\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"import math\n",
|
||
|
|
"\n",
|
||
|
|
"# Real AirfRANS scale, kept only to anchor the toy numbers.\n",
|
||
|
|
"U_INF = (93.00914503999995, 6.161356013756317, 0.0)\n",
|
||
|
|
"NU = 1.56e-5\n",
|
||
|
|
"\n",
|
||
|
|
"# Toy two-cell state. P is the OpenFOAM owner, N is the neighbour.\n",
|
||
|
|
"# Face area vector S points from P to N.\n",
|
||
|
|
"U_P = 92.0 # x-velocity in owner cell P [m/s]\n",
|
||
|
|
"U_N = 90.0 # x-velocity in neighbour cell N [m/s]\n",
|
||
|
|
"p_P = 0.20 # incompressible pressure p/rho in P [m^2/s^2]\n",
|
||
|
|
"p_N = 0.10 # incompressible pressure p/rho in N [m^2/s^2]\n",
|
||
|
|
"S = 1.0 # face area vector x-component, owner-outward [m^2]\n",
|
||
|
|
"d = 1.0 # owner-to-neighbour center distance [m]\n",
|
||
|
|
"Gamma = NU + 2e-2 # nu + toy turbulent viscosity [m^2/s]\n",
|
||
|
|
"\n",
|
||
|
|
"# 1) Face flux: surfaceScalarField phi = U_f · S_f.\n",
|
||
|
|
"U_for_phi = 0.5 * (U_P + U_N)\n",
|
||
|
|
"phi_f = U_for_phi * S\n",
|
||
|
|
"\n",
|
||
|
|
"# 2) Convection: upwind momentum carried through the face.\n",
|
||
|
|
"U_f = U_P if phi_f >= 0 else U_N\n",
|
||
|
|
"convection_P = phi_f * U_f\n",
|
||
|
|
"convection_N = -convection_P\n",
|
||
|
|
"\n",
|
||
|
|
"# 3) Diffusion: - Gamma * grad(U) · S.\n",
|
||
|
|
"grad_U_dot_S_P = (U_N - U_P) / d * S\n",
|
||
|
|
"diffusion_P = -Gamma * grad_U_dot_S_P\n",
|
||
|
|
"diffusion_N = -diffusion_P\n",
|
||
|
|
"\n",
|
||
|
|
"# 4) Pressure-gradient contribution: p_f * S.\n",
|
||
|
|
"p_f = 0.5 * (p_P + p_N)\n",
|
||
|
|
"pressure_P = p_f * S\n",
|
||
|
|
"pressure_N = -pressure_P\n",
|
||
|
|
"\n",
|
||
|
|
"# 5) Residual contribution from this single internal face.\n",
|
||
|
|
"R_P_internal = convection_P + diffusion_P + pressure_P\n",
|
||
|
|
"R_N_internal = convection_N + diffusion_N + pressure_N\n",
|
||
|
|
"\n",
|
||
|
|
"assert math.isclose(R_P_internal + R_N_internal, 0.0, abs_tol=1e-12)\n",
|
||
|
|
"\n",
|
||
|
|
"def row(label, owner, neighbour):\n",
|
||
|
|
" print(f\"{label:<13} owner P {owner:>12.6g} | neighbour N {neighbour:>12.6g}\")\n",
|
||
|
|
"\n",
|
||
|
|
"print(\"internal face: P -> N\")\n",
|
||
|
|
"print(f\"AirfRANS anchor: U_inf=({U_INF[0]:.3f}, {U_INF[1]:.3f}, {U_INF[2]:.1f}) m/s, nu={NU:.3e} m^2/s\")\n",
|
||
|
|
"print(f\"face flux phi_f = {phi_f:.6g}\")\n",
|
||
|
|
"print(f\"upwind U_f = {U_f:.6g}\")\n",
|
||
|
|
"print()\n",
|
||
|
|
"row(\"convection\", convection_P, convection_N)\n",
|
||
|
|
"row(\"diffusion\", diffusion_P, diffusion_N)\n",
|
||
|
|
"row(\"pressure\", pressure_P, pressure_N)\n",
|
||
|
|
"row(\"residual\", R_P_internal, R_N_internal)\n",
|
||
|
|
"print()\n",
|
||
|
|
"print(f\"conservation check: R_P + R_N = {R_P_internal + R_N_internal:.6g}\")\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 2,
|
||
|
|
"id": "7325d390",
|
||
|
|
"metadata": {
|
||
|
|
"execution": {
|
||
|
|
"iopub.execute_input": "2026-07-24T06:22:08.077905Z",
|
||
|
|
"iopub.status.busy": "2026-07-24T06:22:08.077807Z",
|
||
|
|
"iopub.status.idle": "2026-07-24T06:22:08.081575Z",
|
||
|
|
"shell.execute_reply": "2026-07-24T06:22:08.080994Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"boundary face on P: noSlip wall\n",
|
||
|
|
"patch algebra: U_b=0, p_b=p_P, phi_b=0\n",
|
||
|
|
"\n",
|
||
|
|
"convection owner P 0\n",
|
||
|
|
"diffusion owner P 3.68287\n",
|
||
|
|
"pressure owner P -0.2\n",
|
||
|
|
"wall total owner P 3.48287\n",
|
||
|
|
"\n",
|
||
|
|
"R_P before boundary = 8372.19\n",
|
||
|
|
"R_P after boundary = 8375.67\n",
|
||
|
|
"R_N unchanged = -8372.19\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# Add one boundary face on the left side of owner cell P.\n",
|
||
|
|
"# This mimics the airfoil wall idea: noSlip velocity, zeroGradient pressure, no normal flux.\n",
|
||
|
|
"\n",
|
||
|
|
"U_wall = 0.0 # noSlip\n",
|
||
|
|
"p_wall = p_P # zeroGradient pressure -> boundary value equals owner value\n",
|
||
|
|
"S_wall = -1.0 # owner-outward area vector points left\n",
|
||
|
|
"area_wall = 1.0\n",
|
||
|
|
"d_wall = 0.5 # distance from cell center to wall face\n",
|
||
|
|
"phi_wall = 0.0 # impermeable wall\n",
|
||
|
|
"\n",
|
||
|
|
"# Convection through an impermeable wall is zero.\n",
|
||
|
|
"convection_wall_P = phi_wall * U_wall\n",
|
||
|
|
"\n",
|
||
|
|
"# Boundary diffusion uses the patch value instead of a neighbour value.\n",
|
||
|
|
"grad_U_dot_S_wall = (U_wall - U_P) / d_wall * area_wall\n",
|
||
|
|
"diffusion_wall_P = -Gamma * grad_U_dot_S_wall\n",
|
||
|
|
"\n",
|
||
|
|
"# Pressure term uses the boundary pressure and the boundary area vector.\n",
|
||
|
|
"pressure_wall_P = p_wall * S_wall\n",
|
||
|
|
"\n",
|
||
|
|
"R_wall_P = convection_wall_P + diffusion_wall_P + pressure_wall_P\n",
|
||
|
|
"R_P_after_wall = R_P_internal + R_wall_P\n",
|
||
|
|
"R_N_after_wall = R_N_internal\n",
|
||
|
|
"\n",
|
||
|
|
"assert R_N_after_wall == R_N_internal\n",
|
||
|
|
"assert diffusion_wall_P > 0\n",
|
||
|
|
"\n",
|
||
|
|
"def one(label, value):\n",
|
||
|
|
" print(f\"{label:<13} owner P {value:>12.6g}\")\n",
|
||
|
|
"\n",
|
||
|
|
"print(\"boundary face on P: noSlip wall\")\n",
|
||
|
|
"print(\"patch algebra: U_b=0, p_b=p_P, phi_b=0\")\n",
|
||
|
|
"print()\n",
|
||
|
|
"one(\"convection\", convection_wall_P)\n",
|
||
|
|
"one(\"diffusion\", diffusion_wall_P)\n",
|
||
|
|
"one(\"pressure\", pressure_wall_P)\n",
|
||
|
|
"one(\"wall total\", R_wall_P)\n",
|
||
|
|
"print()\n",
|
||
|
|
"print(f\"R_P before boundary = {R_P_internal:.6g}\")\n",
|
||
|
|
"print(f\"R_P after boundary = {R_P_after_wall:.6g}\")\n",
|
||
|
|
"print(f\"R_N unchanged = {R_N_after_wall:.6g}\")\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"id": "ec2cb6d6",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## What this teaches\n",
|
||
|
|
"\n",
|
||
|
|
"The OpenFOAM loop shape is now visible:\n",
|
||
|
|
"\n",
|
||
|
|
"| notebook arithmetic | OpenFOAM idea |\n",
|
||
|
|
"|---|---|\n",
|
||
|
|
"| `convection_*` | `fvm::div(phi,U)` moves momentum with face flux |\n",
|
||
|
|
"| `diffusion_*` | `momentumTransport->divDevSigma(U)` / diffusion uses face area and cell spacing |\n",
|
||
|
|
"| `pressure_*` | `-fvc::grad(p)` appears as face pressure times area vector |\n",
|
||
|
|
"| internal equal/opposite signs | owner/neighbour addressing conserves internal face transfers |\n",
|
||
|
|
"| wall-only contribution | a boundary patch replaces the missing neighbour with boundary-condition algebra |\n",
|
||
|
|
"\n",
|
||
|
|
"This is not yet SIMPLE, turbulence, pressure correction, or a real airfoil mesh. It is the smallest executable picture of what every finite-volume face loop contributes to the residual.\n",
|
||
|
|
"\n",
|
||
|
|
"Good use: change one value at a time — `U_N`, `p_N`, `Gamma`, `U_wall`, or `phi_wall` — predict the sign of the residual change, then rerun.\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"kernelspec": {
|
||
|
|
"display_name": "Python 3",
|
||
|
|
"language": "python",
|
||
|
|
"name": "python3"
|
||
|
|
},
|
||
|
|
"language_info": {
|
||
|
|
"codemirror_mode": {
|
||
|
|
"name": "ipython",
|
||
|
|
"version": 3
|
||
|
|
},
|
||
|
|
"file_extension": ".py",
|
||
|
|
"mimetype": "text/x-python",
|
||
|
|
"name": "python",
|
||
|
|
"nbconvert_exporter": "python",
|
||
|
|
"pygments_lexer": "ipython3",
|
||
|
|
"version": "3.12.12"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 5
|
||
|
|
}
|