Navier-Stokes equations¶
This section repeats the Robust PINN method presented by Peng et.al.
Steady 2D NS equations¶
The prototype problem of incompressible flow past a circular cylinder is considered.
The velocity vector is set to zero at all walls and the pressure is set to p = 0 at the outlet. The fluid density is taken as $\rho = 1kg/m^3$ and the dynamic viscosity is taken as $\mu = 2 · 10^{−2}kg/m^3$ . The velocity profile on the inlet is set as $u(0, y)=4 \frac{U_M}{H^2}(H-y) y$ with $U_M = 1m/s$ and $H = 0.41m$.
The two-dimensional steady-state Navier-Stokes equation is equivalently transformed into the following equations:
$$ \begin{equation} \begin{aligned} \sigma^{11} &=-p+2 \mu u_x \ \sigma^{22} &=-p+2 \mu v_y \ \sigma^{12} &=\mu\left(u_y+v_x\right) \ p &=-\frac{1}{2}\left(\sigma^{11}+\sigma^{22}\right) \ \left(u u_x+v u_y\right) &=\mu\left(\sigma_x^{11}+\sigma_y^{12}\right) \ \left(u v_x+v v_y\right) &=\mu\left(\sigma_x^{12}+\sigma_y^{22}\right) \end{aligned} \end{equation} $$
We construct a neural network with six outputs to satisfy the PDE constraints above:
$$ \begin{equation} u, v, p, \sigma^{11}, \sigma^{12}, \sigma^{22}=net(x, y) \end{equation} $$
Define Symbols and Geometric Objects¶
For the 2d problem, we define two coordinate symbolsx
andy
, six variables$ u, v, p, \sigma^{11}, \sigma^{12}, \sigma^{22}$ are defined.
The geometry object is a simple rectangle and circle with the operator -
.
x = Symbol('x')
y = Symbol('y')
rec = sc.Rectangle((0., 0.), (1.1, 0.41))
cir = sc.Circle((0.2, 0.2), 0.05)
geo = rec - cir
u = sp.Function('u')(x, y)
v = sp.Function('v')(x, y)
p = sp.Function('p')(x, y)
s11 = sp.Function('s11')(x, y)
s22 = sp.Function('s22')(x, y)
s12 = sp.Function('s12')(x, y)
Define Sampling Methods and Constraints¶
For the problem, three boundary conditions , PDE constraint and external data are presented. We use the robust-PINN model inspired by the traditional LAD (Least Absolute Derivation) approach, where the L1 loss replaces the squared L2 data loss.
@sc.datanode
class Inlet(sc.SampleDomain):
def sampling(self, *args, **kwargs):
points = rec.sample_boundary(1000, sieve=(sp.Eq(x, 0.)))
constraints = sc.Variables({'u': 4 * (0.41 - y) * y / (0.41 * 0.41)})
return points, constraints
@sc.datanode
class Outlet(sc.SampleDomain):
def sampling(self, *args, **kwargs):
points = geo.sample_boundary(1000, sieve=(sp.Eq(x, 1.1)))
constraints = sc.Variables({'p': 0.})
return points, constraints
@sc.datanode
class Wall(sc.SampleDomain):
def sampling(self, *args, **kwargs):
points = geo.sample_boundary(1000, sieve=((x > 0.) & (x < 1.1)))
#print("points3", points)
constraints = sc.Variables({'u': 0., 'v': 0.})
return points, constraints
@sc.datanode(name='NS_external')
class Interior_domain(sc.SampleDomain):
def __init__(self):
self.density = 2000
def sampling(self, *args, **kwargs):
points = geo.sample_interior(2000)
constraints = {'f_s11': 0., 'f_s22': 0., 'f_s12': 0., 'f_u': 0., 'f_v': 0., 'f_p': 0.}
return points, constraints
@sc.datanode(name='NS_domain', loss_fn='L1')
class NSExternal(sc.SampleDomain):
def __init__(self):
points = pd.read_csv('NSexternel_sample.csv')
self.points = {col: points[col].to_numpy().reshape(-1, 1) for col in points.columns}
self.constraints = {'u': self.points.pop('u'), 'v': self.points.pop('v'), 'p': self.points.pop('p')}
def sampling(self, *args, **kwargs):
return self.points, self.constraints
Define Neural Networks and PDEs¶
In the PDE definition part, we add these PDE nodes:
net = sc.MLP([2, 40, 40, 40, 40, 40, 40, 40, 40, 6], activation=sc.Activation.tanh)
net = sc.get_net_node(inputs=('x', 'y'), outputs=('u', 'v', 'p', 's11', 's22', 's12'), name='net', arch=sc.Arch.mlp)
pde1 = sc.ExpressionNode(name='f_s11', expression=-p + 2 * nu * u.diff(x) - s11)
pde2 = sc.ExpressionNode(name='f_s22', expression=-p + 2 * nu * v.diff(y) - s22)
pde3 = sc.ExpressionNode(name='f_s12', expression=nu * (u.diff(y) + v.diff(x)) - s12)
pde4 = sc.ExpressionNode(name='f_u', expression=u * u.diff(x) + v * u.diff(y) - nu * (s11.diff(x) + s12.diff(y)))
pde5 = sc.ExpressionNode(name='f_v', expression=u * v.diff(x) + v * v.diff(y) - nu * (s12.diff(x) + s22.diff(y)))
pde6 = sc.ExpressionNode(name='f_p', expression=p + (s11 + s22) / 2)
Define A Solver¶
Direct use of Adam optimization is less effective, so the LBFGS optimization method or a combination of both (Adam+LBFGS) is used for training:
s = sc.Solver(sample_domains=(Inlet(), Outlet(), Wall(), Interior_domain(), NSExternal()),
netnodes=[net],
init_network_dirs=['network_dir_adam'],
pdes=[pde1, pde2, pde3, pde4, pde5, pde6],
max_iter=300,
opt_config = dict(optimizer='LBFGS', lr=1)
)
The result is shown as follows:
Unsteady 2D N-S equations with unknown parameters¶
A two-dimensional incompressible flow and dynamic vortex shedding past a circular cylinder in a steady-state are numerically simulated. Respectively, the Reynolds number of the incompressible flow is $Re = 100$. The kinematic viscosity of the fluid is $\nu = 0.01$. The cylinder diameter D is 1. The simulation domain size is $[-15,25] × [-8,8]$. The computational domain is much smaller: $[1,8] × [-2,2]× [0,20]$.
$$ \begin{equation} \begin{aligned} &u_t+\lambda_1\left(u u_x+v u_y\right)=-p_x+\lambda_2\left(u_{x x}+u_{y y}\right) \ &v_t+\lambda_1\left(u v_x+v v_y\right)=-p_y+\lambda_2\left(v_{x x}+v_{y y}\right) \end{aligned} \end{equation} $$
where $\lambda_1$ and $\lambda_2$ are two unknown parameters to be recovered. We make the assumption that $u=\psi_y, \quad v=-\psi_x$
for some stream function $\psi(x, y)$. Under this assumption, the continuity equation will be automatically satisfied. The following architecture is used in this example,
$$ \begin{equation} \psi, p=net\left(t, x, y, \lambda_1, \lambda_2\right) \end{equation} $$
Define Symbols and Geometric Objects¶
We define three coordinate symbols x
, y
and t
, three variables $u,v,p$ are defined.
x = Symbol('x')
y = Symbol('y')
t = Symbol('t')
geo = sc.Rectangle((1., -2.), (8., 2.))
u = sp.Function('u')(x, y, t)
v = sp.Function('v')(x, y, t)
p = sp.Function('p')(x, y, t)
time_range = {t: (0, 20)}
Define Sampling Methods and Constraints¶
This example has only two equation constraints, while the former has six equation constraints. We also use the LAD-PINN model. Then the PDE constrained optimization model is formulated as:
$$ \min {\theta, \lambda} \frac{1}{# \mathbf{D}u} \sum{\left(t_i, x_i, u_i\right) \in \mathbf{D}u}\left|u_i-u\theta\left(t_i, x_i ; \lambda\right)\right|+\omega \cdot L{p d e} . $$
@sc.datanode(name='NS_domain', loss_fn='L1')
class NSExternal(sc.SampleDomain):
def __init__(self):
points = pd.read_csv('NSexternel_sample.csv')
self.points = {col: points[col].to_numpy().reshape(-1, 1) for col in points.columns}
self.constraints = {'u': self.points.pop('u'), 'v': self.points.pop('v'), 'p': self.points.pop('p')}
def sampling(self, *args, **kwargs):
return self.points, self.constraints
@sc.datanode(name='NS_external')
class NSEq(sc.SampleDomain):
def sampling(self, *args, **kwargs):
points = geo.sample_interior(density=2000, param_ranges=time_range)
constraints = {'continuity': 0, 'momentum_x': 0, 'momentum_y': 0}
return points, constraints
Define Neural Networks and PDEs¶
IDRLnet defines a network node to represent the unknown Parameters.
net = sc.MLP([3, 20, 20, 20, 20, 20, 20, 20, 20, 3], activation=sc.Activation.tanh)
net = sc.get_net_node(inputs=('x', 'y', 't'), outputs=('u', 'v', 'p'), name='net', arch=sc.Arch.mlp)
var_nr = sc.get_net_node(inputs=('x', 'y'), outputs=('nu', 'rho'), arch=sc.Arch.single_var)
pde = sc.NavierStokesNode(nu='nu', rho='rho', dim=2, time=True, u='u', v='v', p='p')
Define A Solver¶
Two nodes trained together
s = sc.Solver(sample_domains=(NSExternal(), NSEq()),
netnodes=[net, var_nr],
pdes=[pde],
network_dir='network_dir',
max_iter=10000)
Finally, the real velocity field and pressure field at t=10s are compared with the predicted results: