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!
No comments:
Post a Comment