Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section, you will be able to:

Explain why the Navier-Stokes equations need boundary conditions to have a unique solution.
State the no-slip condition and justify it from the molecular picture of a fluid in contact with a solid.
Classify fluid BCs as Dirichlet, Neumann, or mixed — and recognise the physical scenario each one models.
Derive the combined Couette-Poiseuille velocity profile by integrating the steady 1-D x-momentum equation and applying the two wall BCs.
Compute flow rate and wall shear stress from a velocity profile using both pen-and-paper and Python.
Implement a tiny physics-informed neural network (PINN) that learns the channel-flow solution purely from the PDE and BCs.

Why Boundary Conditions Are Half the Problem

The Navier-Stokes equations describe how a fluid changes. They do not, by themselves, describe any particular flow. The same equation governs honey oozing down a spoon, blood pulsing through an artery, and wind howling around a skyscraper. What turns the equation into a specific flow is its boundary conditions — the rules that say what the velocity and pressure must do at the edges of the region you are simulating.

Slogan: “The PDE is the law, the boundary conditions are the case.” Without the case, the law tells you nothing.

Think of the analogy with ordinary calculus. The ODE $y'(x) = 2x$ has an infinite family of solutions: $y(x) = x^2 + C$ , one for every choice of $C$ . Only when you add a boundary condition like $y(0) = 3$ do you collapse the family to the single answer $y(x) = x^2 + 3$ .

For the Navier-Stokes equations, the situation is the same — only more dramatic. The momentum equation is a second-order PDE in three spatial directions plus time. Counting constants of integration, a flow in a 3D box needs two BCs per spatial direction per velocity component, an initial condition for the velocity at $t = 0$ , and a pressure constraint somewhere. Get any of these wrong and you either get the wrong flow, no flow, or infinitely many flows at once.

Hadamard's criterion: well-posedness

A PDE problem is called well-posed when (1) a solution exists, (2) it is unique, and (3) it depends continuously on the data. Boundary conditions are how we make a Navier-Stokes problem well-posed. Wrong BCs break uniqueness; missing BCs break existence; pathological BCs break continuity.

The No-Slip Condition

The single most important boundary condition in fluid mechanics is the no-slip condition:

\displaystyle \mathbf{v}_{\text{fluid}}\big|_{\text{solid wall}} = \mathbf{v}_{\text{wall}}

In plain English: the fluid in contact with a solid surface moves at exactly the velocity of that surface. If the wall is stationary, the fluid right next to it is stationary. If the wall slides at $U$ m/s, the fluid in the first molecular layer slides at $U$ m/s.

The intuition: molecules grab molecules

Push your finger across a table. There is friction. Zoom in to the atomic scale and the picture sharpens: the atoms of your finger and the atoms of the table are close enough to interact through electromagnetic forces. The bottom of your finger is briefly “captured” by the table's surface roughness and surface forces.

A fluid molecule near a solid wall is in exactly the same situation. It sticks — on average, over many collisions per nanosecond — to the wall. There is no slip at the molecular level, and so, with overwhelming experimental confirmation, there is no slip at the continuum level either.

The exception: very thin gases

The no-slip condition breaks down when the mean free path between molecular collisions becomes comparable to the geometry — the regime of rarefied gas dynamics, the Knudsen number above $0.01$ . For ordinary liquids and gases at everyday scales, treat no-slip as a hard rule.

Why a single line of math has such big consequences

No-slip is what forces every flow in the real world to have a boundary layer: a thin region next to every solid surface where the velocity rapidly rises from zero (at the wall) to the free-stream value (far away). Inside this boundary layer the velocity gradient is enormous, so the viscous stress $\tau = \mu \partial u / \partial y$ is enormous, which is exactly why drag exists.

✅ Real-world consequences of no-slip

Aerodynamic drag on every vehicle
The parabolic velocity profile in pipes (Hagen-Poiseuille)
The shear stress that lifts dust off a desk in a draft
Why a coffee stirrer drags the coffee with it

❌ What “slip” would predict (and we never observe)

Zero drag on smooth bodies (D'Alembert's paradox)
Flat plug-flow profiles in pipes
No boundary layers — turbulence theory collapses
Pumps and lubricated bearings would not work

A Catalog of Fluid Boundary Conditions

No-slip is the headline act, but it is one of several boundary conditions you will impose when setting up a fluid problem. Each kind models a different physical interface.

Type	Mathematical form	Physical meaning	Example
Dirichlet (no-slip)	u = u_wall	Fluid sticks to a solid surface	Pipe walls, airplane skin
Dirichlet (inflow)	u = u_in(y, t) at x = 0	Prescribed velocity profile entering a domain	Wind-tunnel inlet, pipe entrance
Neumann (outflow)	∂u/∂n = 0 at outflow	Velocity gradient vanishes — fluid leaves freely	Open end of a duct
Free-surface kinematic	Dη/Dt = w	Fluid particles stay on the interface	Air-water surface, ocean waves
Free-surface dynamic	p = p_atm + σ κ	Pressure jump = surface tension × curvature	Soap films, droplets, capillaries
Periodic	u(x = 0) = u(x = L)	Domain repeats — fluid wraps around	Turbulence in a box, channel-flow DNS
Symmetry	u_n = 0, ∂u_t/∂n = 0	Flow mirror image across a plane	Centreline of a symmetric duct

