Damage Modeling using Convection Diffusion Equation

     Inspired from fluid dynamics, where the convection diffusion equation is used in several areas including to model turbulence and shocks (discontinuities), a damage transport model has been encoded in the code based on the convection diffusion equation. The code has been modified to create material damage using the strain energy density. Material damage, d, is then added into the damage transport equation as a source term. The modified stress equations takes the form

σ_d = (1 - d) * σ

     where σ_d is stress with material damage included. The damage transport equation is of the form

d_t + u * d_x + v * d_y = dc * (d_xx + d_yy) + S

     d_t is time derivative of damage, d_x and d_y are spatial derivatives of damage. d_xx and d_yy diffuse the damage in the material with coefficient dc, and S is the strain energy density. The u and v velocities in the damage transport equation are calculated using time rate of change of x and y displacements. The damage variable is shown in Fig. 1, along with resultant displacement superimposed on the undeformed beam. As the discrete formulation of the problem is explicit, the simulation ran on a very coarse mesh to keep the time step large. A failure mask is used to mark high stress regions both in compression and tension. The first and third principal stresses are used as the criteria for fracture propagation.

Fig. 1, Displacement and damage progression

Code

#Copyright <2025> <FAHAD BUTT>

#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

#%% import necessary libraries

import numpy as np

import matplotlib.pyplot as plt

#%% define parameters

L = 0.2 # domain length

D = 0.04 # domain width

h = D / 10 # grid spacing

h2 = 2 * h # multiplication factor

Nx = round(L / h) # grid points in x-axis

Ny = round(D / h) # grid points in y-axis

nu = 0.175 # Poisson's ratio

E = 3e10 # Young's modulus [Pa]

alpha = E / ((1 + nu) * (1 - (2 * nu))) # multiplication parameter

beta = E / (2 * (1 + nu)) # multiplication parameter

rho = 2350 # density [Kg/m3]

rho_2_h = rho * h2 # multiplication factor

cp = np.sqrt(E * (1 - nu) / (rho * (1 + nu) * (1 - 2 * nu))) # wave speed

CFL = 0.9 # CFL number

dt = CFL * h / cp # time step

stress_rate = 1e5 # Pa/s

compressive_strength = -35e6 # maximum compressive strength [MPa]

tensile_strength = -0.1 * compressive_strength # maximum tensile strength [MPa]

s_initial = 0 # initial load [MPa]

travel = 3000000 # times the disturbance travels entire length of computational domain

TT = travel * L / cp # total time

ns = int(TT / dt) # number of time steps

dc = 1e-5 # damage diffusion

print(f"Time step [s]: {dt:.3}") # show time step

print(f"Total time-steps: {ns:.0f}") # show total time steps

#%% initialize fields

u = np.zeros((Nx, Ny)) # u displacement (next time step) u_n+1

v = np.zeros((Nx, Ny)) # v displacement (next time step)

s_xx = np.zeros((Nx, Ny)) # normal stress x

s_yy = np.zeros((Nx, Ny)) # normal stress y

s_xy = np.zeros((Nx, Ny)) # shear stress xy

e_xx = np.zeros((Nx, Ny)) # normal strain x

e_yy = np.zeros((Nx, Ny)) # normal strain y

e_xy = np.zeros((Nx, Ny)) # shear strain xy

sigma_1 = np.zeros((Nx, Ny)) # first principal stress

sigma_3 = np.zeros((Nx, Ny)) # third principal stress

u_vel = np.zeros((Nx, Ny)) # initial velocity in x-direction

v_vel = np.zeros((Nx, Ny)) # initial velocity in y-direction

un = np.copy(u) # u displacement (current time step) u

unn = un - dt * u_vel # u displacement (previous time step) u_n-1

vn = np.copy(v) # v displacement (current time step) v

vnn = vn - dt * v_vel # v displacement (previous time step v_n-1

d = np.zeros((Nx, Ny)) # material damage

H = np.zeros((Nx, Ny)) # strain energy density

#%% solve equations of motion

