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.