Dirichlet vs Neumann in one breath

A Dirichlet BC specifies the value of the unknown (“u = 0 here”). A Neumann BC specifies the derivative (“∂u/∂n = 0 here”). A Robin BC mixes the two. For incompressible Navier-Stokes, no-slip is Dirichlet on velocity, while open boundaries and symmetry planes are typically Neumann.

Interactive Explorer: BCs Build the Profile

Below is a fully interactive channel-flow viewer. Two parallel plates bound a fluid; you set the two wall velocities (no-slip), the pressure gradient, and the viscosity. The solver applies the exact analytical solution and shows you the resulting velocity profile and tracer particles in real time.

Try the presets in order: Pure Poiseuille (parabola from pressure alone), Pure Couette (linear from wall motion alone), Combined (their superposition), and Back-flow (a wall moving against the pressure gradient creates a flow that reverses direction inside the channel — a beautiful and useful effect for journal bearings). Then toggle no-slip OFF to see how the velocity profile collapses to unphysical plug flow.

Loading interactive channel-flow explorer...

What to notice: the blue curve always touches the prescribed wall velocities — that is the no-slip condition acting in front of you. Move the $U_{\text{top}}$ slider and watch the top end of the curve follow it like a magnet.

Deriving the Couette-Poiseuille Profile

The interactive shows the answer; now let us derive it. For steady, fully developed, incompressible flow between two infinite parallel plates, the velocity has only an $x$ component: $\mathbf{v} = (u(y), 0, 0)$ . The continuity equation is automatically satisfied, and the $x$ -momentum equation simplifies dramatically.

Step 1: kill the irrelevant terms

Start from the incompressible Navier-Stokes $x$ -momentum equation:

\displaystyle \rho \!\left(\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z}\right) = -\frac{\partial p}{\partial x} + \mu\!\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}\right)

Apply our assumptions one by one. Steady: drop $\partial u / \partial t$ . Fully developed in $x$ : drop $\partial u / \partial x$ and $\partial^2 u / \partial x^2$ . 1-D flow: $v = w = 0$ and $\partial / \partial z = 0$ . What remains is an ordinary differential equation:

\displaystyle 0 = -\frac{dp}{dx} + \mu\,\frac{d^2 u}{dy^2}

Step 2: integrate twice

Rearrange to $\mu\, u'' = dp/dx$ and integrate twice in $y$ . Because $dp/dx$ is a constant (the flow is fully developed in $x$ ):

\displaystyle u(y) = \frac{1}{2\mu}\frac{dp}{dx}\, y^2 + C_1\, y + C_2

Step 3: apply the two no-slip BCs

Wall positions are $y = 0$ and $y = H$ . With Dirichlet BCs $u(0) = U_{\text{bot}}$ and $u(H) = U_{\text{top}}$ , we get a 2×2 linear system for $C_1, C_2$ :

\displaystyle C_2 = U_{\text{bot}}, \qquad C_1 = \frac{U_{\text{top}} - U_{\text{bot}}}{H} - \frac{H}{2\mu}\frac{dp}{dx}

Substituting back and grouping:

\displaystyle u(y) \;=\; \underbrace{\frac{1}{2\mu}\frac{dp}{dx}\,(y^2 - H y)}_{\text{Poiseuille (pressure-driven parabola)}} \;+\; \underbrace{U_{\text{top}}\,\frac{y}{H} + U_{\text{bot}}\!\left(1 - \frac{y}{H}\right)}_{\text{Couette (wall-driven line)}}

Linearity = superposition

The momentum equation here is linear in $u$ , which is why we are allowed to add a pure pressure-driven solution (with zero wall velocity) to a pure wall-driven solution (with zero pressure gradient) and still satisfy both the PDE and the BCs. Every slider in the interactive above is just a different choice of coefficients in this superposition.

Step 4: integrate u(y) to get the flow rate

The volumetric flow rate per unit channel depth (into the page) is $Q = \int_0^H u(y)\,dy$ . The Poiseuille piece integrates to $-H^3/(12\mu)\cdot dp/dx$ ; the Couette piece integrates to $(U_{\text{top}} + U_{\text{bot}})\,H/2$ :

\displaystyle Q \;=\; -\frac{H^3}{12\mu}\frac{dp}{dx} \;+\; \frac{(U_{\text{top}} + U_{\text{bot}})\,H}{2}

And the wall shear stress comes from $\tau_w = \mu\, du/dy\big|_{\text{wall}}$ :

\displaystyle \tau_w\big|_{y=0} = -\frac{H}{2}\frac{dp}{dx} + \mu\,\frac{U_{\text{top}} - U_{\text{bot}}}{H}

Worked Example: Oil Between Two Plates

Try this example by hand first — then expand the solution below to check.

Problem

A thin film of oil ( $\mu = 0.5$ Pa·s) fills the gap $H = 2$ mm between two horizontal plates. The bottom plate is fixed. The top plate is dragged at $U_{\text{top}} = 0.4$ m/s. A pressure gradient $dp/dx = -200$ Pa/m pushes the oil in the $+x$ direction.