for nt in range(ns):

    e_xx[1:-1, 1:-1] = (u[2:, 1:-1] - u[:-2, 1:-1]) / h2 # normal strain x

    # boundary conditions for e_xx

    e_xx[0, :] = e_xx[1, :] # de_xx/dx = 0 at x = 0

    e_xx[-1, :] = e_xx[-2, :] # de_xx/dx = 0 at x = L

    e_xx[:, 0] = e_xx[:, 1] # de_xx/dy = 0 at y = 0

    e_xx[:, -1] = e_xx[:, -2] # de_xx/dy = 0 at y = D

    

    e_yy[1:-1, 1:-1] = (v[1:-1, 2:] - v[1:-1, :-2]) / h2 # normal strain y

    # boundary conditions for e_yy

    e_yy[0, :] = e_yy[1, :] # de_yy/dx = 0 at x = 0

    e_yy[-1, :] = e_yy[-2, :] # de_yy/dx = 0 at x = L

    e_yy[:, 0] = e_yy[:, 1] # de_yy/dy = 0 at y = 0

    e_yy[:, -1] = e_yy[:, -2] # de_yy/dy = 0 at y = D

    

    e_xy[1:-1, 1:-1] = 0.5 * (u[1:-1, 2:] - u[1:-1, :-2] + v[2:, 1:-1] - v[:-2, 1:-1]) / h2 # shear strain xy

    # boundary conditions for e_xy

    e_xy[0, :] = e_xy[1, :] # de_xy/dx = 0 at x = 0

    e_xy[-1, :] = e_xy[-2, :] # de_xy/dx = 0 at x = L

    e_xy[:, 0] = e_xy[:, 1] # de_xy/dy = 0 at y = 0

    e_xy[:, -1] = e_xy[:, -2] # de_xy/dy = 0 at y = D

    

    s_xx[1:-1, 1:-1] = (1 - d[1:-1, 1:-1]) * alpha * ((1 - nu) * e_xx[1:-1, 1:-1] + nu * e_yy[1:-1, 1:-1]) # normal stress x

    # boundary conditions for s_xy

    s_xx[0, :] = 0 # s_xx = 0 at x = 0

    s_xx[-1, :] = 0 # s_xx = 0 at x = L

    s_xx[:, 0] = s_xx[:, 1] # ds_xx/dy = 0 at y = 0

    s_xx[:, -1] = s_xx[:, -2] # ds_xx/dy = 0 at y = D

    s_yy[1:-1, 1:-1] = (1 - d[1:-1, 1:-1]) * alpha * ((1 - nu) * e_yy[1:-1, 1:-1] + nu * e_xx[1:-1, 1:-1]) # normal stress y

    # boundary conditions for s_yy

    s_yy[0, :] = s_yy[1, :] # ds_yy/dx = 0 at x = 0

    s_yy[-1, :] = s_yy[-2, :] # ds_yy/dx = 0 at x = L

    s_yy[:, 0] = 0 # s_yy = 0 at y = 0

    s_yy[round((L * 0.1) * Nx / L) - 1:round((L * 0.1) * Nx / L) + 1, 0] = s_yy[round((L * 0.1) * Nx / L) - 1:round((L * 0.1) * Nx / L) + 1, 1] # ds_yy/dy = 0 at x = 0.1 and y = 0

    s_yy[round((L * 0.9) * Nx / L) - 1:round((L * 0.9) * Nx / L) + 1, 0] = s_yy[round((L * 0.9) * Nx / L) - 1:round((L * 0.9) * Nx / L) + 1, 1] # ds_yy/dy = 0 at x = 0.9 and y = 0

    s_yy[:, -1] = 0 # s_yy = 0 at y = D

    s_0 = s_initial + stress_rate * nt * dt # applied load

    s_yy[round((L * 0.5) * Nx / L) - 1:round((L * 0.5) * Nx / L) + 1, -1] = -s_0 # s_yy = s_0 at x = 0.5 and y = D

    

    s_xy[1:-1, 1:-1] = (1 - d[1:-1, 1:-1]) * beta * e_xy[1:-1, 1:-1] # shear stress xy

    # boundary conditions for s_xy

    s_xy[0, :] = 0 # s_xy = 0 at x = 0

    s_xy[-1, :] = 0 # s_xy = 0 at x = L

    s_xy[:, 0] = 0 # s_xy = 0 at y = 0

    s_xy[round((L * 0.1) * Nx / L) - 1:round((L * 0.1) * Nx / L) + 1, 0] = s_xy[round((L * 0.1) * Nx / L) - 1:round((L * 0.1) * Nx / L) + 1, 1] # ds_xy/dy = 0 at x = 0.1 and y = 0

    s_xy[round((L * 0.9) * Nx / L) - 1:round((L * 0.9) * Nx / L) + 1, 0] = s_xy[round((L * 0.9) * Nx / L) - 1:round((L * 0.9) * Nx / L) + 1, 1] # ds_xy/dy = 0 at x = 0.9 and y = 0

    s_xy[:, -1] = 0 # s_xy = 0 at y = D

    s_xy[round((L * 0.5) * Nx / L) - 1:round((L * 0.5) * Nx / L) + 1, -1] = s_xy[round((L * 0.5) * Nx / L) - 1:round((L * 0.5) * Nx / L) + 1, -2] # ds_xy/dy = 0 at y = D

    u[1:-1, 1:-1] = 2 * un[1:-1, 1:-1] - unn[1:-1, 1:-1] + (dt**2 / rho_2_h) * (s_xx[2:, 1:-1] - s_xx[:-2, 1:-1] + s_xy[1:-1, 2:] - s_xy[1:-1, :-2]) # x momentum

    # boundary conditions for u

    u[0, :] = u[1, :] # du/dx = 0 at x = 0

    u[-1, :] = u[-2, :] # du/dx = 0 at x = L

    u[:, 0] = u[:, 1] # du/dy = 0 at y = 0

    u[:, -1] = u[:, -2] # du/dy = 0 at y = D

    un[:, :] = u[:, :] # shift u_n+1 to u_n

    unn[:, :] = un[:, :] # shift u_n to u_n-1

    

    v[1:-1, 1:-1] = 2 * vn[1:-1, 1:-1] - vnn[1:-1, 1:-1] + (dt**2 / rho_2_h) * (s_xy[2:, 1:-1] - s_xy[:-2, 1:-1] + s_yy[1:-1, 2:] - s_yy[1:-1, :-2]) # y momentum

    # boundary conditions for v

    v[0, :] = v[1, :] # dv/dx = 0 at x = 0

    v[-1, :] = v[-2, :] # dv/dx = 0 at x = L

    v[:, 0] = v[:, 1] # dv/dy = 0 at y = 0

    v[:, -1] = v[:, -2] # dv/dy = 0 at y = D

    v[round((L * 0.1) * Nx / L) - 1:round((L * 0.1) * Nx / L) + 1, 0] = 0 # v = 0 at x = 0.1 and y = 0

    v[round((L * 0.9) * Nx / L) - 1:round((L * 0.9) * Nx / L) + 1, 0] = 0 # v = 0 at x = 0.9 and y = 0

    vn[:, :] = v[:, :] # shift u_n+1 to u_n

    vnn[:, :] = vn[:, :] # shift u^n to u_n-1

    sigma_1[1:-1, 1:-1] = 0.5 * (s_xx[1:-1, 1:-1] + s_yy[1:-1, 1:-1]) + np.sqrt((0.5 * (s_xx[1:-1, 1:-1] - s_yy[1:-1, 1:-1]))**2 + s_xy[1:-1, 1:-1]**2)  # first principal stress

    # boundary conditions for sigma_1

    sigma_1[0, :] = sigma_1[1, :] # dsigma_1/dx = 0 at x = 0

    sigma_1[-1, :] = sigma_1[-2, :] # dsigma_1/dx = 0 at x = L

    sigma_1[:, 0] = sigma_1[:, 1] # dsigma_1/dy = 0 at y = 0

    sigma_1[:, -1] = sigma_1[:, -2] # dsigma_1/dy = 0 at y = D

    

    sigma_3[1:-1, 1:-1] = 0.5 * (s_xx[1:-1, 1:-1] + s_yy[1:-1, 1:-1]) - np.sqrt((0.5 * (s_xx[1:-1, 1:-1] - s_yy[1:-1, 1:-1]))**2 + s_xy[1:-1, 1:-1]**2)  # third principal stress

    # boundary conditions for sigma_3

    sigma_3[0, :] = sigma_3[1, :] # dsigma_3/dx = 0 at x = 0

    sigma_3[-1, :] = sigma_3[-2, :] # dsigma_3/dx = 0 at x = L

    sigma_3[:, 0] = sigma_3[:, 1] # dsigma_3/dy = 0 at y = 0

    sigma_3[:, -1] = sigma_3[:, -2] # dsigma_3/dy = 0 at y = D

    H[1:-1, 1:-1] = 0.5 * (s_xx[1:-1, 1:-1] * e_xx[1:-1, 1:-1] + s_yy[1:-1, 1:-1] * e_yy[1:-1, 1:-1] + 2 * s_xy[1:-1, 1:-1] * e_xy[1:-1, 1:-1]) # energy density (damage creation)

    # boundary conditions for H

    H[0, :] = H[1, :] # dH/dx = 0 at x = 0

    H[-1, :] = H[-2, :] # dH/dx = 0 at x = L

    H[:, 0] = H[:, 1] # dH/dy = 0 at y = 0

    H[:, -1] = H[:, -2] # dH/dy = 0 at y = D

    

    failure_mask = (sigma_1 > tensile_strength) | (sigma_3 < compressive_strength) # mark high stress regions

    dn = d.copy() # shift d_n+1 to d_n

    d[1:-1, 1:-1] = np.where(failure_mask[1:-1, 1:-1], dn[1:-1, 1:-1] - (1 / h2) * ((u[1:-1, 1:-1] - un[1:-1, 1:-1]) * (dn[2:, 1:-1] - dn[:-2, 1:-1]) + (v[1:-1, 1:-1] - vn[1:-1, 1:-1]) * (dn[1:-1, 2:] - dn[1:-1, :-2])) + dt * dc / h**2 * (dn[2:, 1:-1] + dn[:-2, 1:-1] + dn[1:-1, 2:] + dn[1:-1, :-2] - 4 * dn[1:-1, 1:-1]) + dt * (1 - dn[1:-1, 1:-1])**2 * H[1:-1, 1:-1], d[1:-1, 1:-1]) # damage

    d[0, :] = d[1, :] # dd/dx = 0 at x = 0

    d[-1, :] = d[-2, :] # dd/dx = 0 at x = L

    d[:, 0] = d[:, 1] # dd/dy = 0 at y = 0

    d[:, -1] = d[:, -2] # dd/dy = 0 at y = D

    d = np.clip(d, 0, 1) # keep damage between 0 and 1

