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!