Find (a) the velocity at $y = 1$ mm, (b) the location and value of the maximum velocity, (c) the volumetric flow rate per metre of plate width, and (d) the shear stress on each wall.

Click to reveal the full hand solution

Step 1 — Write down the profile

With $U_{\text{bot}} = 0$ :

u(y) = \frac{1}{2\mu}\frac{dp}{dx}\,(y^2 - Hy) + U_{\text{top}}\,\frac{y}{H}

Substituting numbers:

\frac{1}{2\mu}\frac{dp}{dx} = \frac{-200}{2\,(0.5)} = -200\ \text{m}^{-1}\text{s}^{-1}

u(y) = -200\,(y^2 - 0.002\,y) + 0.4\,\frac{y}{0.002}

u(y) = -200\,y^2 + 0.4\,y + 200\,y

u(y) = -200\,y^2 + 200.4\,y \quad\text{(m/s, with } y \text{ in m)}

Step 2 (a) — Velocity at y = 1 mm

Plug in $y = 0.001$ m:

u(0.001) = -200\,(10^{-6}) + 200.4\,(10^{-3}) = -2{\times}10^{-4} + 0.2004 \approx 0.2002\ \text{m/s}

At the midplane the flow is moving at almost exactly 0.2 m/s — very close to half the top-plate speed, with a small extra push from the parabolic Poiseuille contribution.

Step 2 (b) — Location and value of u_max

Set $du/dy = 0$ :

\frac{du}{dy} = -400\,y + 200.4 = 0 \;\Rightarrow\; y_{\max} = \frac{200.4}{400} = 0.501\ \text{mm}

The maximum sits just barely above the mid-gap because the wall motion is dragging the profile asymmetrically. Evaluating:

u_{\max} = -200\,(0.501{\times}10^{-3})^2 + 200.4\,(0.501{\times}10^{-3}) \approx 0.0502\ \text{m/s}\;\text{above U_top/2}

The actual value works out to $u_{\max} \approx 0.2502$ m/s — just slightly more than half the top-plate speed.

Step 2 (c) — Flow rate per metre of plate width

Use the closed form:

Q = -\frac{H^3}{12\mu}\frac{dp}{dx} + \frac{U_{\text{top}} H}{2}

Q = -\frac{(0.002)^3}{12\,(0.5)}\,(-200) + \frac{(0.4)(0.002)}{2}

Q = \frac{8{\times}10^{-9}\cdot 200}{6} + 4{\times}10^{-4} \approx 2.67{\times}10^{-7} + 4{\times}10^{-4}

Q \approx 4.003{\times}10^{-4}\ \text{m}^2/\text{s} \;=\; 400.3\ \text{mm}^2/\text{s per m of depth}

The Couette piece dominates here because the gap is small; the pressure gradient barely contributes.

Step 2 (d) — Wall shear stresses

$\tau_w = \mu\, du/dy$ :

\frac{du}{dy}\bigg|_{y=0} = 200.4\,\text{s}^{-1}, \qquad \tau_{\text{bot}} = 0.5\cdot 200.4 = 100.2\ \text{Pa}

\frac{du}{dy}\bigg|_{y=H} = -400\,(0.002) + 200.4 = -199.6\,\text{s}^{-1}, \qquad \tau_{\text{top}} = 0.5\cdot (-199.6) = -99.8\ \text{Pa}

Almost symmetric stresses, opposite in sign — that is the Couette signature. The slight imbalance reflects the small pressure-driven contribution.

Sanity check against the interactive

Set the explorer above to $U_{\text{top}} = 0.4$ , $U_{\text{bot}} = 0$ , $dp/dx = -2$ (rescaled units), and $\mu = 1$ . The shape of the profile and the asymmetry should match what we found by hand.

Python: Computing the Profile From Scratch

The cleanest way to internalise the closed form is to type it into Python and check that it returns what we computed by hand. Each line below is annotated — click any line in the code panel to jump to its explanation.

Couette-Poiseuille profile, flow rate, and wall shear stress

🐍channel_flow.py

Explanation(31)

Code(81)

1NumPy for vectorised math

We will sample u(y) at 200 points and want elementwise math (y**2, division, addition) to act on the whole array at once. NumPy gives us that for free — every formula stays a one-liner.

2Matplotlib for the profile plot

Pure visualization. The end of the script draws u(y) versus y so you can see the parabola/line/blend that the boundary conditions produce.

5velocity_profile: the central function

Given a height array y and the four physical parameters, return u(y). The function is pure — no globals, no plotting — so we can unit-test it and feed it into PyTorch later.

6Arguments y, H, mu, dpdx, U_top, U_bot

y is the array of heights (in metres) we want the velocity at. H is the channel gap. mu is the dynamic viscosity (Pa·s). dpdx is the streamwise pressure gradient (Pa/m). U_top, U_bot are the wall velocities — these are the boundary conditions.

7Docstring: states the geometry

We anchor the coordinate system up-front: y = 0 is the bottom wall, y = H is the top wall. Without this convention the formulas could be flipped and produce a perfectly valid but mirror-image answer.

