A GPU Accelerated Immersed Boundary Method for Linear Elasticity

     Immersed Boundary Methods (IBMs) have been extensively used for fluid-structure interaction (FSI) and fluid dynamics simulations. IBMs can also be implemented to solve linear elasticity problems. A SIBM is already developed by yours truly in abundant spare time, of course 😼. The developed SIBM is  now adapted to solve the linear-elasticity problems 〰️.

     SIBM offers several advantages compared to legacy finite element methods. With SIBM, there is no need to create expensive 🤑 and time-consuming ⏱︎  meshes for deforming boundaries. SIBM embeds objects of any arbitrary shapes on a simple Cartesian grid. The SIBM works well for domains with irregular or those boundaries which are evolving with time. In such cases legacy meshes become impractical.

     The ray casting algorithm🛸; a fundamental 🧱 technique used in video game 🎮 development and computer graphics, has been implemented in a similar manner to the fluid dynamics version. The detailed description of the method is already presented here.

     NOTE: As is the case with the fluid dynamics version, this code also requires a GPU to run. If dear readers don't have a GPU; a solution is to please stop being peasants... 🙉

Brazilian Test

     The first case presented is the case of indirect tensile test also referred to as the Brazilian test. The boundary conditions are described in Fig. 1. The simple boundary conditions include a load at the top in compression with a vertical support at the bottom. Due to the geometry of the specimen, the stress field center portion comes under tension. Meanwhile, the results from the simulations are presented in Fig. 2. It is clear that the SIBM can handle curved boundary problems with relative ease on a simple cartesian grid, without complicated the meshing.

Fig. 1, The boundary conditions

Fig. 2, The simulation results

Plate in Tension with a Hole

     This famous 📰 benchmark case of linear-elasticity 〰️ is now solvable with high accuracy. The challenge in this simple case is about application of the boundary conditions on the curved portion i.e. the hole. This is where IBM really shines. As the name suggests, the IBM immerses the hole inside grid. This makes application of the boundary 📏 conditions simple and straight forward. The boundary conditions are shown in Fig. 3. In the SIBM code, full plate is simulated. Note that this case has already been validated using a meshless 🕸️ method by yours truly 🤓. More details 📚 and free 💸 code  is available here.

Fig. 3, The boundary conditions

     The validation of this particular case has been extensively published. A not to be named 🤫, pre-validated commercial 🤑 code is used to validate the SIBM. The result from the commercial solver predict the von-Mises stress to be at 3.39 MPa. The SIBM indicate the von-Mises stress to be at, 🥁, 3.49 MPa ❗The resultant displacement from the 💰 code are at 0.000554 mm as compared to the SIBM value of 0.0006 mm. The resulting equivalent strain from the commercial code is at 1.31e-5. The strain from the SIBM is at 1.78e-5. The results from SIBM are plotted in Fig. 4. Within Fig. 4, the deformed plate is shown using red color while the undeformed plate is shown using black color. Please note that in the actual simulation, there are no moving nodes. This is another important advantage of SIBM which results in a much simple algorithm.

     The von-Mises stress, equivalent strain and resultant displacement (bottom) are shown in Fig. 5. The results shown in Fig. 5 are on an undeformed plate, for clarity. Within Fig. 5, the color red 🩷 means maximum and the color blue 💙 means minimum value. It can be seen 👓 that the code captures stress concentrations 🧐 with high accuracy.

Fig. 4, Deformed VS undeformed plate

Fig. 5, The results from post-processing

NAFEMS Elliptical Membrane

     The National Agency for Finite Element Methods and Standards is a long name 😫. Anyways 🙈, this case is the benchmark for testing irregular and curved 〜 boundaries domains, mixed 🥘 boundary conditions and symmetry 🪞. The SIBM  avoids time-consuming and complex meshes to capture curved 〰️ boundaries by embedding ☂️ the geometry in a regular cartesian grid. The result is accurate simulation of complex cases such as the symmetric elliptical ⬭ membrane. The peak stress of 98.05 MPa is predicted from the SIBM as compared to 92.7 MPa from the published literature, on a very coarse mesh 💠. The boundary conditions are shown in Fig. 6. The predicted normal stress σᵧᵧ is shown in Fig. 7. The deformed membrane (red) is also shown within Fig.7.

Fig. 6, The boundary conditions

Fig. 7, The results

     Thank you for reading! If you want to hire me as a post-doc researcher in the fields of thermo-fluids and / or fracture mechanics, do reach out!

GPU Accelerated Linear Elasticity Solver

      In abundant spare time ⏳, yours truly has implemented GPU accelerated linear elasticity solver [1]. The solver uses finite difference method on a cartesian grid. The first case solved is of  2D cylinder under compressive loading. The results are already validated and can be viewed here. The results are shown in Fig. 1. The code is provided to reproduce the results.

Fig. 1, post-processing of results

