AI for Mechanical Engineering: Fluid Mechanics
http://iailab.kaist.ac.kr/
Industrial AI Lab at KAIST
We want to design a Venturi meter using cvxpy to maximize the pressure at Point 1. Given the constraints on a, b, and c, find the maximum pressure at Point 1 and the corresponding values of a, b, and c. (Define the specific weight of JP-4 fuel as 9.807)
\begin{array}{I} & 10 \leq a-c \leq 30 \\ & 15 \leq a+b \leq 25 \\ & 25 \leq 3b+c \leq 45 \\ & a, b \geq 1 \\ & c \geq 6 \\ \end{array}
$$ \begin{array}{I} \max \quad &\gamma a + \gamma (b sin(30°)) + \frac{\gamma V^2 (1-\frac{c^2}{16})}{2g} \\ \text{subject to} \quad & 10 \leq a-c \leq 30 \\ & 15 \leq a+b \leq 25 \\ & 25 \leq 3b+c \leq 45 \\ & a, b \geq 1 \\ & c \geq 6 \\ \end{array} \quad\implies\quad \begin{array}{I} \min \quad &-\gamma(a + (b sin(30°)) + \frac{ V^2 (1-\frac{c^2}{16})}{2g}) \\ \text{subject to} \quad & 10 \leq a-c \leq 30 \\ & 15 \leq a+b \leq 25 \\ & 25 \leq 3b+c \leq 45 \\ & a, b \geq 1 \\ & c \geq 6 \\ \end{array} $$
$$ \begin{array}{I} \min \quad &-\gamma(a + (b sin(30°)) + \frac{ V^2 (1-\frac{c^2}{16})}{2g}) \\ \text{subject to} \quad & 10 \leq a-c \leq 30 \\ & 15 \leq a+b \leq 25 \\ & 25 \leq 3b+c \leq 45 \\ & a, b \geq 1 \\ & c \geq 6 \\ \end{array} \quad\implies\quad \begin{array}{I} \min_{X} \quad & -\gamma(X^T HX + f^T X) \\ \text{subject to} \quad & AX \leq b \\ & AX \geq b_{lower} \\ & X \geq lb \end{array} $$
import cvxpy as cvx
import numpy as np
x = cvx.Variable(3)
f = np.array([[1], [np.sin(np.deg2rad(30))], [0]])
H = np.array([
[0, 0, 0],
[0, 0, 0],
[0, 0, -10**2 / (16 * 2 * 9.8)]
])
objective = cvx.Minimize(-(cvx.quad_form(x, H) + f.T @ x + 10**2/(2*9.8)))
A = np.array([
[1, 0, -1],
[1, 1, 0],
[0, 3, 1]
])
b = np.array([30, 25, 45])
b_lower = np.array([10, 15, 25])
lb = np.array([1, 1, 6])
constraints = [A @ x <= b, A @ x >= b_lower, lb <= x ]
prob = cvx.Problem(objective, constraints)
result = prob.solve()
print("The optimal value is", result * - 9.807)
print("A solution [a, b, c] is")
print(x.value)
2. Dimensionality Reduction¶
2.1. SVD¶
Data Desciription
ALL.mat - 151 snapshots of u and v velocity as well as vorticity (UALL, VALL, VORTALL).
Datat set for fluid flow past at cylinder at RE = 100
Download data here
Given the data above, perform SVD to obtain U, S, and V. Then, perform truncated SVD to obtain a 5-rank approximation and compare the results with the original data. The data to be used for this SVD is 'flows_mat', and the code for this is provided below.
- Input shape: (89351 x 151)
- 449 x 199 dimensions for 151 time steps
Data Preprocessing and Visualization
from google.colab import drive
drive.mount('/content/drive')
from scipy import io
import numpy as np
flows = io.loadmat('/content/drive/MyDrive/Data/CYLINDER_ALL.mat')
flows_mat_u = flows['UALL']
flows_mat_v = flows['VALL']
flows_mat_vort = flows['VORTALL']
flows_mat_u.shape
nx = 449
ny = 199
flows_mat_u[:,0].reshape(nx, ny),
import matplotlib.pyplot as plt
plt.subplot(1,3,1)
plt.imshow(flows_mat_u[:,150].reshape(nx, ny), 'gray')
plt.title('u velocity')
plt.axis('off')
plt.subplot(1,3,2)
plt.imshow(flows_mat_v[:,150].reshape(nx, ny), 'gray')
plt.title('v velocity')
plt.axis('off')
plt.subplot(1,3,3)
plt.imshow(flows_mat_vort[:,150].reshape(nx, ny), 'gray')
plt.title('vorticity')
plt.axis('off')
plt.show()
from matplotlib.animation import FuncAnimation
import time
from IPython import display
fig = plt.figure()
ax1 = fig.add_subplot(1, 3, 1)
ax2 = fig.add_subplot(1, 3, 2)
ax3 = fig.add_subplot(1, 3, 3)
def animate(i):
ax1.clear()
ax1.imshow(flows_mat_u[:,i].reshape(nx, ny), 'gray')
ax1.set_title('u velocity')
ax1.axis('off')
ax2.clear()
ax2.imshow(flows_mat_v[:,i].reshape(nx, ny), 'gray')
ax2.set_title('v velocity')
ax2.axis('off')
ax3.clear()
ax3.imshow(flows_mat_vort[:,i].reshape(nx, ny), 'gray')
ax3.set_title('vorticity')
ax3.axis('off')
ani = FuncAnimation(fig, animate, frames = int(150/1), interval = 10)
video = ani.to_html5_video()
html = display.HTML(video)
display.display(html)
plt.close()
flows_mat_vort_mean = flows_mat_vort.mean(axis = 1)
flows_mat_vort_centered = flows_mat_vort - flows_mat_vort_mean[:, None]
flows_mat_vort_centered_normalized = (flows_mat_vort_centered - flows_mat_vort_centered.min()) / (flows_mat_vort_centered.max() - flows_mat_vort_centered.min())
plt.imshow(flows_mat_vort_mean.reshape(nx, ny), 'gray')
plt.title('vorticity')
plt.axis('off')
plt.show()
flows_mat = flows_mat_vort_centered_normalized
Perform SVD
U, S, VT = np.linalg.svd(flows_mat, full_matrices = False)
print(U.shape)
print(S.shape)
print(VT.shape)
plt.semilogy(S, '-o')
plt.ylabel('Singular Value', fontsize = 15)
plt.xlabel('Index', fontsize = 15)
plt.show()
plt.subplot(1,4,1)
plt.imshow(U[:,0].reshape(nx, ny), 'gray')
plt.title('POD Mode 1')
plt.axis('off')
plt.subplot(1,4,2)
plt.imshow(U[:,1].reshape(nx, ny), 'gray')
plt.title('POD Mode 2')
plt.axis('off')
plt.subplot(1,4,3)
plt.imshow(U[:,4].reshape(nx, ny), 'gray')
plt.title('POD Mode 5')
plt.axis('off')
plt.subplot(1,4,4)
plt.imshow(U[:,12].reshape(nx, ny), 'gray')
plt.title('POD Mode 13')
plt.axis('off')
plt.show()
U = np.asmatrix(U)
diagS = np.asmatrix(np.diag(S))
VT = np.asmatrix(VT)
print(U.shape, diagS.shape, VT.shape, '\n')
Perform Truncated SVD with Five Rank Approximation
k = 5
low_U = U[:, :k]
low_S = diagS[:k, :k]
low_VT = VT[:k, :]
low_flows_mat = low_U*low_S*low_VT
print(low_U.shape, low_S.shape, low_VT.shape, '\n')
from matplotlib.animation import FuncAnimation
import time
from IPython import display
fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)
def animate(i):
ax1.clear()
ax1.imshow(flows_mat[:,i].reshape(nx, ny), 'gray')
ax1.axis('off')
ax2.clear()
ax2.imshow(low_flows_mat[:,i].reshape(nx, ny), 'gray')
ax2.axis('off')
ani = FuncAnimation(fig, animate, frames = int(150/1), interval = 10)
video = ani.to_html5_video()
html = display.HTML(video)
display.display(html)
plt.close()