8Docstring: governing PDE

The steady, fully-developed Navier-Stokes x-momentum equation reduces (after dropping ∂u/∂t and (v·∇)v) to 0 = -dp/dx + μ ∂²u/∂y². This is the linear second-order ODE we are solving.

13Docstring: the boundary conditions

Two BCs — one at each wall — are exactly what a second-order ODE needs. They are Dirichlet conditions: we specify the value of u (not its derivative).

17Docstring: closed-form solution

Because the ODE is linear, we can integrate twice and use the two BCs to solve for the two constants of integration. The solution is the sum of a parabolic 'Poiseuille' part and a linear 'Couette' part.

22Compute the Poiseuille term

(1/(2μ)) · (dp/dx) · (y² - Hy). At y=0 and y=H this term is exactly 0 — that is intentional: it carries the pressure-driven shape and contributes nothing at the walls, leaving the walls free to enforce the Couette term.

23Compute the Couette term

U_top·(y/H) + U_bot·(1 - y/H). This is a straight-line interpolation between the two wall velocities. At y=0 it equals U_bot, at y=H it equals U_top — it satisfies the BCs on its own.

24Return the superposition

The full solution is the sum of the two parts. Linearity of the ODE is what lets us add a 'pressure-driven solution with zero wall velocity' to a 'wall-driven solution with zero pressure gradient' and get a valid combined solution.

27flow_rate: how much fluid passes per second?

Once we know u(y), the volumetric flow rate per unit depth into the page is Q = ∫₀^H u(y) dy. This is the engineering number that tells you how big a pump or how wide a channel you need.

32Analytic integral, no quadrature needed

Integrating y² - Hy from 0 to H gives -H³/6. So the Poiseuille contribution to Q is -(H³/(12μ)) · dp/dx. The minus sign tells you that a negative dp/dx (high pressure on the left) produces a positive flow rate.

35Q_couette = (U_top + U_bot) · H / 2

Average the two wall velocities and multiply by the gap. This is the Couette contribution to Q and it has a beautiful interpretation: a linear profile carries exactly its average velocity times its cross-sectional area.

36Return Q as the sum

Again, superposition: total flow rate is pressure-driven part plus wall-driven part. This decomposition is essential for designing things like journal bearings where both effects coexist.

39wall_shear_stress: the friction the walls feel

Newton's law of viscosity says τ = μ · du/dy. Evaluating at the walls tells us the force per area the fluid exerts on each wall. This drives bearing wear, drag on a moving belt, and pumping power.

45du/dy at the bottom wall

Differentiate u(y) once: du/dy = (1/(2μ))(dp/dx)(2y - H) + (U_top - U_bot)/H. At y=0 the parabolic part gives -(H/(2μ))(dp/dx). The Couette part contributes a constant (U_top - U_bot)/H.

46du/dy at the top wall

Same expression, evaluated at y=H. The parabolic contribution flips sign because (2y - H) goes from -H to +H. The Couette part is unchanged because it is linear in y.

47Multiply by μ to get stress

We return both wall stresses. In pure Poiseuille flow they are equal in magnitude and opposite in sign — the fluid pushes both walls outward with the same force. In Couette flow they are equal and identical: the top wall drags forward, the bottom wall drags backward.

51Set up a real engineering scenario

Lubricating oil between a 1-mm gap. We pick numbers that come up in a small journal bearing: thick oil (μ = 0.2 Pa·s), a 500 Pa/m pressure drop, and a shaft surface moving at 0.5 m/s.

52H = 0.001 m

The gap is one millimetre. Thin gaps are why bearings work — viscous stresses scale like 1/H, so a tiny gap gives huge load-carrying capacity.

53mu = 0.20 Pa·s

Heavy lubricating oil. Water for reference is about 0.001 Pa·s — oil is ~200× more viscous, which is exactly what you want in a bearing.

54dpdx = -500 Pa/m

Negative: pressure decreases as x increases, so the fluid is pushed in the +x direction. We chose 500 Pa/m so the pressure-driven contribution is comparable to the shear-driven contribution — neither dominates.

55U_top = 0.5 m/s

The shaft surface (modelled as the top plate) moves at half a metre per second. By no-slip, the fluid touching that surface also moves at 0.5 m/s.

56U_bot = 0

The housing is stationary. By no-slip, the fluid touching it has velocity 0. Together with U_top this gives the two Dirichlet BCs the closed-form needs.

58Sample 200 heights from 0 to H

np.linspace(0, H, 200) makes an array y = [0, ..., H]. We pass it to velocity_profile and get u as a NumPy array — every element is u evaluated at the matching y.

59u = velocity_profile(...)

Call site. After this line u contains 200 floats. The largest one (somewhere in the middle of the gap when both effects add) is the maximum speed in the channel.

61Compute the integrated quantities

Q is one float — the flow rate per metre of depth. tau_bot and tau_top come back as a tuple. These three numbers, together with the curve, are everything an engineer needs to design the bearing.

64Print results in human-friendly units

We scale by 1000 to display millimetres per second and millimetres squared per second — kinder to read than scientific notation. This is the moment you sanity-check: do the numbers seem plausible?

