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)

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

  • 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

In [ ]:
from google.colab import drive
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

flows = io.loadmat('/content/drive/MyDrive/Data/CYLINDER_ALL.mat')
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()
No description has been provided for this image
In [ ]:
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()
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()
No description has been provided for this image
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()
No description has been provided for this image
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()
No description has been provided for this image
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()
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()