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