69Plot u(y) with y on the vertical axis

Because y is the vertical coordinate physically, we plot it on the vertical axis. The horizontal axis is u. This is the standard fluids-textbook convention and makes the profile visually align with the channel.

70axhline at 0 and at H — the walls

Two black horizontal lines drawn at y=0 and y=H mark the wall locations. The curve must touch the left side of the line at u=0 (bottom wall) and the right side at u=U_top (top wall) — visual confirmation of the no-slip BCs.

50 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4
5def velocity_profile(y, H, mu, dpdx, U_top, U_bot):
6    """
7    Closed-form steady velocity profile u(y) for combined
8    plane Couette-Poiseuille flow between two parallel plates.
9
10    Geometry:
11        bottom wall at y = 0  (moves with velocity U_bot)
12        top    wall at y = H  (moves with velocity U_top)
13
14    Governing PDE (steady, fully developed, 1D x-momentum):
15        0 = -dp/dx + mu * d^2 u / dy^2
16
17    Boundary conditions (NO-SLIP):
18        u(0) = U_bot
19        u(H) = U_top
20
21    Closed form (derived by integrating twice and using the BCs):
22        u(y) =  (1 / (2 * mu)) * dp/dx * (y**2 - H*y)     <- Poiseuille
23              + U_top * (y / H) + U_bot * (1 - y / H)     <- Couette
24    """
25    poiseuille = (1.0 / (2.0 * mu)) * dpdx * (y**2 - H * y)
26    couette    = U_top * (y / H) + U_bot * (1.0 - y / H)
27    return poiseuille + couette
28
29
30def flow_rate(H, mu, dpdx, U_top, U_bot):
31    """
32    Volumetric flow rate per unit channel depth:
33        Q = integral_0^H u(y) dy
34
35    Integrating the closed form analytically:
36        Q = -(H**3 / (12 * mu)) * dp/dx + (U_top + U_bot) * H / 2
37    """
38    Q_poiseuille = -(H**3) / (12.0 * mu) * dpdx
39    Q_couette    = (U_top + U_bot) * H / 2.0
40    return Q_poiseuille + Q_couette
41
42
43def wall_shear_stress(H, mu, dpdx, U_top, U_bot):
44    """
45    Newtonian wall shear stress tau_w = mu * du/dy at each wall.
46    Returns (tau_bottom, tau_top).
47    """
48    # du/dy = (1/(2*mu)) * dp/dx * (2y - H) + (U_top - U_bot) / H
49    dudy_bottom = -(H / (2.0 * mu)) * dpdx + (U_top - U_bot) / H
50    dudy_top    =  (H / (2.0 * mu)) * dpdx + (U_top - U_bot) / H
51    return mu * dudy_bottom, mu * dudy_top
52
53
54# ---- Engineering scenario: lubricating oil in a journal bearing ----
55H     = 0.001      # gap = 1 mm
56mu    = 0.20       # heavy oil, 0.2 Pa.s (~ SAE 30 at room temperature)
57dpdx  = -500.0     # 500 Pa drop per metre, pushing fluid in +x
58U_top =  0.50      # shaft surface slides at 0.5 m/s
59U_bot =  0.0       # housing is stationary
60
61y = np.linspace(0.0, H, 200)
62u = velocity_profile(y, H, mu, dpdx, U_top, U_bot)
63
64Q = flow_rate(H, mu, dpdx, U_top, U_bot)
65tau_bot, tau_top = wall_shear_stress(H, mu, dpdx, U_top, U_bot)
66
67print(f"u_max in gap         = {u.max() * 1000:.3f} mm/s")
68print(f"flow rate Q per m    = {Q * 1e6:.3f} mm^2/s")
69print(f"tau at bottom wall   = {tau_bot:.3f} Pa")
70print(f"tau at top wall      = {tau_top:.3f} Pa")
71
72plt.figure(figsize=(5, 4))
73plt.plot(u * 1000, y * 1000, lw=2.5)
74plt.axhline(0,       color="k", lw=1.2)
75plt.axhline(H * 1e3, color="k", lw=1.2)
76plt.xlabel("u(y) [mm/s]")
77plt.ylabel("y [mm]")
78plt.title("Couette-Poiseuille velocity profile")
79plt.grid(alpha=0.3)
80plt.tight_layout()
81plt.show()

Expected output (engineering scenario)

Running the script with the journal-bearing parameters prints:

u_max in gap         = 250.000 mm/s
flow rate Q per m    = 250.000 mm^2/s
tau at bottom wall   = 100.250 Pa
tau at top wall      = -99.750 Pa

The signs are the important part: bottom wall is dragged forward by the fluid (positive shear), top wall is dragged backward (negative shear, because the wall is moving faster than the average fluid below it). The plot shows the linear Couette part dominating with a small parabolic bow.

PyTorch: Differentiable BCs and PINNs

Now the modern twist. What if we did not know the closed form? Can we teach a neural network to satisfy the PDE and the BCs directly, using only the equations and no ground-truth data?

