{"id":399350,"date":"2026-03-23T07:42:08","date_gmt":"2026-03-23T07:42:08","guid":{"rendered":"https:\/\/www.europesays.com\/ie\/399350\/"},"modified":"2026-03-23T07:42:08","modified_gmt":"2026-03-23T07:42:08","slug":"building-a-navier-stokes-solver-in-python-from-scratch-simulating-airflow","status":"publish","type":"post","link":"https:\/\/www.europesays.com\/ie\/399350\/","title":{"rendered":"Building a Navier-Stokes Solver in Python from Scratch: Simulating Airflow"},"content":{"rendered":"<p class=\"wp-block-paragraph\"> (CFD) is often seen as a black box of complex commercial software. However, implementing a solver \u201cfrom scratch\u201d is one of the most powerful ways to learn the physics of fluid motion. I started this as a personal project, and as a part of a course on Biophysics, I took it as an opportunity to finally understand how these beautiful simulations work.<\/p>\n<p class=\"wp-block-paragraph\">This guide is designed for data scientists and engineers who want to move beyond high-level libraries and understand the underlying mechanics of numerical simulations by translating partial differential equations into discretized Python code. We will also explore fundamental programming concepts like vectorized operations with NumPy and stochastic convergence, which are essential skills for everyone interested in broader scientific computing and machine learning architectures.<\/p>\n<p class=\"wp-block-paragraph\">We will walk through the derivation and Python implementation of a simple incompressible Navier-Stokes (NS) solver. And then, we will then apply this solver to simulate airflow around a bird\u2019s wing profile.<\/p>\n<p>The Physics: Incompressible Navier-Stokes<\/p>\n<p class=\"wp-block-paragraph\">The fundamental equations of CFD are the Navier-Stokes equations, which describe how velocity and pressure evolve in a fluid. For steady flight (like a bird gliding), we assume the air is incompressible (constant density) and laminar. They can be understood as just Newton\u2019s motion law, but for an infinitesimal element of fluid, with the forces that affect it. That\u2019s mostly pressure and viscosity, but depending on the context you could add in gravity, mechanical stresses, or even electromagnetism if you\u2019re feeling like it. The author can attest to this being very much not recommend for a first project.\u00a0<\/p>\n<p class=\"wp-block-paragraph\">The equations in vector form are:<\/p>\n<p class=\"wp-block-shortcode\">\\[<br \/>\\frac{\\partial \\mathbf{v}}{\\partial t} + (\\mathbf{v} \\cdot \\nabla)\\mathbf{v} = -\\frac{1}{\\rho}\\nabla p + \\nu \\nabla^2 \\mathbf{v} \\\\<br \/>\\nabla \\cdot \\mathbf{v} = 0<br \/>\\]<\/p>\n<p class=\"wp-block-paragraph\">Where:<\/p>\n<ul class=\"wp-block-list\">\n<li class=\"wp-block-list-item\"><strong>v<\/strong>: Velocity field (u,v)<\/li>\n<li class=\"wp-block-list-item\"><strong>p<\/strong>: Pressure<\/li>\n<li class=\"wp-block-list-item\"><strong>\u03c1<\/strong>: Fluid density<\/li>\n<li class=\"wp-block-list-item\"><strong>\u03bd<\/strong>: Kinematic viscosity<\/li>\n<\/ul>\n<p class=\"wp-block-paragraph\">The first equation (Momentum) balances inertia against pressure gradients and viscous diffusion. The second equation (Continuity) enforces that the fluid density remains constant.<\/p>\n<p>The Pressure Coupling\u00a0Problem<\/p>\n<p class=\"wp-block-paragraph\">A major challenge in CFD is that pressure and velocity are coupled: the pressure field must adjust constantly to ensure the fluid remains incompressible.<\/p>\n<p class=\"wp-block-paragraph\">To solve this, we derive a <strong>Pressure-Poisson equation<\/strong> by taking the divergence of the momentum equation. In a discretized solver, we solve this Poisson equation at every single timestep to update the pressure, ensuring the velocity field remains divergence-free.<\/p>\n<p>Discretization: From Math to\u00a0Grid<\/p>\n<p class=\"wp-block-paragraph\">To solve these equations on a computer, we use <strong>Finite Difference<\/strong> schemes on a uniform grid.<\/p>\n<ul class=\"wp-block-list\">\n<li class=\"wp-block-list-item\"><strong>Time:<\/strong> Forward difference (Explicit Euler).<\/li>\n<li class=\"wp-block-list-item\"><strong>Advection (Nonlinear terms):<\/strong> Backward\/Upwind difference (for stability).<\/li>\n<li class=\"wp-block-list-item\"><strong>Diffusion &amp; Pressure:<\/strong> Central difference.<\/li>\n<\/ul>\n<p class=\"wp-block-paragraph\">For example, the update formula for the u (x-velocity) component looks like this in finite-difference form<\/p>\n<p class=\"wp-block-shortcode\">\\[<br \/>\\begin{aligned}<br \/>u_{i,j}^{n+1} = \\; &amp; u_{i,j}^n<br \/>\u2013 u_{i,j}^n \\frac{\\Delta t}{\\Delta x}\\left(u_{i,j}^n \u2013 u_{i-1,j}^n\\right)<br \/>\u2013 v_{i,j}^n \\frac{\\Delta t}{\\Delta y}\\left(u_{i,j}^n \u2013 u_{i,j-1}^n\\right) \\\\<br \/>&amp; \u2013 \\frac{\\Delta t}{2\\rho \\Delta x}\\left(p_{i+1,j}^n \u2013 p_{i-1,j}^n\\right)<br \/>\\\\ &amp; + \\nu \\left(\\frac{\\Delta t}{\\Delta x^2}\\left(u_{i+1,j}^n \u2013 2u_{i,j}^n + u_{i-1,j}^n\\right)<br \/>+ \\frac{\\Delta t}{\\Delta y^2}\\left(u_{i,j+1}^n \u2013 2u_{i,j}^n + u_{i,j-1}^n\\right)\\right)<br \/>\\end{aligned}<br \/>\\]<\/p>\n<p class=\"wp-block-paragraph\">In code, the advection term u\u2202x\u2202u\u200b uses a backward difference:\u00a0<\/p>\n<p class=\"wp-block-shortcode\">\\[<br \/>u_{i,j}^n \\frac{u_{i,j}^n\u200a-\u200au_{i-1,j}^n}{\\Delta x}<br \/>\\]<\/p>\n<p><img decoding=\"async\" src=\"https:\/\/www.europesays.com\/ie\/wp-content\/uploads\/2026\/03\/image-276-1024x370.png\" alt=\"\" class=\"wp-image-651268\"\/><\/p>\n<p>The Python Implementation<\/p>\n<p class=\"wp-block-paragraph\">The implementation proceeds in four distinct steps using NumPy arrays.<\/p>\n<p>1. Initialization<\/p>\n<p class=\"wp-block-paragraph\">We define the grid size (nx, ny), time step (dt), and physical parameters (rho, nu). We initialize velocity fields (u,v) and pressure (p) to zeros or a uniform flow.<\/p>\n<p>2. The Wing Geometry (Immersed Boundary)<\/p>\n<p class=\"wp-block-paragraph\">To simulate a wing on a Cartesian grid, we need to mark which grid points lie inside the solid wing.<\/p>\n<ul class=\"wp-block-list\">\n<li class=\"wp-block-list-item\">We load a wing mesh (e.g., from an STL file).<\/li>\n<li class=\"wp-block-list-item\">We create a Boolean mask array where True indicates a point inside the wing.<\/li>\n<li class=\"wp-block-list-item\">During the simulation, we force velocity to zero at these masked points (no-slip\/no-penetration condition).<\/li>\n<\/ul>\n<p>3. The Main Solver\u00a0Loop<\/p>\n<p class=\"wp-block-paragraph\">The core loop repeats until the solution reaches a steady state. The steps are:<\/p>\n<ol class=\"wp-block-list\">\n<li class=\"wp-block-list-item\"><strong>Build the Source Term (b):<\/strong> Calculate the divergence of the velocity terms.<\/li>\n<li class=\"wp-block-list-item\"><strong>Solve Pressure:<\/strong> Solve the Poisson equation for p using Jacobi iteration.<\/li>\n<li class=\"wp-block-list-item\"><strong>Update Velocity:<\/strong> Use the new pressure to update u and v.<\/li>\n<li class=\"wp-block-list-item\"><strong>Apply Boundary Conditions:<\/strong> Enforce inlet velocity and zero velocities inside the wing.<\/li>\n<\/ol>\n<p>The Code<\/p>\n<p class=\"wp-block-paragraph\">Here is how the core mathematical updates look in Python (vectorized for performance).<\/p>\n<p class=\"wp-block-paragraph\"><strong>Step A: Building the Pressure Source Term<\/strong> This represents the Right-Hand Side (RHS) of the Poisson equation based on current velocities.<\/p>\n<p># b is the source term<br \/>\n# u and v are current velocity arrays<br \/>\nb[1:-1, 1:-1] = (rho * (<br \/>\n    1 \/ dt * ((u[1:-1, 2:] &#8211; u[1:-1, 0:-2]) \/ (2 * dx) +<br \/>\n              (v[2:, 1:-1] &#8211; v[0:-2, 1:-1]) \/ (2 * dy)) &#8211;<br \/>\n    ((u[1:-1, 2:] &#8211; u[1:-1, 0:-2]) \/ (2 * dx))**2 &#8211;<br \/>\n    2 * ((u[2:, 1:-1] &#8211; u[0:-2, 1:-1]) \/ (2 * dy) *<br \/>\n         (v[1:-1, 2:] &#8211; v[1:-1, 0:-2]) \/ (2 * dx)) &#8211;<br \/>\n    ((v[2:, 1:-1] &#8211; v[0:-2, 1:-1]) \/ (2 * dy))**2<br \/>\n))<\/p>\n<p class=\"wp-block-paragraph\"><strong>Step B: Solving for Pressure (Jacobi Iteration)<\/strong> We iterate to smooth out the pressure field until it balances the source term.<\/p>\n<p>for _ in range(nit):<br \/>\n    pn = p.copy()<br \/>\n    p[1:-1, 1:-1] = (<br \/>\n        (pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +<br \/>\n        (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2 &#8211;<br \/>\n        b[1:-1, 1:-1] * dx**2 * dy**2<br \/>\n    ) \/ (2 * (dx**2 + dy**2))<br \/>\n# Boundary conditions: p=0 at edges (gauge pressure)<br \/>\n    p[:, -1] = 0; p[:, 0] = 0; p[-1, :] = 0; p[0, :] = 0<\/p>\n<p class=\"wp-block-paragraph\"><strong>Step C: Updating Velocity<\/strong> Finally, we update the velocity using the explicit discretized momentum equations.<\/p>\n<p>un = u.copy()<br \/>\nvn = v.copy()<br \/>\n# Update u (x-velocity)<br \/>\nu[1:-1, 1:-1] = (un[1:-1, 1:-1] &#8211;<br \/>\n                 un[1:-1, 1:-1] * dt \/ dx * (un[1:-1, 1:-1] &#8211; un[1:-1, 0:-2]) &#8211;<br \/>\n                 vn[1:-1, 1:-1] * dt \/ dy * (un[1:-1, 1:-1] &#8211; un[0:-2, 1:-1]) &#8211;<br \/>\n                 dt \/ (2 * rho * dx) * (p[1:-1, 2:] &#8211; p[1:-1, 0:-2]) +<br \/>\n                 nu * (dt \/ dx**2 * (un[1:-1, 2:] &#8211; 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +<br \/>\n                       dt \/ dy**2 * (un[2:, 1:-1] &#8211; 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])))<br \/>\n# Update v (y-velocity)<br \/>\nv[1:-1, 1:-1] = (vn[1:-1, 1:-1] &#8211;<br \/>\n                 un[1:-1, 1:-1] * dt \/ dx * (vn[1:-1, 1:-1] &#8211; vn[1:-1, 0:-2]) &#8211;<br \/>\n                 vn[1:-1, 1:-1] * dt \/ dy * (vn[1:-1, 1:-1] &#8211; vn[0:-2, 1:-1]) &#8211;<br \/>\n                 dt \/ (2 * rho * dy) * (p[2:, 1:-1] &#8211; p[0:-2, 1:-1]) +<br \/>\n                 nu * (dt \/ dx**2 * (vn[1:-1, 2:] &#8211; 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +<br \/>\n                       dt \/ dy**2 * (vn[2:, 1:-1] &#8211; 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))<\/p>\n<p>Results: Does it\u00a0Fly?<\/p>\n<p class=\"wp-block-paragraph\">We ran this solver on a rigid wing profile with a constant far-field inflow.<\/p>\n<p><img decoding=\"async\" src=\"https:\/\/www.europesays.com\/ie\/wp-content\/uploads\/2026\/03\/Wing_speed2D.png\" alt=\"\" class=\"wp-image-650645\" style=\"border-style:none;border-width:0px;border-top-left-radius:0px;border-top-right-radius:0px;border-bottom-left-radius:0px;border-bottom-right-radius:0px\"\/><\/p>\n<p class=\"wp-block-paragraph\"><strong>Qualitative Observations<\/strong> The results align with physical expectations. The simulations show high pressure beneath the wing and low pressure above it, which is exactly the mechanism that generates lift. Velocity vectors show the airflow accelerating over the top surface (Bernoulli\u2019s principle).<\/p>\n<p class=\"wp-block-paragraph\"><strong>Forces: Lift vs. Drag<\/strong> By integrating the pressure field over the wing surface, we can calculate lift.<\/p>\n<ul class=\"wp-block-list\">\n<li class=\"wp-block-list-item\">The solver demonstrates that <strong>pressure forces dominate<\/strong> viscous friction forces by a factor of nearly 1000x in air.<\/li>\n<li class=\"wp-block-list-item\">As the angle of attack increases (from 0\u2218 to \u221220\u2218), the lift-to-drag ratio rises, matching trends seen in wind tunnels and professional CFD packages like OpenFOAM.<\/li>\n<\/ul>\n<p><img decoding=\"async\" src=\"https:\/\/www.europesays.com\/ie\/wp-content\/uploads\/2026\/03\/LiftDrag1.png\" alt=\"\" class=\"wp-image-650646\"\/><\/p>\n<p>Limitations &amp; Next\u00a0Steps<\/p>\n<p class=\"wp-block-paragraph\">While making this solver was great for learning, the tool itself has its limitations:<\/p>\n<ul class=\"wp-block-list\">\n<li class=\"wp-block-list-item\"><strong>Resolution:<\/strong> 3D simulations on a Cartesian grid are computationally expensive and require coarse grids, making quantitative results less reliable.<\/li>\n<li class=\"wp-block-list-item\"><strong>Turbulence:<\/strong> The solver is laminar; it lacks a turbulence model (like k\u2212\u03f5) required for high-speed or complex flows.<\/li>\n<li class=\"wp-block-list-item\"><strong>Diffusion:<\/strong> Upwind differencing schemes are stable but numerically diffusive, potentially \u201csmearing\u201d out fine flow details.<\/li>\n<\/ul>\n<p class=\"wp-block-paragraph\"><strong>Where to go from here?<\/strong> This project serves as a starting point. Future improvements could include implementing higher-order advection schemes (like WENO), adding turbulence modeling, or moving to Finite Volume methods (like OpenFOAM) for better mesh handling around complex geometries. There are lots of clever techniques to get around the plethora of scenarios that you may want to implement. This is just a first step on the direction of really understanding CFD!<\/p>\n","protected":false},"excerpt":{"rendered":"(CFD) is often seen as a black box of complex commercial software. However, implementing a solver \u201cfrom scratch\u201d&hellip;\n","protected":false},"author":2,"featured_media":399351,"comment_status":"","ping_status":"","sticky":false,"template":"","format":"standard","meta":{"footnotes":"","_share_on_mastodon":"0"},"categories":[271],"tags":[52137,53041,18,19,17,452,59808,133,35058],"class_list":["post-399350","post","type-post","status-publish","format-standard","has-post-thumbnail","category-physics","tag-aerodynamics","tag-cfd","tag-eire","tag-ie","tag-ireland","tag-physics","tag-python","tag-science","tag-simulation"],"share_on_mastodon":{"url":"https:\/\/pubeurope.com\/@ie\/116277369052530335","error":""},"_links":{"self":[{"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/posts\/399350","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/comments?post=399350"}],"version-history":[{"count":0,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/posts\/399350\/revisions"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/media\/399351"}],"wp:attachment":[{"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/media?parent=399350"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/categories?post=399350"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/tags?post=399350"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}