AI for Mechanical Engineering: Fluid Mechanics

By Jaejung Park
http://iailab.kaist.ac.kr/
Industrial AI Lab at KAIST

# 1. OptimizationÂ¶

## 1.1. Convex Optimization with CVXPYÂ¶

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}$$

InÂ [Â ]:
import cvxpy as cvx
import numpy as np

x = cvx.Variable(3)

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)
The optimal value is 151.57485714285716
A solution [a, b, c] is
[18.66666667  6.33333333  6.        ]

# 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

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

InÂ [Â ]:
drive.mount('/content/drive')
Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount("/content/drive", force_remount=True).
InÂ [Â ]:
from scipy import io
import numpy as np

InÂ [Â ]:
flows_mat_u = flows['UALL']
flows_mat_v = flows['VALL']
flows_mat_vort = flows['VORTALL']

flows_mat_u.shape
Out[Â ]:
(89351, 151)
InÂ [Â ]:
nx = 449
ny = 199

flows_mat_u[:,0].reshape(nx, ny),
Out[Â ]:
(array([[1.01243 , 1.01186 , 1.01119 , ..., 1.01816 , 1.01862 , 1.019   ],
[1.01394 , 1.01323 , 1.01255 , ..., 1.01956 , 1.02001 , 1.02046 ],
[1.01516 , 1.01458 , 1.01393 , ..., 1.02097 , 1.02139 , 1.02175 ],
...,
[1.01202 , 1.01168 , 1.01129 , ..., 0.888106, 0.891065, 0.894012],
[1.00918 , 1.00873 , 1.00829 , ..., 0.890758, 0.893679, 0.896936],
[1.00618 , 1.00579 , 1.00532 , ..., 0.893477, 0.896224, 0.898832]]),)
InÂ [Â ]:
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()
InÂ [Â ]:
from matplotlib.animation import FuncAnimation
import time
from IPython import display

fig = plt.figure()

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()
InÂ [Â ]:
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())
InÂ [Â ]:
plt.imshow(flows_mat_vort_mean.reshape(nx, ny), 'gray')
plt.title('vorticity')
plt.axis('off')
plt.show()
InÂ [Â ]:
flows_mat = flows_mat_vort_centered_normalized

Perform SVD

InÂ [Â ]:
U, S, VT = np.linalg.svd(flows_mat, full_matrices = False)
print(U.shape)
print(S.shape)
print(VT.shape)
(89351, 151)
(151,)
(151, 151)
InÂ [Â ]:
plt.semilogy(S, '-o')
plt.ylabel('Singular Value', fontsize = 15)
plt.xlabel('Index', fontsize = 15)
plt.show()
InÂ [Â ]:
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()
InÂ [Â ]:
U = np.asmatrix(U)
diagS = np.asmatrix(np.diag(S))
VT = np.asmatrix(VT)

print(U.shape, diagS.shape, VT.shape, '\n')
(89351, 151) (151, 151) (151, 151)

Perform Truncated SVD with Five Rank Approximation

InÂ [Â ]:
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')
(89351, 5) (5, 5) (5, 151)

InÂ [Â ]:
from matplotlib.animation import FuncAnimation
import time
from IPython import display

fig = plt.figure()

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()