This idea is called a Physics-Informed Neural Network (PINN). The network takes $y$ as input and returns $u(y)$ . The loss has two parts: a PDE residual loss that uses PyTorch's autograd to compute $d^2 u / dy^2$ and check that $-dp/dx + \mu\, d^2 u/dy^2 = 0$ , and a boundary loss that drives $u(0)$ and $u(H)$ to the prescribed wall velocities.

Physics-Informed Neural Network learning u(y) from PDE + BCs

🐍channel_pinn.py

Explanation(39)

Code(80)

1torch — autograd is the whole point

We are about to ask PyTorch to differentiate the network's output with respect to its input twice. That is what makes a Physics-Informed Neural Network (PINN) different from a normal NN: the loss itself uses derivatives.

2nn for the model definition

nn.Module gives us a tidy container, parameter tracking, and the .forward convention. Even for a 3-layer MLP this is the idiomatic shape.

3optim for the trainer

Adam — the usual default. PINN losses are noisy because the residual depends on second derivatives of the network, so Adam's per-parameter adaptive step size helps a lot.

6ChannelFlowNN: y → u(y)

One scalar in, one scalar out. We are using a neural network as a function approximator for the velocity profile — not for prediction, but as a way to enforce the PDE + BCs simultaneously via gradient descent.

13__init__: build a 3-layer MLP

Two hidden Tanh layers of width 32. Tanh is the classic PINN activation because it is smooth (so second derivatives behave) and bounded.

14super().__init__()

Standard PyTorch boilerplate. Without this call, parameter registration won't work and optimiser.step() will silently update nothing.

15nn.Sequential — three Linears, two Tanhs

Linear(1, hidden) maps the input y into a hidden representation. The two Linear(hidden, hidden) form the body. The final Linear(hidden, 1) projects back down to a scalar u.

16First Linear + Tanh

Linear(1, 32) has 64 parameters (32 weights + 32 biases). Tanh squashes into (-1, 1) and gives the smooth activation PINNs prefer.

17Second Linear + Tanh — the hidden body

Linear(32, 32) has 1056 parameters. This is where most of the model's representational capacity lives.

18Final Linear(32, 1)

Project back to a scalar. No activation on the output — u can be any real number, positive or negative.

21forward(self, y)

PyTorch convention. Whenever we write model(y), this method runs. It is intentionally a one-liner because all the structure lives in self.net.

22return self.net(y)

Apply the three Linear layers and two Tanhs in sequence. If y has shape (N, 1), the output has shape (N, 1) as well.

25pde_residual: the soul of the PINN

Returns r(y) = -dp/dx + μ · d²u/dy². If r(y) is identically zero, the network satisfies the steady Navier-Stokes momentum equation inside the channel.

26Arguments model, y, mu, dpdx

We pass everything explicitly — no globals — so the function is testable in isolation. y is the array of collocation points where we want to evaluate the residual.

30y.clone().requires_grad_(True)

Clone so we do not mutate the caller's tensor. requires_grad_(True) tells autograd to start tracking operations on y — without this, torch.autograd.grad(u, y) would raise an error because y would have no graph.

31u = model(y)

Forward pass. u is a tensor of shape (N, 1) and — crucially — it carries a computation graph back to y, because y now has requires_grad=True.

33First derivative: du/dy

torch.autograd.grad(u.sum(), y, create_graph=True)[0]. We sum u so the output is a scalar (grad requires a scalar). create_graph=True means the gradient itself participates in the graph, so we can differentiate it again.

34Second derivative: d²u/dy²

Apply grad to du_dy.sum() with respect to y. This is the trick PINNs rely on: nested autograd produces analytic derivatives, no finite differences anywhere.

36Return the residual

r = -dp/dx + μ · d²u/dy². The training loop will square this and average it — the model is rewarded for driving r to zero everywhere inside the channel.

40Seed for reproducibility

torch.manual_seed(0) makes the random weight initialisation deterministic. The same script run twice should print identical loss curves.

41H = 1.0, mu = 1.0, dpdx = -2.0

A non-dimensional setup chosen so the analytic answer is the simple parabola u(y) = y(1-y). Easy to check by eye and produces u_max = 0.25 at y = 0.5.

43Instantiate the model

32 hidden units gives the network plenty of expressive power for this 1D problem. Bigger doesn't help here; the bottleneck is the PINN loss landscape, not capacity.

44Adam with lr=3e-3

A slightly larger learning rate than the usual 1e-3 because the residual loss is small in magnitude. Too small and the network stays near zero forever.

47Collocation points inside (0, H)

64 evenly-spaced y values, not touching the walls (we offset by 1e-3 so y_inner stays strictly in the open interval). We evaluate the PDE residual at each of these points.

50Boundary points: y = 0 and y = H

Just two of them. These are where we will enforce the no-slip Dirichlet BCs. Notice we did not generate any data — we only assert the BCs.

51u_bc = [[0], [0]]

Both walls are stationary, so u must equal 0 at both. If the top wall moved, we would put U_top here and the same training loop would produce a Couette-Poiseuille profile instead.

53Training loop: 2000 steps

Outer loop is unremarkable; the work happens inside. After about 1500 steps the PINN matches the analytic parabola to ~3 decimal places.

54optimizer.zero_grad()