NOTE: This method requires a GPU, if dear readers don't have a GPU then please stop being  peasants... 🙉

     If you plan to use these codes in your scholarly work, do cite this blog as:

     Fahad Butt (2025). Accelerated Linear Elasticity (https://crackedfoundations.blogspot.com/2025/07/gpu-accelerated-linear-elasticity-solver.html), Blogger. Retrieved Month Date, Year

Code

#%% import necessary libraries

import cupy as cp

import matplotlib.pyplot as plt

#%% define parameters

L = 0.075 # domain length

D = 0.15 # domain width

h = L / 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

#%% initialize fields

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

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

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

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

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

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

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

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

#%% stability paremeters

epsilon_d = 1e-18 # 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]) / 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] = 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] = 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] = s_yy[:, 1] # ds_yy/dy = 0 at y = 0

    s_yy[:, -1] = -s_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/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

    # 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

    

    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

    # 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] = 0 # v = 0 at y = 0

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

   

#%% plotting

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

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

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

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

plt.contourf(X.get(), Y.get(), (u.get() * 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.get(), Y.get(), (v.get() * 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.get(), Y.get(), (s_xx.get() * 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.get(), Y.get(), (s_yy.get() * 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.get(), Y.get(), (s_xy.get() * 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.get(), Y.get(), e_xx.get().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.get(), Y.get(), e_yy.get().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.get(), Y.get(), e_xy.get().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()

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

References

[1] Ryosuke Okuta, Yuya Unno, Daisuke Nishino, Shohei Hido and Crissman Loomis. “CuPy: A NumPy-Compatible Library for NVIDIA GPU Calculations”. Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS), 2017

ncorr_2D_matlab (Compatibility with newer versions 2023+)

     This post is about running ncorr with newer versions of MATLAB. With help from [1], I edited the ncorr.m. The following section of the ncorr.m should be edited if MATLAB version is newer than 2023a. The code to be replaced starts at line 176 in the ncorr.m and continues till line 185. ncorr can be downloaded from [2].

Replace

% Initialize opengl ------------------------------------------%

% In earlier versions of Matlab this will fix inverted plots

% Plotting tools are also run based on opengl

if (ispc)

data_opengl = opengl('data');

if (~data_opengl.Software)

% Set opengl to software

opengl('software'); 

end

end

Replacement Code

if ispc

if verLessThan('matlab', '9.14') % R2023a or earlier

data_opengl = opengl('data');

if ~data_opengl.Software

opengl('software');

end

else

% In newer versions, opengl('data') is deprecated

info = rendererinfo(gca);

% Only set to software if info contains field 'GraphicsRenderer'

if isfield(info, 'GraphicsRenderer')

renderer_type = lower(info.GraphicsRenderer);

if ~contains(renderer_type, 'software')

opengl('software');

end

else

% If field is missing, assume hardware and force software

opengl('software');

end

end

end

References

[1] OpenAI. (2025). ChatGPT (June 2025 version) [Large language model].  https://chat.openai.com/chat

[2] https://github.com/justinblaber/ncorr_2D_matlab

2D Wave Equation Code (Finite Difference Method)

     After 🎉 successfully teaching 👨‍🏫 the neural network about the 〰️ wave equation, yours truly thought it is the time ⏱️ to write a simple code 🖥️ for discretized domain as well. A code yours truly created in abundant spare time is mentioned in this blog 📖.

     A complex (sinusoidal) 🌊 boundary condition is implemented. Dirichlet and Neumann boundary conditions are also implemented. The code is vectorized ↗ so there is only one loop 😁.  CFL condition is implemented in the time-step calculation so time-step ⏳ is adjusted based on the mesh size automatically 😇.

     For simple geometries, traditional numerical methods are is still better than PINNs. 🧠

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 = 1 # domain width

h = D / 100 # grid spacing

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

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

nu = 0.175 # Poisson's ratio

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

rho = 2350 # density [Kg/m3]

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

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

TT = travel * L / cp # total time

CFL = 0.1 # CFL number

dt = CFL * h / cp # time step

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

#%% initialize fields

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

u_vel = np.zeros((Nx, Ny)) # initial velocity

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

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

#%% pre calculate for speed

P1 = dt**2 * cp**2 / h**2

f = 2 / TT # frequency

omega = 2 * np.pi * f * dt

#%% 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

#%% solve 2D-transient wave equation

for nt in range(ns):

    u[1:-1, 1:-1] = 2 * un[1:-1, 1:-1] - unn[1:-1, 1:-1] + (dt**2 * cp**2 / h**2) * (un[2:, 1:-1] + un[:-2, 1:-1] + un[1:-1, 2:] + un[1:-1, :-2] - 4 * un[1:-1, 1:-1])

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

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

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

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

    u[int(0.5 * Nx / L) - 1:int(0.5 * Nx / L) + 2, int(0.5 * Ny / D) - 1:int(0.5 * Ny / D) + 2] = np.sin(omega * nt) # u = sin(2*pi*f*t) at x = 0.5 and y = 0.5

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

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

    if nt % 1000 == 0:

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

        plt.contourf(X, Y, u.T, cmap = 'prism', levels = 256, vmin = -0.25, vmax = 0.25) # contour plot

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

        plt.xticks([])

        plt.yticks([])

        plt.show()

     The code creates a pulse, as shown in Fig. 1. Within Fig. 1, 1 of several snapshots is shown. The material properties used are that of concrete.

Fig. 1, A nice wave

     Thank you for reading! If you want to hire me as a post-doc researcher in the fields of thermo-fluids and / or fracture mechanics, do reach out!

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.