AI for Mechanical Engineering: Solid Mechanics
http://iailab.kaist.ac.kr/
Industrial AI Lab at KAIST
1. Stress-strain Curve¶
Tensile Test
Tensile testing is a destructive test process that provides information about the tensile strength, yield strength, and ductility of the metallic material.
Stress-strain Curve
The tensile test result can be represented as a stress-strain curve, which shows the relationship between stress and strain.
Some of the major properties in stress-strain curve
- Young’s modulus: Relationship between tensile stress and axial strain in the linear elastic region
- Yield strength: Stress corresponding to the yield point
- Ultimate tensile strength: The maximum stress that a material can withstand while being stretched
How to find Young's modulus from the stress-strain curve
- Young's modulus is the slope of the stress-strain curve in the elastic region
- We can caculate the slope using linear regression
1.1. Linear Regression¶
Recall linear regression
$\text{Given} \; \begin{cases} x_{i} \; \text{: inputs} \\ y_{i} \; \text{: outputs} \end{cases}$ , Find $\theta_{0}$ and $\theta_{1}$
$$x= \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{m} \end{bmatrix}, \qquad y= \begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{m} \end{bmatrix} \approx \hat{y}_{i} = \theta_{0} + \theta_{1}x_{i}$$
$ \hat{y}_{i} $ : predicted output
in many cases, a linear model is used to predict $y_{i}$
$$ \hat{y}_{i} = \theta_{0} + \theta_{1}x_{i} \; \quad \text{ such that }\quad \min\limits_{\theta_{0}, \theta_{1}}\sum\limits_{i = 1}^{m} (\hat{y}_{i} - y_{i})^2$$
- Optimization problem
$$\min\limits_{\theta_{0}, \theta_{1}}\sum\limits_{i = 1}^{m} (\hat{y}_{i} - y_{i})^2 =\min\limits_{\theta}\lVert\Phi\theta-y\rVert^2_2 \qquad \qquad \left(\text{same as} \; \min_{x} \lVert Ax-b \rVert_2^2 \right)$$
$$ \text{solution} \; \theta^* = (\Phi^{T}\Phi)^{-1}\Phi^{T} y $$
1.2. Young’s Modulus with Linear Regression¶
Find Young's modulus from tensile test data using linear regression
Dataset Preparation
Download tensile test data for aluminum 6061-T651
import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model
data = np.loadtxt('T_020_A_1_001_022_03.csv', delimiter=',', skiprows=1)
data = data[:data.shape[0]-1,:]
plt.plot(data[:,0],data[:,1])
plt.xlabel('strain')
plt.ylabel('stress [MPa]')
plt.title('Stress–strain curve')
plt.grid('on')
Use linear algebra for linear regression
elastic_data = data[:101,:]
x = elastic_data[:,0].reshape(-1,1)
y = elastic_data[:,1].reshape(-1,1)
A = np.hstack([x**0, x])
A = np.asmatrix(A)
theta = (A.T*A).I*A.T*y
xp = np.arange(data[0,0], data[yield_index,0], 0.0001).reshape(-1, 1)
yp = theta[0,0] + theta[1,0]*xp
plt.plot(data[:,0],data[:,1])
plt.plot(xp, yp)
plt.xlabel('strain')
plt.ylabel('stress [MPa]')
plt.title('Stress–strain curve')
plt.grid('on')
plt.show()
print('Youngs modulus: {:.2f}[GPa]'.format(theta[1,0]/1000))
Use Scikit-Learn module for linear regression
elastic_data = data[:yield_index,:]
x = elastic_data[:,0].reshape(-1,1)
y = elastic_data[:,1].reshape(-1,1)
reg = linear_model.LinearRegression()
reg.fit(x, y)
xp = np.arange(data[0,0], data[yield_index,0], 0.0001).reshape(-1, 1)
plt.plot(data[:,0],data[:,1])
plt.plot(xp, reg.predict(xp))
plt.xlabel('strain')
plt.ylabel('stress [MPa]')
plt.title('Stress–strain curve')
plt.grid('on')
plt.show()
print('Youngs modulus: {:.2f}[GPa]'.format(reg.coef_[0][0]/1000))
1.3. Young’s Modulus and Temperature¶
data = np.loadtxt('T_200_A_1_094_041_14.csv', delimiter=',', skiprows=1)
data = data[:data.shape[0]-1,:]
plt.plot(data[:,0],data[:,1])
plt.xlabel('strain')
plt.ylabel('stress [MPa]')
plt.title('Stress–strain curve')
plt.grid('on')
Use linear algebra for linear regression
elastic_data = data[:37,:]
x = elastic_data[:,0].reshape(-1,1)
y = elastic_data[:,1].reshape(-1,1)
A = np.hstack([x**0, x])
A = np.asmatrix(A)
theta = (A.T*A).I*A.T*y
xp = np.arange(data[0,0], data[yield_index,0], 0.0001).reshape(-1, 1)
yp = theta[0,0] + theta[1,0]*xp
plt.plot(data[:,0],data[:,1])
plt.plot(xp, yp)
plt.xlabel('strain')
plt.ylabel('stress [MPa]')
plt.title('Stress–strain curve')
plt.grid('on')
plt.show()
print('Youngs modulus: {:.2f}[GPa]'.format(theta[1,0]/1000))
Use Scikit-Learn module for linear regression
elastic_data = data[:yield_index,:]
x = elastic_data[:,0].reshape(-1,1)
y = elastic_data[:,1].reshape(-1,1)
reg = linear_model.LinearRegression()
reg.fit(x, y)
xp = np.arange(data[0,0], data[yield_index,0], 0.0001).reshape(-1, 1)
plt.plot(data[:,0],data[:,1])
plt.plot(xp, reg.predict(xp))
plt.xlabel('strain')
plt.ylabel('stress [MPa]')
plt.title('Stress–strain curve')
plt.grid('on')
plt.show()
print('Youngs modulus: {:.2f}[GPa]'.format(reg.coef_[0][0]/1000))
2. Vibration Signal of Beam¶
2.1. Data Analysis¶
2.1.1. Data in Industry¶
Detection of abnormal data prevents major accidents by detecting defects in machines or structures in the industrial field
Non-destructive testing finds internal defects or cracks without damaging the inspected system
- Ultrasonic Testing
- Magnetic Particle Testing
- Eddy Current Testing
2.1.2. Vibration Signals of Beam¶
Consider a system composed of a beam and masses
The beam is reinforced with masses under four different conditions
(a) All six masses attached
(b) The mass at position L1 is removed
(c) Masses at positions L1 and L3 are removed
(d) Masses at positions L1, L3 and L5 are removed
Vibration Signal
- 70 experiments for each condition
- 6400 data points for each experiment
- Data size (70*4, 6400)
This high dimensional data requires dimension reduction
Let's use PCA dimension reduction!
2.2. PCA Dimension Reduction¶
Recall PCA dimension reduction
Dimension reduction without losing too much information
- Highly correlated data contains redundant features
- Find a low-dimensional, yet useful representation of the data
- Data $\rightarrow$ projection onto unit vector $\hat{u}_1$
Roughly speaking, PCA does a change of axes that can represent the data in a succinct manner
- Maximize variance (most separable)
- Minimize the sum-of-squares (minimum squared error)
Pictorial summary of PCA
from google.colab import drive
drive.mount('/content/drive')
import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
class_a = np.load('/content/drive/MyDrive/Colab Notebooks/AI for ME/수업자료/class_a.npy')
class_b = np.load('/content/drive/MyDrive/Colab Notebooks/AI for ME/수업자료/class_b.npy')
class_c = np.load('/content/drive/MyDrive/Colab Notebooks/AI for ME/수업자료/class_c.npy')
class_d = np.load('/content/drive/MyDrive/Colab Notebooks/AI for ME/수업자료/class_d.npy')
# class_a = np.load('class_a.npy')
# class_b = np.load('class_b.npy')
# class_c = np.load('class_c.npy')
# class_d = np.load('class_d.npy')
total_data = np.concatenate((class_a, class_b, class_c, class_d), axis = 0)
print(total_data.shape)
- Note that the dimension of the data is very high
label = np.zeros(280)
label[70:140] = 1
label[140:210] = 2
label[210:] = 3
print(label.shape)
# Plot sample data
i = np.random.randint(280)
plt.plot(total_data[i])
plt.xlabel('Frequency [Hz]')
plt.ylabel('Magnitude [dB]')
plt.title('Beam vibration data #{}'.format(i))
2.3.2. PCA: 1D¶
pca = PCA(n_components=1)
pca_data = pca.fit_transform(total_data)
plt.title('PCA Visualization')
plt.hist(pca_data[:70], bins=20, label='a', alpha=0.7)
plt.hist(pca_data[70:140], bins=20, label='b', alpha=0.7)
plt.hist(pca_data[140:210], bins=20, label='c', alpha=0.7)
plt.hist(pca_data[210:], bins=20, label='d', alpha=0.7)
plt.legend()
2.3.3. PCA: 2D¶
pca = PCA(n_components=2)
pca_data = pca.fit_transform(total_data)
plt.title('PCA Visualization')
plt.scatter(pca_data[:70,0], pca_data[:70,1], label='a')
plt.scatter(pca_data[70:140,0], pca_data[70:140,1], label='b')
plt.scatter(pca_data[140:210,0], pca_data[140:210,1], label='c')
plt.scatter(pca_data[210:,0], pca_data[210:,1], label='d')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend()
eigenvalues = pca.explained_variance_
print('Eigen value for PC1: {:.2f}'.format(eigenvalues[0]))
print('Eigen value for PC2: {:.2f}'.format(eigenvalues[1]))
3. PINN Beam¶
3.1. Deflection of Beams¶
From solid mechanics, beam deflection follows the second-order differential equation below
$${\partial^2 y \over \partial x^2} = {M(x) \over EI}$$
Where
- M: Bending moment
- E: Modulus of elasticity
- I : Moment of inertia
Under some boundary conditions, we can accurately calculate the deflection of the beam
3.2. Euler Beam Example¶
Recall PINN example about Euler Beam
- Problem properties
$$E = 1 \operatorname{pa}, \quad I = 1 \operatorname{kg\cdot m^2}, \quad q = 1 \operatorname{N/m}, \quad l = 1 \operatorname{m}$$
Make a neural network and loss functions like below :
PINN prediction
- Partial differential equations
$${\partial^2 y \over \partial x^2} = {M(x) \over EI}$$
- Exact solution
$$ y(x) = \begin{cases} -{2 \over 3}x^3+{1 \over 6}x, & x < {1 \over 2} \newline {2 \over 3}x^3-2x^2+{7 \over 6}x, & x \geq {1 \over 2} \end{cases} $$
3.3.2. Dynamic Analysis of Beam¶
- Force distribution
- Torque at $ x=0 $ is 0
- Bending moment
Can be written as
$$ M(x) = -4x + 8 (x-1/2) u(x-1/2) $$
Where $ u(x) $ is a unit step function.
- PDE for PINN
$${\partial^2 y \over \partial x^2} = {M(x) \over EI} = - 4x + 8 (x-1/2) u(x-1/2) $$
$${\partial^2 y \over \partial x^2} + 4x - 8 (x-1/2) u(x-1/2) = 0 $$
!pip install deepxde
import deepxde as dde
import numpy as np
from deepxde.backend.set_default_backend import set_default_backend
set_default_backend("tensorflow")
import tensorflow as tf
3.4.2. Define PDE System¶
# Define unit step function
def step(x,t):
step_1 = tf.constant(1., dtype=tf.float32)
step_0 = tf.constant(0., dtype=tf.float32)
step = tf.where(x>t, step_1, step_0)
return step
def ddy(x, y):
return dde.grad.hessian(y, x)
def pde(x, y):
dy_xx = ddy(x, y)
return dy_xx + 4*x*step(x,0) - 8*(x-1/2)*step(x,1/2)
3.4.3. Define Boundary Condition¶
def boundary_left(x, on_boundary):
return on_boundary and np.isclose(x[0], 0)
def boundary_right(x, on_boundary):
return on_boundary and np.isclose(x[0], 1)
# The fixed point (x=1/2) should be defined as an anchor point
observe_x = np.array([1/2])
observe_y = dde.icbc.PointSetBC(observe_x, 0, component=0)
geom = dde.geometry.Interval(0, 1)
bc1 = dde.DirichletBC(geom, lambda x: 0, boundary_left) # u(0) = 0
def true_solution(x):
return np.where(x < 0.5, -2*(x ** 3)/3 + x/6, 2*(x ** 3)/3 - 2*(x ** 2) + 7*x/6 - 1/6 )
data = dde.data.PDE(geom,
pde,
[bc1, observe_y],
num_domain = 10,
num_boundary = 2,
anchors=observe_x,
solution = true_solution,
num_test = 100)
3.4.4. Define and Train Model¶
layer_size = [1] + [20] * 3 + [1]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr=0.001)
losshistory, train_state = model.train(iterations = 5000)
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
5. PINN Catenary Cable¶
5.1. Catenary Cable¶
Catenary is a curve of an ideally hanging cable under it’s weight
5.1.1. Dynamic Analysis of Catenary¶
- For small interval $ds$
$$ \Sigma F_x = T(x+dx)\cos(\theta(x+dx)) - T(x)\cos(\theta(x)) = 0$$
$$ \Sigma F_y = T(x+dx)\sin(\theta(x+dx)) - T(x)\sin(\theta(x)) - \omega ds = 0$$
- From x-direction component
$$ T(x)\cos(\theta(x)) = T_0 $$
- From y-direction component
$$ T(x)\sin(\theta(x)) = \omega ds $$
- Therefore,
$$ {dy \over dx} = \tan(\theta(x)) = {\omega ds \over T_0} = {ds \over a} $$
- By substituting this into the formula for arc length
$$ {dy \over dx} = \sinh {x-x_0 \over a} $$
$$ y(x) = a \cosh {x-x_0 \over a} + C $$
5.1.2. Catenary Cable Example¶
Let's assume an example of catenary cable
- Cable properties $$ a = {T_0 \over \omega} = 0 $$
- Ragne of $x$ $$ -{1 \over 2} \leq x \leq 1 $$
- Boundary condition $$ y(-{1 \over 2}) = \cosh (-{1 \over 2}), \qquad y(1) = \cosh (1)$$
- PDE system $$ {dy \over dx} - \sinh x = 0 $$
- Exact solution $$ y(x) = \cosh (x) $$
!pip install deepxde
import deepxde as dde
import numpy as np
from deepxde.backend.set_default_backend import set_default_backend
set_default_backend("tensorflow")
import tensorflow as tf
import matplotlib.pyplot as plt
5.2.2. Define PDE System and Boundary Condition¶
def dy(x, y):
return dde.grad.jacobian(y, x)
def pde(x, y):
dy_x = dy(x, y)
return dy_x - tf.math.sinh(x)
def boundary_left(x, on_boundary):
return on_boundary and np.isclose(x[0], -1/2)
def boundary_right(x, on_boundary):
return on_boundary and np.isclose(x[0], 1)
geom = dde.geometry.Interval(-1/2, 1)
bc1 = dde.DirichletBC(geom, lambda x: np.cosh(-1/2), boundary_left) # u(0) = 0
bc2 = dde.DirichletBC(geom, lambda x: np.cosh(1), boundary_right) # u(0) = 0
def true_solution(x):
return np.cosh(x)
data = dde.data.PDE(geom,
pde,
[bc1, bc2],
num_domain = 10,
num_boundary = 2,
solution = true_solution,
num_test = 100)
5.2.3. Define and Train Model¶
layer_size = [1] + [20] * 3 + [1]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr=0.001)
losshistory, train_state = model.train(iterations = 5000)
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
5.2.4. Model Evaluation¶
test_x = np.array([[1/2]])
true_y = np.cosh(test_x)
pred_y = model.predict(test_x)
print('True: {0:.5f}'.format(true_y[0][0]))
print('Predict: {0:.5f}'.format(pred_y[0][0]))
print('Error: {0:.5f}'.format(true_y[0][0]-pred_y[0][0]))