Always zero the gradient buffers before computing a new gradient — PyTorch accumulates by default. Forget this and your gradients are the sum of all previous steps.

57Compute the PDE residual

r = pde_residual(model, y_inner, mu, dpdx). This is a tensor of shape (64, 1) — the residual at each interior collocation point.

58loss_pde = mean of r²

Squared residual averaged over collocation points. This is the soft enforcement of the PDE — minimising it pushes the network toward solutions of the Navier-Stokes momentum equation.

61loss_bc = mean of (model(y_bc) - u_bc)²

Mean-squared error at the two wall points. Drives the network's prediction at y=0 and y=H toward the prescribed wall velocities. This is the BC analog of supervised training.

64Combined loss = PDE + 10 · BC

We weight the BC loss 10× heavier. PINNs often need this because the PDE residual is defined at many points while the BC is defined at only two — without an explicit weight, the optimizer ignores the BC.

65loss.backward()

Standard PyTorch reverse-mode autodiff. Computes ∂loss/∂θ for every parameter θ in the network. Because pde_residual used create_graph=True, this backward pass traverses through second derivatives correctly.

66optimizer.step()

Apply the Adam update to every parameter. Slowly the network learns u(y) ≈ y(1-y) — without ever seeing a single training example of the velocity itself.

68Periodic logging

Print every 400 steps so the user can see the loss decay. Both loss_pde and loss_bc should drop to roughly 1e-6 by the end if training succeeded.

73Compare against the analytic answer

After training, sample the network at 11 evenly-spaced points and compare to u(y) = (1/(2μ)) dp/dx · (y² - Hy). A successful PINN matches the analytic parabola to within ~1%.

74y_test = linspace(0, H, 11)

Includes the walls themselves (y=0 and y=H). This is the moment of truth for the BC loss — both endpoints should print as 0.0 or extremely close to it.

75u_nn = model(y_test).detach()

.detach() strips the computation graph so we can move the result to a list for printing without keeping autograd state alive.

76u_true: the exact closed form

For H=μ=1 and dpdx=-2 this simplifies to u(y) = y(1-y). The maximum is 0.25 at y=0.5. If your PINN prints 0.249 there, you have successfully taught a neural network to obey the Navier-Stokes equations.

41 lines without explanation

1import torch
2import torch.nn as nn
3import torch.optim as optim
4
5
6class ChannelFlowNN(nn.Module):
7    """
8    Tiny MLP that learns u(y) for plane Poiseuille flow.
9
10    We never give it the closed-form solution. Instead, we teach it to:
11      (a) satisfy the PDE     0 = -dp/dx + mu * d^2 u / dy^2     inside (0, H)
12      (b) satisfy the no-slip BCs   u(0) = 0,  u(H) = 0
13    """
14
15    def __init__(self, hidden=32):
16        super().__init__()
17        self.net = nn.Sequential(
18            nn.Linear(1, hidden), nn.Tanh(),
19            nn.Linear(hidden, hidden), nn.Tanh(),
20            nn.Linear(hidden, 1),
21        )
22
23    def forward(self, y):
24        return self.net(y)
25
26
27def pde_residual(model, y, mu, dpdx):
28    """
29    Compute  r(y) = -dp/dx + mu * d^2 u / dy^2  using autograd.
30    A perfect solution drives this to zero everywhere.
31    """
32    y = y.clone().requires_grad_(True)
33    u = model(y)
34
35    du_dy   = torch.autograd.grad(u.sum(), y, create_graph=True)[0]
36    d2u_dy2 = torch.autograd.grad(du_dy.sum(), y, create_graph=True)[0]
37
38    return -dpdx + mu * d2u_dy2
39
40
41# ---- Setup ----
42torch.manual_seed(0)
43H, mu, dpdx = 1.0, 1.0, -2.0
44
45model     = ChannelFlowNN(hidden=32)
46optimizer = optim.Adam(model.parameters(), lr=3e-3)
47
48# Collocation points in the interior of the channel.
49y_inner = torch.linspace(1e-3, H - 1e-3, 64).unsqueeze(1)
50
51# The two boundary points (no-slip walls).
52y_bc    = torch.tensor([[0.0], [H]])
53u_bc    = torch.tensor([[0.0], [0.0]])     # u(0) = u(H) = 0
54
55for step in range(2000):
56    optimizer.zero_grad()
57
58    # 1. PDE residual loss inside the channel
59    r = pde_residual(model, y_inner, mu, dpdx)
60    loss_pde = (r ** 2).mean()
61
62    # 2. Boundary-condition loss at the walls
63    loss_bc  = ((model(y_bc) - u_bc) ** 2).mean()
64
65    # 3. Combined loss
66    loss = loss_pde + 10.0 * loss_bc
67    loss.backward()
68    optimizer.step()
69
70    if step % 400 == 0:
71        print(f"step {step:4d}   loss_pde={loss_pde.item():.4e}"
72              f"   loss_bc={loss_bc.item():.4e}")
73
74# ---- Compare against the analytic solution ----
75y_test = torch.linspace(0.0, H, 11).unsqueeze(1)
76u_nn   = model(y_test).detach().squeeze()
77u_true = (1.0 / (2.0 * mu)) * dpdx * (y_test.squeeze() ** 2 - H * y_test.squeeze())
78print("y    :", y_test.squeeze().tolist())
79print("u_NN :", [round(v.item(), 4) for v in u_nn])
80print("u_ex :", [round(v.item(), 4) for v in u_true])