#%% plotting

x = np.linspace(0, L, Nx) # initialize x

y = np.linspace(0, D, Ny) # initialize y

X, Y = np.meshgrid(x, y) # create combinations of sample points along the domain

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, (u * 1000).T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('u')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, (v * 1000).T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('v')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, (s_xx * 1e-6).T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('s_xx')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, (s_yy * 1e-6).T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('s_yy')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, (s_xy * 1e-6).T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('s_xy')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, e_xx.T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('e_xx')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, e_yy.T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('e_yy')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, e_xy.T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('e_xy')

plt.xticks([])

plt.yticks([])

plt.show()

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, H.T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('H')

plt.xticks([])

plt.yticks([])

plt.show()

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, (sigma_1 * 1e-6).T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('sigma_1')

plt.xticks([])

plt.yticks([])

plt.show()

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, (sigma_3 * 1e-6).T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('sigma_3')

plt.xticks([])

plt.yticks([])

plt.show()

plt.show()

plt.figure(dpi = 200) # make a nice crisp image :)

plt.contourf(X, Y, d.T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('d')

plt.xticks([])

plt.yticks([])

plt.show()

plt.show()

     Thank you for reading. If you want to hire me as your shining new post-doc, do reach out!

Multi-Physics Simulations using in-House Code

     One fine morning, in abundant spare time, yours truly decided to adapt the CFD and LE code that has been developed over the years in to sort of a code capable of simulating one-way FSI. In this post details about Fluid Structure Interaction (FSI) code are made available.

     The code is very simple and basic form is less than 150 lines. A total of 50 lines for the CFD portion and around about 75 lines for the LE (Linear Elasticity) portion. The both CFD and LE codes have be thoroughly validated undependably. For CFD validation, please refer to here and here. Meanwhile the validation for LE portion is available here.

     As of now, the work process is CFD→LE. Which means the flow-field has to be calculated first. More details about the CFD portion are available here, here, and here. The pressure and shear stress on the walls of interest are extracted. This pressure and shear stress is then applied as stress loads on the a boundaries of the object in the LE code.

     The flow-field is simulated around a thin slander beam which is fixed with the floor. Within Fig. 1, From left to right u and v components of velocity are shown. Meanwhile, pressure and streamlines are shown on the bottom row. Within Fig. 2, the deflected thin beam is shown. The deflection is consistent with the intuition i.e. towards the flow and slightly below.

Fig. 1, The CFD results

Fig. 2, The results of LE simulation!

     Thank you for reading. If you want to hire me as your shining new post-doc, do reach out!

Linear Elasticity in 75 Lines!

      In abundant spare time ⏳, yours truly has implemented the linear elasticity 〰️ equations in a finite difference method code. The code 🖳 in it's simplest form is less than 75 lines including importing libraries and plotting! 😲 For validation, refer here. More validation examples will be shared in this blog post as these become available. The purpose of this code is not to run linear elasticity simulations but this code is supposed to be used for teaching and has another higher purpose, which will be made available shortly 😁. Happy codding❗❕

Code

#Copyright <2025> <FAHAD BUTT>

#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

import numpy as np

import matplotlib.pyplot as plt

L = 0.5 # domain length

D = 0.5 # domain width

h = D / 250 # grid spacing

h2 = 2 * h # multiplication factor

Nx = round(L / h) + 1 # grid points in x-axis

Ny = round(D / h) + 1 # grid points in y-axis

nu = 0.3 # Poisson's ratio

s_0 = 1e6 # applied load [Pa]

E = 2e11 # Young's modulus [Pa]

alpha = E / ((1 + nu) * (1 - (2 * nu))) # non dimensional parameter

beta = E / (2 * (1 + nu)) # non dimensional parameter

u = np.zeros((Nx, Ny)) # u displacement

v = np.zeros((Nx, Ny)) # v displacement

s_xx = np.zeros((Nx, Ny)) # normal stress x

s_yy = np.zeros((Nx, Ny)) # normal stress y

s_xy = np.zeros((Nx, Ny)) # shear stress xy

e_xx = np.zeros((Nx, Ny)) # normal strain x

e_yy = np.zeros((Nx, Ny)) # normal strain y

e_xy = np.zeros((Nx, Ny)) # shear strain xy

epsilon_d = 1e-17 # stability factor

for nt in range(100000): # solve equations of motion

    e_xx[1:-1, 1:-1] = (u[2:, 1:-1] - u[:-2, 1:-1]) / h2 # normal strain x

    e_xx[0, :] = e_xx[1, :] # de_xx/dx = 0 at x = 0

    e_xx[-1, :] = e_xx[-2, :] # de_xx/dx = 0 at x = L

    e_xx[:, 0] = e_xx[:, 1] # de_xx/dy = 0 at y = 0

    e_xx[:, -1] = e_xx[:, -2] # de_xx/dy = 0 at y = D    

    e_yy[1:-1, 1:-1] = (v[1:-1, 2:] - v[1:-1, :-2]) / h2 # normal strain y

    e_yy[0, :] = e_yy[1, :] # de_yy/dx = 0 at x = 0

    e_yy[-1, :] = e_yy[-2, :] # de_yy/dx = 0 at x = L

    e_yy[:, 0] = e_yy[:, 1] # de_yy/dy = 0 at y = 0

    e_yy[:, -1] = e_yy[:, -2] # de_yy/dy = 0 at y = D 

    e_xy[1:-1, 1:-1] = 0.5 * (u[1:-1, 2:] - u[1:-1, :-2] + v[2:, 1:-1] - v[:-2, 1:-1]) / h2 # shear strain xy

    e_xy[0, :] = e_xy[1, :] # de_xy/dx = 0 at x = 0

    e_xy[-1, :] = e_xy[-2, :] # de_xy/dx = 0 at x = L

    e_xy[:, 0] = e_xy[:, 1] # de_xy/dy = 0 at y = 0

    e_xy[:, -1] = e_xy[:, -2] # de_xy/dy = 0 at y = D

    s_xx[1:-1, 1:-1] = alpha * ((1 - nu) * e_xx[1:-1, 1:-1] + nu * e_yy[1:-1, 1:-1]) # normal stress x

    s_xx[0, :] = s_xx[1, :] # ds_xx/dx = 0 at x = 0

    s_xx[-1, :] = s_0 # s_xx = s_0 at x = L

    s_xx[:, 0] = s_xx[:, 1] # ds_xx/dy = 0 at y = 0

    s_xx[:, -1] = s_xx[:, -2] # ds_xx/dy = 0 at y = D

    s_yy[1:-1, 1:-1] = alpha * ((1 - nu) * e_yy[1:-1, 1:-1] + nu * e_xx[1:-1, 1:-1]) # normal stress y

    s_yy[0, :] = s_yy[1, :] # ds_yy/dx = 0 at x = 0

    s_yy[-1, :] = s_yy[-2, :] # ds_yy/dx = 0 at x = L

    s_yy[:, 0] = s_yy[:, 1] # ds_yy/dy = 0 at y = 0

    s_yy[:, -1] = 0 # s_yy = 0 at y = D

    s_xy[1:-1, 1:-1] = beta * e_xy[1:-1, 1:-1] # shear stress xy

    s_xy[0, :] = s_xy[1, :] # ds_xy/dx = 0 at x = 0

    s_xy[-1, :] = 0 # s_xy = 0 at x = L

    s_xy[:, 0] = s_xy[:, 1] # ds_xy/dy = 0 at y = 0

    s_xy[:, -1] = 0 # s_xy = 0 at y = D

    u[1:-1, 1:-1] += (epsilon_d / h2) * (s_xx[2:, 1:-1] - s_xx[:-2, 1:-1] + s_xy[1:-1, 2:] - s_xy[1:-1, :-2]) # x momentum

    u[0, :] = 0 # u = 0 at x = 0

    u[-1, :] = u[-2, :] # du/dx = 0 at x = L

    u[:, 0] = u[:, 1] # du/dy = 0 at y = 0

    u[:, -1] = u[:, -2] # du/dy = 0 at y = D    

    v[1:-1, 1:-1] += (epsilon_d / h2) * (s_xy[2:, 1:-1] - s_xy[:-2, 1:-1] + s_yy[1:-1, 2:] - s_yy[1:-1, :-2]) # y momentum

    v[0, :] = v[1, :] # v = 0 at x = 0

    v[-1, :] = v[-2, :] # dv/dx = 0 at x = L

    v[:, 0] = 0 # dv/dy = 0 at y = 0

    v[:, -1] = v[:, -2] # dv/dy = 0 at y = D

x = np.linspace(0, L, Nx) # initialize x

y = np.linspace(0, D, Ny) # initialize y

X, Y = np.meshgrid(x, y) # create mesh

plt.figure(dpi = 500) # make a nice crisp image :)

plt.contourf(X, Y, (u * 1000).T, cmap = 'jet', levels = 256) # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('u')

plt.xticks([])

plt.yticks([])

plt.show()

Flat Plate in Tension

     The case of flat plate in tension can be solved in very short time. The resultant displacement is shown in Fig. 1. The simulation takes advantage of symmetry and only a quarter of the plate simulated and shown within Fig. 1. The black plate represents The results of the code presented in this blog are compared with a not to be named commercial software 🤫 and the results are in close agreement with each other 🤝. The displacement 📏 from commercial 💰 software is at 0.00235 mm as compared to 0.00247 mm from the current code 🖳! As can be seen in Fig. 1, the plate expands towards the load and contracts in the other direction, as expected.

Fig. 1, The animation of plate displacement

Dog-Bone Specimen in Tension

     It is now possible to simulate a dog-bone shape specimens in pure tension. As is the case with this entire LE code, the code for simulation of stresses in the dog-bone specimen is also based on the FOANSS, the CFD 🌬 code developed by yours truly. The details on how to create the dog-bone shape are mentioned here. The results are compared with commercial, not to be named 🤫, code and the results are in very close agreement with each other 🤗. For example, the plate's maximum deflection in horizontal axis is at 0.00828 from the code as compared to a reading of 0.00895 mm from the commercial code❗ The von-mises stress and horizontal displacement is shown in Fig. 2.

Fig. 2, The results!

     If you want to hire me as you new shining post-doc or collaborate in research, do reach out.

Finite Difference Method for Linear Elasticity Update 03: Plates with Holes

      Welcome to my third blog! 👋 In this new adventure 👣, I try to solve linear elasticity 🩹 and fracture mechanics ⚡ problems and using the finite difference method. As with the other two, this blog is also created for my digital CV. The code shared here is heavily influenced by FOANSS 🌬️ and FOAMNE 🧠, the CFD and PINN codes 🖥️ developed by yours truly, in abundant spare time ⏳ of course. 😇 This finite difference based code can run easily on computers with even the modest specifications!

     The first case shared here is an example of a flat plate in plain strain under tensile loading. The boundary conditions are mentioned here. While using the finite difference method, artificial diffusion had to added to the iterative loop to stabilize the solution. As is the case with FOANSS, this code is also fully vectorized. Fig. 1 shows the results. The results can be compared to the results shown in Fig. 1 here. The results from a leading linear elasticity solver are also shown within Fig. 1, for a comparison. It can be seen that results are in excellent agreement with each other. For example, the maximum u displacement corresponding to the applied load is at 0.0043 mm as compared to 0.0044 mm from the commercial software that shall not be named 😁. 

Fig. 1, Current code VS commercial software

Code

#Copyright <2025> <FAHAD BUTT>

#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

#%% import necessary libraries

import numpy as np

import matplotlib.pyplot as plt

#%% define parameters

L = 1 # domain length

D = 0.5 # domain width

h = 1 / 50 # grid spacing

Nx = round(L / h) + 1 # grid points in x-axis

Ny = round(D / h) + 1 # grid points in y-axis

nu = 0.3 # Poisson's ratio

s_0 = 1e6 # applied load [Pa]

E = 2e11 # Young's modulus [Pa]

E_0 = E / s_0 # non dimensional Young's modulus

U = s_0 * L / E # scaled u displacement

V = s_0 * D / E # scaled v displacement

e_0 = s_0 / E # scaled strain

alpha = 1 / ((1 + nu) * (1 - (2 * nu))) # non dimensional parameter

beta = 1 / (2 * (1 + nu)) # non dimensional parameter

#%% initialize fields

u = np.zeros((Nx, Ny)) # u displacement

v = np.zeros((Nx, Ny)) # v displacement

s_xx = np.zeros((Nx, Ny)) # normal stress x

s_yy = np.zeros((Nx, Ny)) # normal stress y

s_xy = np.zeros((Nx, Ny)) # shear stress xy

e_xx = np.zeros((Nx, Ny)) # normal strain x

e_yy = np.zeros((Nx, Ny)) # normal strain y

e_xy = np.zeros((Nx, Ny)) # shear strain xy

#%% stability paremeters

relaxation = 0.1 # stabilization factor displacement

relaxation_I = relaxation / 10 # stabilization factor compatibility

epsilon_s = 1e-4 # relaxation factor artifitial diffusion (strain)

epsilon_d = 1e-7 # relaxation factor artifitial diffusion (displacement)

epsilon_r = 1e-10 # relaxation factor artifitial diffusion (compatibility)

#%% solve equations of motion

for nt in range(1000000):

    u[1:-1, 1:-1] += relaxation * h**2 * (D * ((s_xx[2:, 1:-1] - s_xx[:-2, 1:-1]) / (2 * h)) + L * ((s_xy[1:-1, 2:] - s_xy[1:-1, :-2]) / (2 * h))) # x momentum

    u[1:-1, 1:-1] += epsilon_d * (u[2:, 1:-1] - 2 * u[1:-1, 1:-1] + u[:-2, 1:-1]) / h**2  # artificial diffusion

    # boundary conditions for u

    u[0, :] = 0 # u = 0 at x = 0

    u[-1, :] = u[-2, :] # du/dx = 0 at x = L

    u[:, 0] = u[:, 1] # du/dy = 0 at y = 0

    u[:, -1] = u[:, -2] # du/dy = 0 at y = D

    v[1:-1, 1:-1] += relaxation * h**2 * (D * ((s_xy[2:, 1:-1] - s_xy[:-2, 1:-1]) / (2 * h)) + L * ((s_yy[1:-1, 2:] - s_yy[1:-1, :-2]) / (2 * h))) # y momentum

    v[1:-1, 1:-1] += epsilon_d * (v[2:, 1:-1] - 2 * v[1:-1, 1:-1] + v[:-2, 1:-1]) / h**2  # artificial diffusion

    # boundary conditions for v

    v[0, :] = 0 # v = 0 at x = 0

    v[-1, :] = v[-2, :] # dv/dx = 0 at x = L

    v[:, 0] = v[:, 1] # dv/dy = 0 at y = 0

    v[:, -1] = v[:, -2] # dv/dy = 0 at y = D

    e_xx[1:-1, 1:-1] = (u[2:, 1:-1] - u[:-2, 1:-1]) / (2 * h) # normal strain x

    e_xx[1:-1, 1:-1] += epsilon_s * (e_xx[2:, 1:-1] - 2 * e_xx[1:-1, 1:-1] + e_xx[:-2, 1:-1]) / h**2  # artificial diffusion

    # boundary conditions for e_xx

    e_xx[0, :] = e_xx[1, :] # de_xx/dx = 0 at x = 0

    e_xx[-1, :] = e_xx[-2, :] # de_xx/dx = 0 at x = L

    e_xx[:, 0] = e_xx[:, 1] # de_xx/dy = 0 at y = 0

    e_xx[:, -1] = e_xx[:, -2] # de_xx/dy = 0 at y = D

    exx_yy = (e_xx[1:-1, 2:] - 2 * e_xx[1:-1, 1:-1] + e_xx[1:-1, :-2]) / h**2

    e_yy[1:-1, 1:-1] = (v[1:-1, 2:] - v[1:-1, :-2]) / (2 * h) # normal strain y

    e_yy[1:-1, 1:-1] += epsilon_s * (e_yy[2:, 1:-1] - 2 * e_yy[1:-1, 1:-1] + e_yy[:-2, 1:-1]) / h**2  # artificial diffusion

    # boundary conditions for e_yy

    e_yy[0, :] = e_yy[1, :] # de_yy/dx = 0 at x = 0

    e_yy[-1, :] = e_yy[-2, :] # de_yy/dx = 0 at x = L

    e_yy[:, 0] = e_yy[:, 1] # de_yy/dy = 0 at y = 0

    e_yy[:, -1] = e_yy[:, -2] # de_yy/dy = 0 at y = D

    eyy_xx = (e_yy[2:, 1:-1] - 2 * e_yy[1:-1, 1:-1] + e_yy[:-2, 1:-1]) / h**2

    e_xy[1:-1, 1:-1] = (u[1:-1, 2:] - u[1:-1, :-2] + v[2:, 1:-1] - v[:-2, 1:-1]) / (2 * h) # shear strain xy

    e_xy[1:-1, 1:-1] += epsilon_s * (e_xy[2:, 2:] - 2 * e_xy[1:-1, 1:-1] + e_xy[:-2, :-2]) / h**2  # artificial diffusion

    # boundary conditions for e_xy

    e_xy[0, :] = e_xy[1, :] # de_xy/dx = 0 at x = 0

    e_xy[-1, :] = e_xy[-2, :] # de_xy/dx = 0 at x = L

    e_xy[:, 0] = e_xy[:, 1] # de_xy/dy = 0 at y = 0

    e_xy[:, -1] = e_xy[:, -2] # de_xy/dy = 0 at y = D

    exy_xy = (e_xy[2:, 2:] - e_xy[2:, :-2] - e_xy[:-2, 2:] + e_xy[:-2, :-2]) / (4 * h**2)

    R = (L / D) * exx_yy - 2 * exy_xy + (D / L) * eyy_xx # compatibility factor

    R[1:-1, 1:-1] += epsilon_r * (R[2:, 2:] - 2 * R[1:-1, 1:-1] + R[:-2, :-2]) / h**2 # artificial diffusion

    e_xx[1:-1, 1:-1] -= relaxation_I * h**2 * R * (D / L) # corrected normal strain x 

    e_yy[1:-1, 1:-1] -= relaxation_I * h**2 * R * (L / D) # corrected normal strain y

    e_xy[1:-1, 1:-1] -= relaxation_I * h**2 * R / 2 # corrected shear strain xy

    s_xx[1:-1, 1:-1] = alpha * ((1 - nu) * e_xx[1:-1, 1:-1] + nu * e_yy[1:-1, 1:-1]) # normal stress x

    # boundary conditions for s_xy

    s_xx[0, :] = s_xx[1, :] # ds_xx/dx = 0 at x = 0

    s_xx[-1, :] = 1 # s_xx = 1 at x = L

    s_xx[:, 0] = s_xx[:, 1] # ds_xx/dy = 0 at y = 0

    s_xx[:, -1] = s_xx[:, -2] # ds_xx/dy = 0 at y = D

    s_yy[1:-1, 1:-1] = alpha * ((1 - nu) * e_yy[1:-1, 1:-1] + nu * e_xx[1:-1, 1:-1]) # normal stress y

    # boundary conditions for s_xy

    s_yy[0, :] = s_yy[1, :] # ds_yy/dy = 0 at x = 0

    s_yy[-1, :] = s_yy[-2, :] # ds_yy/dy = 0 at x = L

    s_yy[:, 0] = 0 # s_yy = 0 at y = 0

    s_yy[:, -1] = 0 # s_yy = 0 at y = D

    s_xy[1:-1, 1:-1] = beta * e_xy[1:-1, 1:-1] # shear stress xy

    # boundary conditions for s_xy

    s_xy[0, :] = s_xy[1, :] # ds_xy/dx = 0 at x = 0

    s_xy[-1, :] = s_xy[-2, :] # ds_xy/dx = 0 at x = L

    s_xy[:, 0] = s_xy[:, 1] # ds_xy/dx = 0 at y = 0

    s_xy[:, -1] = s_xy[:, -2] # ds_xy/dx = 0 at y = D

#%% plotting

x = np.linspace(0, L, Nx) # initialize x

y = np.linspace(0, D, Ny) # initialize y

X, Y = np.meshgrid(x, y) # create combinations of sample points along the domain

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (u * U * 1000).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('u')

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (v * V * 1000).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('v')

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (s_xx * s_0 * 1e-6).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('s_xx')

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (s_yy * s_0 * 1e-6).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('s_yy')

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (s_xy * s_0 * 1e-6).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('s_xy')

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (e_xx * e_0).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('e_xx')

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (e_yy * e_0).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('e_yy')

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (e_xy * e_0).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('e_xy')

plt.show()

Update 01

     The code 💻 shared in the original post had different tricks added to stabilize the solution. These tricks allowed the code to run on relatively coarse meshes. The code as a result became to complicated with too many tunable parameters. Update 01 has the vanilla 🍦 finite difference code to solve linear elasticity problems ⭕. It takes slightly more time to run, and is more accurate. The results are shown in Fig. 2. Within Fig. 2, left to right in top row is u and v displacements, middle row is composed of stresses i.e. σ_xx, σ_yy and σ_xy. While the bottom row shows strains i.e. ε_xx, ε_yy and ε_xy.

     I completely understand that the finite element method is more suitable for linear elasticity problems, but this code has, well, a higher purpose 😱. The said purpose shall be revealed soon! Happy coding! 😎

Fig. 2, Results of the vanilla code

#Copyright <2025> <FAHAD BUTT>

#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

#%% import necessary libraries

import numpy as np

import matplotlib.pyplot as plt

#%% define parameters

L = 1 # domain length

D = 0.5 # domain width

h = 1 / 250 # grid spacing

Nx = round(L / h) + 1 # grid points in x-axis

Ny = round(D / h) + 1 # grid points in y-axis

nu = 0.3 # Poisson's ratio

s_0 = 1e6 # applied load [Pa]

E = 2e11 # Young's modulus [Pa]

E_0 = E / s_0 # non dimensional Young's modulus

U = s_0 * L / E # scaled u displacement

V = s_0 * D / E # scaled v displacement

e_0 = s_0 / E # scaled strain

alpha = 1 / ((1 + nu) * (1 - (2 * nu))) # non dimensional parameter

beta = 1 / (2 * (1 + nu)) # non dimensional parameter

#%% initialize fields

u = np.zeros((Nx, Ny)) # u displacement

v = np.zeros((Nx, Ny)) # v displacement

s_xx = np.zeros((Nx, Ny)) # normal stress x

s_yy = np.zeros((Nx, Ny)) # normal stress y

s_xy = np.zeros((Nx, Ny)) # shear stress xy

e_xx = np.zeros((Nx, Ny)) # normal strain x

e_yy = np.zeros((Nx, Ny)) # normal strain y

e_xy = np.zeros((Nx, Ny)) # shear strain xy

#%% stability paremeters

epsilon_d = 1e-5 # stability factor

#%% solve equations of motion

for nt in range(500000):

    e_xx[1:-1, 1:-1] = (u[2:, 1:-1] - u[:-2, 1:-1]) / (2 * h) # normal strain x

    # boundary conditions for e_xx

    e_xx[0, :] = e_xx[1, :] # de_xx/dx = 0 at x = 0

    e_xx[-1, :] = e_xx[-2, :] # de_xx/dx = 0 at x = L

    e_xx[:, 0] = e_xx[:, 1] # de_xx/dy = 0 at y = 0

    e_xx[:, -1] = e_xx[:, -2] # de_xx/dy = 0 at y = D

    e_yy[1:-1, 1:-1] = (v[1:-1, 2:] - v[1:-1, :-2]) / (2 * h) # normal strain y

    # boundary conditions for e_yy

    e_yy[0, :] = e_yy[1, :] # de_yy/dx = 0 at x = 0

    e_yy[-1, :] = e_yy[-2, :] # de_yy/dx = 0 at x = L

    e_yy[:, 0] = e_yy[:, 1] # de_yy/dy = 0 at y = 0

    e_yy[:, -1] = e_yy[:, -2] # de_yy/dy = 0 at y = D

    e_xy[1:-1, 1:-1] = 0.5 * (u[1:-1, 2:] - u[1:-1, :-2] + v[2:, 1:-1] - v[:-2, 1:-1]) / (2 * h) # shear strain xy

    # boundary conditions for e_xy

    e_xy[0, :] = e_xy[1, :] # de_xy/dx = 0 at x = 0

    e_xy[-1, :] = e_xy[-2, :] # de_xy/dx = 0 at x = L

    e_xy[:, 0] = e_xy[:, 1] # de_xy/dy = 0 at y = 0

    e_xy[:, -1] = e_xy[:, -2] # de_xy/dy = 0 at y = D

    

    s_xx[1:-1, 1:-1] = alpha * ((1 - nu) * e_xx[1:-1, 1:-1] + nu * e_yy[1:-1, 1:-1]) # normal stress x

    # boundary conditions for s_xy

    s_xx[0, :] = s_xx[1, :] # ds_xx/dx = 0 at x = 0

    s_xx[-1, :] = 1 # s_xx = 1 at x = L

    s_xx[:, 0] = s_xx[:, 1] # ds_xx/dy = 0 at y = 0

    s_xx[:, -1] = s_xx[:, -2] # ds_xx/dy = 0 at y = D

    s_yy[1:-1, 1:-1] = alpha * ((1 - nu) * e_yy[1:-1, 1:-1] + nu * e_xx[1:-1, 1:-1]) # normal stress y

    # boundary conditions for s_xy

    s_yy[0, :] = s_yy[1, :] # ds_yy/dy = 0 at x = 0

    s_yy[-1, :] = s_yy[-2, :] # ds_yy/dy = 0 at x = L

    s_yy[:, 0] = 0 # s_yy = 0 at y = 0

    s_yy[:, -1] = 0 # s_yy = 0 at y = D

    s_xy[1:-1, 1:-1] = beta * e_xy[1:-1, 1:-1] # shear stress xy

    # boundary conditions for s_xy

    s_xy[0, :] = s_xy[1, :] # ds_xy/dx = 0 at x = 0

    s_xy[-1, :] = s_xy[-2, :] # ds_xy/dx = 0 at x = L

    s_xy[:, 0] = s_xy[:, 1] # ds_xy/dx = 0 at y = 0

    s_xy[:, -1] = s_xy[:, -2] # ds_xy/dx = 0 at y = D

    

    u[1:-1, 1:-1] += epsilon_d * (D * ((s_xx[2:, 1:-1] - s_xx[:-2, 1:-1]) / (2 * h)) + L * ((s_xy[1:-1, 2:] - s_xy[1:-1, :-2]) / (2 * h))) # x momentum

    # boundary conditions for u

    u[0, :] = 0 # u = 0 at x = 0

    u[-1, :] = u[-2, :] # du/dx = 0 at x = L

    u[:, 0] = u[:, 1] # du/dy = 0 at y = 0

    u[:, -1] = u[:, -2] # du/dy = 0 at y = D

    v[1:-1, 1:-1] += epsilon_d * (D * ((s_xy[2:, 1:-1] - s_xy[:-2, 1:-1]) / (2 * h)) + L * ((s_yy[1:-1, 2:] - s_yy[1:-1, :-2]) / (2 * h))) # y momentum

    # boundary conditions for v

    v[0, :] = 0 # v = 0 at x = 0

    v[-1, :] = v[-2, :] # dv/dx = 0 at x = L

    v[:, 0] = v[:, 1] # dv/dy = 0 at y = 0

    v[:, -1] = v[:, -2] # dv/dy = 0 at y = D

#%% plotting

x = np.linspace(0, L, Nx) # initialize x

y = np.linspace(0, D, Ny) # initialize y

X, Y = np.meshgrid(x, y) # create combinations of sample points along the domain

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (u * U * 1000).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('u')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (v * V * 1000).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('v')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (s_xx * s_0 * 1e-6).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('s_xx')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (s_yy * s_0 * 1e-6).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('s_yy')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (s_xy * s_0 * 1e-6).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('s_xy')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (e_xx * e_0).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('e_xx')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (e_yy * e_0).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('e_yy')

plt.xticks([])

plt.yticks([])

plt.show()

plt.figure(dpi = 200)  # make a nice crisp image :)

plt.contourf(X, Y, (e_xy * e_0).T, cmap = 'jet', levels = 256)  # contour plot

plt.colorbar(orientation = 'vertical')

plt.gca().set_aspect('equal', adjustable = 'box')

plt.title('e_xy')

plt.xticks([])

plt.yticks([])

plt.show()

Update 02: Bending Problems

     The three point bending test is now possible. The results are shown in Fig. 3. Once again, the results are compared to a not be named commercial software 👻 and are in very close agreement 🤝.

Fig. 3, The results!!

     The three point bending test is now possible. The results are shown in Fig. 4. The results in Fig. 4 show the black undeformed beam along with a superimposed beam in deformed state. The applied boundary conditions are shown in Fig. 5. The code is also made available!

Fig. 4, The results!

Fig. 5, The boundary conditions!

#Copyright <2025> <FAHAD BUTT>

#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

import numpy as np

L = 1 # domain length

D = 0.2 # domain width

h = D / 100 # grid spacing

h2 = 2 * h # multiplication factor

Nx = round(L / h) # grid points in x-axis

Ny = round(D / h) # grid points in y-axis

nu = 0.3 # Poisson's ratio

s_0 = 1e6 # applied load [Pa]

E = 2e11 # Young's modulus [Pa]

alpha = E / ((1 + nu) * (1 - (2 * nu))) # non dimensional parameter

beta = E / (2 * (1 + nu)) # non dimensional parameter

u = np.zeros((Nx, Ny)) # u displacement

v = np.zeros((Nx, Ny)) # v displacement

s_xx = np.zeros((Nx, Ny)) # normal stress x

s_yy = np.zeros((Nx, Ny)) # normal stress y

s_xy = np.zeros((Nx, Ny)) # shear stress xy

e_xx = np.zeros((Nx, Ny)) # normal strain x

e_yy = np.zeros((Nx, Ny)) # normal strain y

e_xy = np.zeros((Nx, Ny)) # shear strain xy

epsilon_d = 1e-17 # stability factor

for nt in range(1000000):

    e_xx[1:-1, 1:-1] = (u[2:, 1:-1] - u[:-2, 1:-1]) / h2 # normal strain x

    e_xx[0, :] = e_xx[1, :] # de_xx/dx = 0 at x = 0

    e_xx[-1, :] = e_xx[-2, :] # de_xx/dx = 0 at x = L

    e_xx[:, 0] = e_xx[:, 1] # de_xx/dy = 0 at y = 0

    e_xx[:, -1] = e_xx[:, -2] # de_xx/dy = 0 at y = D    

    e_yy[1:-1, 1:-1] = (v[1:-1, 2:] - v[1:-1, :-2]) / h2 # normal strain y

    e_yy[0, :] = e_yy[1, :] # de_yy/dx = 0 at x = 0

    e_yy[-1, :] = e_yy[-2, :] # de_yy/dx = 0 at x = L

    e_yy[:, 0] = e_yy[:, 1] # de_yy/dy = 0 at y = 0

    e_yy[:, -1] = e_yy[:, -2] # de_yy/dy = 0 at y = D 

    e_xy[1:-1, 1:-1] = 0.5 * (u[1:-1, 2:] - u[1:-1, :-2] + v[2:, 1:-1] - v[:-2, 1:-1]) / h2 # shear strain xy

    e_xy[0, :] = e_xy[1, :] # de_xy/dx = 0 at x = 0

    e_xy[-1, :] = e_xy[-2, :] # de_xy/dx = 0 at x = L

    e_xy[:, 0] = e_xy[:, 1] # de_xy/dy = 0 at y = 0

    e_xy[:, -1] = e_xy[:, -2] # de_xy/dy = 0 at y = D

    s_xx[1:-1, 1:-1] = alpha * ((1 - nu) * e_xx[1:-1, 1:-1] + nu * e_yy[1:-1, 1:-1]) # normal stress x

    s_xx[0, :] = 0 # s_xx = 0 at x = 0

    s_xx[-1, :] = 0 # s_xx = 0 at x = L

    s_xx[:, 0] = s_xx[:, 1] # ds_xx/dy = 0 at y = 0

    s_xx[:, -1] = s_xx[:, -2] # ds_xx/dy = 0 at y = D

    s_yy[1:-1, 1:-1] = alpha * ((1 - nu) * e_yy[1:-1, 1:-1] + nu * e_xx[1:-1, 1:-1]) # normal stress y

    s_yy[0, :] = s_yy[1, :] # ds_yy/dx = 0 at x = 0

    s_yy[-1, :] = s_yy[-2, :] # ds_yy/dx = 0 at x = L

    s_yy[:, 0] = 0 # s_yy = 0 at y = 0

    s_yy[round(0.125 * Nx / L) - 1:round(0.125 * Nx / L) + 1, 0] = s_yy[round(0.125 * Nx / L) - 1:round(0.125 * Nx / L) + 1, 1] # ds_yy/dy = 0 at x = 0.125 and y = 0

    s_yy[round(0.875 * Nx / L) - 1:round(0.875 * Nx / L) + 1, 0] = s_yy[round(0.875 * Nx / L) - 1:round(0.875 * Nx / L) + 1, 1] # ds_yy/dy = 0 at x = 0.875 and y = 0

    s_yy[:, -1] = 0 # s_yy = 0 at y = D

    s_yy[round(0.5 * Nx / L) - 1:round(0.5 * Nx / L) + 1, -1] = -s_0 # s_yy = s_0 at y = D

    s_xy[1:-1, 1:-1] = beta * e_xy[1:-1, 1:-1] # shear stress xy

    s_xy[0, :] = 0 # s_xy = 0 at x = 0

    s_xy[-1, :] = 0 # s_xy = 0 at x = L

    s_xy[:, 0] = 0 # s_xy = 0 at y = 0

    s_xy[round(0.125 * Nx / L) - 1:round(0.125 * Nx / L) + 1, 0] = s_xy[round(0.125 * Nx / L) - 1:round(0.125 * Nx / L) + 1, 1] # ds_xy/dy = 0 at x = 0.125 and y = 0

    s_xy[round(0.875 * Nx / L) - 1:round(0.875 * Nx / L) + 1, 0] = s_xy[round(0.875 * Nx / L) - 1:round(0.875 * Nx / L) + 1, 1] # ds_xy/dy = 0 at x = 0.875 and y = 0

    s_xy[:, -1] = 0 # s_xy = 0 at y = D

    s_xy[round(0.5 * Nx / L) - 1:round(0.5 * Nx / L) + 1, -1] = s_xy[round(0.5 * Nx / L) - 1:round(0.5 * Nx / L) + 1, -2] # ds_xy/dy = s_0 at y = D

    u[1:-1, 1:-1] += (epsilon_d / h2) * (s_xx[2:, 1:-1] - s_xx[:-2, 1:-1] + s_xy[1:-1, 2:] - s_xy[1:-1, :-2]) # x momentum

    u[0, :] = u[1, :] # du/dx = 0 at x = 0

    u[-1, :] = u[-2, :] # du/dx = 0 at x = L

    u[:, 0] = u[:, 1] # du/dy = 0 at y = 0

    u[:, -1] = u[:, -2] # du/dy = 0 at y = D

    v[1:-1, 1:-1] += (epsilon_d / h2) * (s_xy[2:, 1:-1] - s_xy[:-2, 1:-1] + s_yy[1:-1, 2:] - s_yy[1:-1, :-2]) # y momentum

    v[0, :] = v[1, :] # dv/dx = 0 at x = 0

    v[-1, :] = v[-2, :] # dv/dx = 0 at x = L

    v[:, 0] = v[:, 1] # dv/dy = 0 at y = 0

    v[:, -1] = v[:, -2] # dv/dy = 0 at y = D

    v[round(0.125 * Nx / L) - 1:round(0.125 * Nx / L) + 1, 0] = 0 # v = 0 at x = 0.125 and y = 0

    v[round(0.875 * Nx / L) - 1:round(0.875 * Nx / L) + 1, 0] = 0 # v = 0 at x = 0.875 and y = 0

Update 03: Objects with Holes

     This case is about the solution of benchmark problem in solid mechanics, i.e. a think plate with hole. The code and methodology is inspired by CFD 🌬 of flow around a square cylinder 📦. The CFD code is called FOANSS, developed by yours truly over the years. The results are shown in Fig. 6. The plate has a square hole and the plate is being put in tension on the horizontal axis. The code predicts plate shrinking in vertical axis and extending in horizontal quite well. The results are, yet again, very close to a commercial 🤑 software that shall not be named 🤐❗❗. From a background in fluid dynamics, I plotted streamlines using displacements. The arrows point towards compression of the hole and vice-versa. The curved streamlines indicate plate bending.

Fig. 6, The results!

     If you want to hire me as you new shining post-doc or collaborate in research, do reach out.