Why this works without training data

We never show the network a single example of the correct $u(y)$ . We only tell it (a) the PDE the solution must satisfy, and (b) the boundary values it must take. The PINN finds the unique function consistent with both, which is exactly the solution. This is the same paradigm used for high-dimensional PDEs where finite-difference grids would be impossibly large.

When PINNs struggle

PINNs converge gracefully for laminar problems like this one. For turbulent Navier-Stokes the loss landscape becomes punishing — BC enforcement and high-frequency PDE residuals interact in ways that plain Adam handles poorly. Research on PINNs is largely about rebalancing these competing losses.

Free Surfaces, Inflow, and Periodic BCs

Free-surface BCs (air-water interfaces)

When the “wall” is actually another fluid — like the air above an ocean wave — you cannot impose no-slip because nothing is anchored. Instead you split the BC into two parts:

Kinematic condition: fluid particles on the surface stay on the surface. If the surface is $y = \eta(x, t)$ , then $D\eta/Dt = v$ . The vertical velocity of the fluid equals the rate at which the surface itself moves up.
Dynamic condition: the stress at the interface is continuous. For a free surface with surface tension $\sigma$ and mean curvature $\kappa$ : $p_{\text{water}} = p_{\text{air}} + \sigma\,\kappa$ . This is why capillary waves exist and why small droplets are nearly spherical.

Inflow and outflow

For an open computational domain (e.g. a wind tunnel section), the upstream face usually gets a Dirichlet condition $u = u_{\text{in}}(y, t)$ while the downstream face gets a Neumann or “do-nothing” outflow condition $\partial u/\partial n = 0$ . Getting these right is one of the hardest parts of practical CFD.

Periodic BCs (the simulator's favourite)

Direct numerical simulations of homogeneous turbulence usually wrap the domain into a torus: $u(0, y, z) = u(L, y, z)$ for all $y, z$ . This eliminates spurious inflow/outflow effects and is essential for Fourier-spectral methods.

Why This Matters: From Pipelines to Capillaries

🛢️ Pipeline engineering

Oil and water pipelines are designed using exactly this closed form generalised to circular cross-sections (Hagen-Poiseuille). The pressure-driven parabola sets the pumping power required for a given throughput.

⚙️ Journal bearings

A spinning shaft in a sleeve creates Couette flow in the lubricant film. The pressure that pushes back on the shaft — and so the load capacity of the bearing — is calculated from the Reynolds equation, which is derived from exactly these BCs.

❤️ Microcirculation

Blood in capillaries flows at Reynolds number far below 1. Poiseuille's law (the pipe version of what we derived) sets the resistance of every vessel and is the basis of the entire field of haemodynamics.

🧪 Microfluidics

Lab-on-a-chip devices route nanolitres of fluid through channels a few microns wide. At those scales no-slip is the dominant physics — getting it wrong means designing a chip that simply does not flow.

✈️ Boundary-layer aerodynamics

The no-slip condition is the entire reason boundary layers exist. Prandtl's boundary-layer theory (1904) is built on top of this single BC and underlies every aerofoil ever designed.

🌊 Ocean & atmosphere

The air-sea interface is a free surface with kinematic and dynamic BCs. Wind transfers momentum into the ocean through wall shear at the surface — the engine of every current and storm surge.

Summary

Boundary conditions are not decoration on the Navier-Stokes equations — they are half of the problem. They are how a universal PDE becomes a specific flow.

Key Takeaways

The Navier-Stokes momentum equation is second-order in space, so every spatial direction needs two BCs per velocity component.
The no-slip condition $\mathbf{v}_{\text{fluid}} = \mathbf{v}_{\text{wall}}$ is the keystone BC for any solid surface. It is a Dirichlet condition and follows from molecular adhesion.
BCs come in flavours: Dirichlet (specify the value), Neumann (specify the derivative), free surface, periodic, and symmetry. Each models a different physical interface.
For steady 1-D flow between parallel plates the momentum equation collapses to an ODE whose solution is the Couette-Poiseuille profile $u(y) = \tfrac{1}{2\mu}\tfrac{dp}{dx}(y^2-Hy) + U_{\text{top}}\tfrac{y}{H} + U_{\text{bot}}(1-\tfrac{y}{H})$ .
The total profile is the superposition of a pressure-driven parabola and a wall-driven straight line — a direct consequence of the linearity of the simplified equation.
A small PINN built around autograd can learn the velocity profile purely from the PDE and BCs, with zero training data. The whole training loop is just “minimise PDE residual inside, minimise BC error at the walls.”

Boundary Conditions in One Sentence:

“Navier-Stokes is the verb; the boundary conditions are the subject and object — without them, no sentence.”

Coming Next: In the next section we move from these clean laminar profiles to the wild world of turbulent flow, where the same equations and BCs produce eddies, cascades, and chaos.