Surrogate Modeling for Fluid Flows Based on Physics-Constrained Deep Learning Without Simulation Data

TL;DR

This paper introduces a physics-constrained deep learning approach to create surrogate models for fluid flows that do not require any simulation data, enabling efficient and accurate predictions in complex fluid dynamics problems.

Contribution

It proposes a novel structured deep neural network architecture that enforces physical laws and boundary conditions without relying on simulation data, advancing surrogate modeling techniques.

Findings

Accurately predicts fluid flow fields without simulation data.

Effectively propagates uncertainties in fluid properties and geometry.

Achieves excellent agreement with first-principle numerical simulations.

Abstract

Numerical simulations on fluid dynamics problems primarily rely on spatially or/and temporally discretization of the governing equation into the finite-dimensional algebraic system solved by computers. Due to complicated nature of the physics and geometry, such process can be computational prohibitive for most real-time applications and many-query analyses. Therefore, developing a cost-effective surrogate model is of great practical significance. Deep learning (DL) has shown new promises for surrogate modeling due to its capability of handling strong nonlinearity and high dimensionality. However, the off-the-shelf DL architectures fail to operate when the data becomes sparse. Unfortunately, data is often insufficient in most parametric fluid dynamics problems since each data point in the parameter space requires an expensive numerical simulation based on the first principle, e.g., Naiver--Stokes equations. In this paper, we provide a physics-constrained DL approach for surrogate modeling of fluid flows without relying on any simulation data. Specifically, a structured deep neural network (DNN) architecture is devised to enforce the initial and boundary conditions, and the governing partial differential equations are incorporated into the loss of the DNN to drive the training. Numerical experiments are conducted on a number of internal flows relevant to hemodynamics applications, and the forward propagation of uncertainties in fluid properties and domain geometry is studied as well. The results show excellent agreement on the flow field and forward-propagated uncertainties between the DL surrogate approximations and the first-principle numerical simulations.

Figures

Tables

Table 1. Table 1: Summary of learning/prediction performance of vascular flows with 20 20 20 parameter collocation points in parameter spaces for training

	Stenosis		Aneurysm
	Viscosity	Geometry	Viscosity	Geometry
Training Loss	$8\times {10}^{-5}$	$1\times {10}^{-3}$	$5.5\times {10}^{-5}$	$2.9\times {10}^{-5}$
Test error $e_u$	$8.18\times {10}^{-5}$	$9.61\times {10}^{-4}$	$9.2\times {10}^{-5}$	$1.38\times {10}^{-4}$
Test error $e_v$	$7.14\times {10}^{-8}$	$1.76\times {10}^{-6}$	$1.33\times {10}^{-7}$	$1.15\times {10}^{-6}$
Test error $e_p$	$2.33\times {10}^{-5}$	$2.23\times {10}^{-3}$	$5.81\times {10}^{-6}$	$2.09\times {10}^{-5}$

Table 2. Table 2: Converged loss with different penalty coefficient λ 𝜆 \lambda

Loss	$\lambda =1$	$\lambda =10$	$\lambda =100$	$\lambda =1000$
Boundary Condition	$1\times {10}^{-3}$	$2\times {10}^{-4}$	$7\times {10}^{-6}$	$1\times {10}^{-7}$
$x$-Momentum Equation	$1\times {10}^{-4}$	$3\times {10}^{-3}$	$8\times {10}^{-3}$	$1\times {10}^{-2}$

Table 3. Table 3: Training and testing performance for stenosis case with varying viscosity/geometry

Stenosis with varying viscosity (physics-constrained data-free)
Training viscosity points	10	20	30	40	50
Training cost on GPU	7598s	12977s	25988s	28543s	31996s
Training minibatch loss	$3.11\times {10}^{-4}$	$8\times {10}^{-5}$	$7.64\times {10}^{-6}$	$4.6\times {10}^{-5}$	$5.8\times {10}^{-5}$
Test error $e_u$	$4.13\times {10}^{-4}$	$8.18\times {10}^{-5}$	$6.02\times {10}^{-5}$	$3.29\times {10}^{-5}$	$7.27\times {10}^{-5}$
Test error $e_v$	$3.32\times {10}^{-7}$	$7.14\times {10}^{-8}$	$6.43\times {10}^{-8}$	$3.86\times {10}^{-8}$	$1.41\times {10}^{-7}$
Test error $e_p$	$6.08\times {10}^{-5}$	$2.33\times {10}^{-5}$	$9.82\times {10}^{-6}$	$6.87\times {10}^{-6}$	$1.41\times {10}^{-5}$
Stenosis with varying viscosity (purely data-driven)
Training viscosity points	10	20	30	40	50
Training cost on GPU	6900s	11967s	19982s	22497s	30893s
Simulation cost on CPU	400s	800s	1200s	1600s	2000s
Training minibatch loss	$5\times {10}^{-7}$	$5\times {10}^{-7}$	$1\times {10}^{-7}$	$1\times {10}^{-7}$	$1\times {10}^{-7}$
Test error $e_u$	$1.49\times {10}^{-5}$	$6.83\times {10}^{-6}$	$2.43\times {10}^{-6}$	$1.33\times {10}^{-6}$	$2.64\times {10}^{-6}$
Test error $e_v$	$7.84\times {10}^{-8}$	$2.37\times {10}^{-7}$	$8.84\times {10}^{-8}$	$3.09\times {10}^{-8}$	$3.01\times {10}^{-8}$
Test error $e_p$	$6.42\times {10}^{-7}$	$7.13\times {10}^{-7}$	$1.56\times {10}^{-6}$	$3.72\times {10}^{-7}$	$1.65\times {10}^{-6}$
Stenosis with varying geometry (physics-constrained data-free)
Training viscosity points	10	20	30	40	50
Training cost on GPU	6458s	11688s	18239s	24674s	42814s
Training minibatch loss	$5.0\times {10}^{-4}$	$1.0\times {10}^{-3}$	$3.0\times {10}^{-4}$	$7.0\times {10}^{-4}$	$1.0\times {10}^{-4}$
Test error $e_u$	$1.54\times {10}^{-3}$	$9.61\times {10}^{-4}$	$5.24\times {10}^{-4}$	$3.46\times {10}^{-4}$	$4.81\times {10}^{-4}$
Test error $e_v$	$2.93\times {10}^{-6}$	$1.76\times {10}^{-6}$	$1.26\times {10}^{-6}$	$7.43\times {10}^{-7}$	$9.53\times {10}^{-7}$
Test error $e_p$	$3.85\times {10}^{-3}$	$2.23\times {10}^{-3}$	$1.89\times {10}^{-3}$	$1.57\times {10}^{-3}$	$1.33\times {10}^{-3}$
Stenosis with varying geometry (purely data-driven)
Training viscosity points	10	20	30	40	50
Training cost on GPU	5498s	11238s	16285s	21848s	26396s
Simulation cost on CPU	400s	800s	1200s	1600s	2000s
Training minibatch loss	$3\times {10}^{-7}$	$1\times {10}^{-7}$	$2\times {10}^{-7}$	$1\times {10}^{-7}$	$1\times {10}^{-7}$
Test error $e_u$	$3.21\times {10}^{-5}$	$2.18\times {10}^{-5}$	$2.4\times {10}^{-5}$	$2.74\times {10}^{-5}$	$3.12\times {10}^{-5}$
Test error $e_v$	$2.11\times {10}^{-7}$	$1.37\times {10}^{-7}$	$6.92\times {10}^{-8}$	$9.33\times {10}^{-8}$	$5.21\times {10}^{-8}$
Test error $e_p$	$2.72\times {10}^{-6}$	$2.88\times {10}^{-7}$	$1.22\times {10}^{-6}$	$3.41\times {10}^{-7}$	$6.40\times {10}^{-7}$

Table 4. Table 4: Training and testing performance for aneurysm case with varying viscosity/geometry

Aneurysm with varying viscosity (physics-constrained data-free)
Training viscosity points	10	20	30	40	50
Training cost on GPU	5731s	14490s	25350s	24810s	31595s
Training minibatch loss	$9.71\times {10}^{-5}$	$5.50\times {10}^{-5}$	$1.07\times {10}^{-5}$	$2.33\times {10}^{-5}$	$2.22\times {10}^{-5}$
Test error $e_u$	$2.72\times {10}^{-4}$	$9.2\times {10}^{-5}$	$8.33\times {10}^{-5}$	$6.4\times {10}^{-5}$	$4.67\times {10}^{-5}$
Test error $e_v$	$2.22\times {10}^{-7}$	$1.33\times {10}^{-7}$	$1.31\times {10}^{-7}$	$5.96\times {10}^{-8}$	$7.46\times {10}^{-8}$
Test error $e_p$	$6.14\times {10}^{-7}$	$1.18\times {10}^{-5}$	$5.81\times {10}^{-6}$	$1.12\times {10}^{-5}$	$4.75\times {10}^{-6}$
Aneurysm with varying viscosity (purely data-driven)
Training viscosity points	10	20	30	40	50
Training cost on GPU	6851s	11861s	24580s	24353s	31116s
Simulation cost on CPU	400s	800s	1200s	1600s	2000s
Training minibatch loss	$5\times {10}^{-7}$	$5\times {10}^{-7}$	$1\times {10}^{-7}$	$1\times {10}^{-7}$	$1\times {10}^{-7}$
Test error $e_u$	$1.76\times {10}^{-4}$	$5.11\times {10}^{-6}$	$8.51\times {10}^{-6}$	$5.30\times {10}^{-7}$	$3.20\times {10}^{-6}$
Test error $e_v$	$1.07\times {10}^{-7}$	$6.31\times {10}^{-8}$	$1.56\times {10}^{-7}$	$6.95\times {10}^{-8}$	$7.42\times {10}^{-8}$
Test error $e_p$	$6.14\times {10}^{-7}$	$2.38\times {10}^{-7}$	$1.81\times {10}^{-6}$	$1.06\times {10}^{-6}$	$3.28\times {10}^{-7}$
Aneurysm with varying geometry (physics-constrained data-free)
Training viscosity points	10	20	30	40	50
Training cost on GPU	7755s	12904s	24645s	26724s	29515s
Training minibatch loss	$8.8\times {10}^{-5}$	$2.9\times {10}^{-5}$	$1.8\times {10}^{-5}$	$7.8\times {10}^{-6}$	$1.9\times {10}^{-5}$
Test error $e_u$	$3.07\times {10}^{-4}$	$1.38\times {10}^{-4}$	$1.29\times {10}^{-4}$	$1.31\times {10}^{-4}$	$1.19\times {10}^{-4}$
Test error $e_v$	$1.67\times {10}^{-6}$	$1.15\times {10}^{-6}$	$1.01\times {10}^{-6}$	$9.95\times {10}^{-7}$	$9.06\times {10}^{-7}$
Test error $e_p$	$3.77\times {10}^{-5}$	$2.09\times {10}^{-5}$	$1.28\times {10}^{-5}$	$8.99\times {10}^{-6}$	$1.02\times {10}^{-5}$
Aneurysm with varying geometry (purely data-driven)
Training viscosity points	10	20	30	40	50
Training cost on GPU	6322s	11607s	18582s	23309s	32124s
Simulation cost on CPU	400s	800s	1200s	1600s	2000s
Training minibatch loss	$3\times {10}^{-7}$	$2\times {10}^{-7}$	$4\times {10}^{-7}$	$1\times {10}^{-7}$	$8\times {10}^{-8}$
Test error $e_u$	$7.9\times {10}^{-5}$	$8.63\times {10}^{-5}$	$7.74\times {10}^{-5}$	$8.69\times {10}^{-5}$	$6.73\times {10}^{-5}$
Test error $e_v$	$8.77\times {10}^{-7}$	$5.34\times {10}^{-7}$	$5.03\times {10}^{-7}$	$6.24\times {10}^{-7}$	$6.97\times {10}^{-7}$
Test error $e_p$	$1.06\times {10}^{-5}$	$8.94\times {10}^{-6}$	$8.96\times {10}^{-6}$	$9.24\times {10}^{-6}$	$1.26\times {10}^{-5}$

Equations

$$\mathcal{F}(\mathbf{u},p)=0:=\{\begin{array}{rlr}& \nabla \cdot \mathbf{u}=0,& \mathbf{x},t\in {\Omega }_{f,t},\mathbf{\theta }\in {\mathbb{R}}^d,\\ & \frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}+\frac{1}{\rho }\nabla p-\nu {\nabla }^2\mathbf{u}+{\mathbf{b}}_f=0,& \mathbf{x},t\in {\Omega }_{f,t},\mathbf{\theta }\in {\mathbb{R}}^d\end{array}$$

$$\mathcal{I}(\mathbf{x},p,\mathbf{u},\mathbf{\theta })$$

$$\mathcal{B}(t,\mathbf{x},p,\mathbf{u},\mathbf{\theta })$$

$${\mathbf{z}}_l={\sigma }_l({\mathbf{W}}_l^T{\mathbf{z}}_{l-1}+{\mathbf{b}}_l),$$

$$\mathbf{f}(t,\mathbf{x},\mathbf{\theta })\approx {\mathbf{f}}^p(t,\mathbf{x},\mathbf{\theta })\triangleq {\mathbf{z}}_L(t,\mathbf{x},\mathbf{\theta };\mathbf{W},\mathbf{b}),$$

$${\mathcal{L}}_{data}(\mathbf{W},\mathbf{b})$$

$${\mathbf{W}}^{\ast },{\mathbf{b}}^{\ast }$$

$${\mathcal{L}}_{phy}(\mathbf{W},\mathbf{b})$$

$${\mathbf{W}}^{\ast }=\underset{\mathbf{W}}{\mathrm{argmin}}{\mathcal{L}}_{phy}(\mathbf{W}),$$

$$\displaystyle\mathrm{s.t.}\left\{\begin{split}&\mathcal{I}(\mathbf{x},p^{p},\mathbf{u}^{p},\boldsymbol{\theta})=0,\qquad&t=0,\mathrm{in}\ \Omega_{f},\\ &\mathcal{B}(t,\mathbf{x},p^{p},\mathbf{u}^{p},\boldsymbol{\theta})=0,&\mathrm{on}\ \partial\Omega_{f,t},\end{split}\right.$$

$${\mathcal{L}}_{phy}^c(\mathbf{W},\mathbf{b},{\lambda }_i,{\lambda }_b)=\underset{\mathrm{Equation}\mathrm{loss}}{\underbrace{{\mathcal{L}}_{phy}(\mathbf{W},\mathbf{b})}}+\underset{\mathrm{Initial}\mathrm{loss}}{\underbrace{{\lambda }_i{\Vert \mathcal{I}(\mathbf{x},p^p,{\mathbf{u}}^p,\mathbf{\theta })\Vert }_{{\Omega }_{f,t}}}}+\underset{\mathrm{Boundary}\mathrm{loss}}{\underbrace{{\lambda }_b{\Vert \mathcal{B}(t,\mathbf{x},p^p,{\mathbf{u}}^p,\mathbf{\theta })\Vert }_{\partial {\Omega }_{f,t}}}},$$

$$\begin{array}{rl}{\hat{\mathbf{u}}}^p(t,\mathbf{x},\mathbf{\theta };\mathbf{W},\mathbf{b})& ={\mathbf{u}}_{particular}(t,\mathbf{x},\theta )+D(t,\mathbf{x},\mathbf{\theta }){\mathbf{u}}^p(t,\mathbf{x},\mathbf{\theta };\mathbf{W},\mathbf{b}),\\ {\hat{p}}^p(t,\mathbf{x},\mathbf{\theta };\mathbf{W},\mathbf{b})& =p_{particular}(t,\mathbf{x},\theta )+D(t,\mathbf{x},\mathbf{\theta })p^p(t,\mathbf{x},\mathbf{\theta };\mathbf{W},\mathbf{b}),\end{array}$$

$$\begin{array}{rl}{\mathcal{L}}_{phy}^c(\mathbf{W},\mathbf{b})=& \underset{\mathrm{Mass}\mathrm{conservation}}{\underbrace{{\Vert \nabla \cdot {\hat{\mathbf{u}}}^p\Vert }_{{\Omega }_f}}}+\underset{\mathrm{Momentum}\mathrm{conservation}}{\underbrace{\Vert \frac{\partial {\hat{\mathbf{u}}}^p}{\partial t}+({\hat{\mathbf{u}}}^p\cdot \nabla ){\hat{\mathbf{u}}}^p+\frac{1}{\rho }\nabla {\hat{p}}^p-\nu {\nabla }^2{\hat{\mathbf{u}}}^p+\mathbf{b}{\Vert }_{{\Omega }_f}}}\\ & +\underset{\mathrm{Newman}\mathrm{boundary}\mathrm{loss}}{\underbrace{{\Vert {\mathcal{B}}_N(t,\mathbf{x},{\hat{p}}^p,{\hat{\mathbf{u}}}^p,\mathbf{\theta })\Vert }_{\partial {\Omega }_f}}},\end{array}$$

$$\begin{array}{rl}{\mathbf{W}}^{\ast },{\mathbf{b}}^{\ast }& =\underset{\mathbf{W},\mathbf{b}}{\mathrm{argmin}}{\mathcal{L}}_{phy}^c(\mathbf{W},\mathbf{b}).\end{array}$$

$${\hat{u}}^p=\left(\frac{d^2}{4}-y^2\right)u^p,{\hat{v}}^p=\left(\frac{d^2}{4}-y^2\right)v^p,$$

$${\hat{p}}^p=\frac{x-x_{in}}{x_{out}-x_{in}}{p}_{out}+\frac{x_{out}-x}{x_{out}-x_{in}}{p}_{in}+(x-x_{in})(x_{out}-x)p^p,$$

$$u_a=\frac{\Delta p}{2\nu \rho L}\left(\frac{d^2}{4}-y^2\right),$$

$$R(x)=R_0-A\frac{1}{\sqrt{(2\pi {\sigma }^2)}}\exp (-\frac{(x-\mu {)}^2}{2{\sigma }^2})$$

$$e_{\mathbf{u}}=\frac{||{\mathbf{u}}_{DNN}-{\mathbf{u}}_{CFD}|{|}_2}{\Delta p}$$

$$e_p=\frac{||p_{DNN}-p_{CFD}|{|}_2}{(\Delta p{)}^2}$$

$${\mathcal{L}}_{phy}^c={\mathcal{L}}_{phy}+\lambda {\mathcal{L}}_B,$$

Full text

††journal: Elsevier

1 Introduction

Complex fluids are ubiquitous in natural and industrial processes, and accurately simulating the fluid flows is indispensable in many disciplines, e.g., aerospace, civil, and biomedical engineering. A fluid system is typically governed by the Navier–Stokes equations, which is a highly nonlinear partial differential equation (PDE) system. Numerical simulation on fluid dynamics problems primarily relies on solving the PDE systems in a discretized form using, e.g., finite difference (FD), finite volume (FV), or finite element (FE) methods, which is known as the computational fluid dynamics (CFD) approach [1]. However, CFD simulations are often computationally cumbersome, especially for the flows with turbulence and complex geometries. Moreover, mesh generation also usually incurs a huge burden, in particular when moving boundary or large geometric variation is considered. The considerable computational expense greatly limits the use of principled CFD model to real-time predictions and many-query analysis, which are highly relevant to many scientific problems and real-life applications, e.g., timely clinical diagnosis and surgery planning for cardiovascular diseases, optimization design of aircraft with large parameter variations, and forward/inverse uncertainty quantification (UQ) in high-consequence systems. As an alternative, a cost-effective surrogate model is a computationally feasible way to tackle the aforementioned challenges.

A surrogate model only approximates the input-output relation of a system, which can be evaluated efficiently. Namely, given input parameters, e.g., initial/boundary/operational conditions, the quantities of interest (QoIs), such as velocity, pressure, shear stress, and their integrals can be obtained rapidly without conducting the principled CFD simulations. The existing surrogate modeling approaches can be roughly categorized into two classes: projection-based reduced order models (ROMs) and data-fit models [2]. In projection-based ROMs, a reduced basis is extracted from the simulation data using an unsupervised learning technique, e.g., proper orthogonal decomposition (POD, also known as principal component analysis) [3], and the full-order PDE operator is projected onto the subspace spanned by the reduced basis. As a result, the degrees of freedom of the system can be significantly reduced, and meanwhile, the underlying structure of the full-order model can be retained to a certain extent. Although holding some promises, the current projection-based ROM techniques have had limited impact on complex fluid dynamic problems mainly because of the stability and robustness issues [4, 5]. Moreover, projection-based ROMs are highly code-intrusive and their speedup potential is limited when strong nonlinearity exists [6] though several remedies such as sparse sampling [7, 8] exist. Another way to enable rapid simulations is to build a data-fit model, where a response surface of the system is learned from the simulation data in a supervised manner. Namely, a deterministic or probabilistic input-output mapping is constructed using, e.g., polynomial basis functions [9], radial basis functions [10], Gaussian process (GP) [11, 12], and stochastic polynomial chaos expansion (PCE) [13, 14, 15], among others. All these models are built upon the CFD solutions of selected collocation points in parameter space without the need to modify the codes of the CFD solver. Because of the non-intrusive feature and ease of implementation, data-fit surrogates have been used for a wide range of forward and inverse uncertainty quantification (UQ) problems in fluid dynamics [14, 16, 17]. However, traditional data-fit models have a hard time handling the problems with strong nonlinearities and high dimensionality. Deep learning (DL), in particular, the deep neural network (DNN) has become a popular surrogate modeling approach and has shown great potential to deal with high-dimensional nonlinear UQ problems [18, 19, 20]. It has been shown that DNN as a universal function approximator [21] can overcome the curse of dimensionality in certain problems [22, 23, 24]. In broader scientific computing and physical modeling communities, machine learning (ML) has been receiving a lot of attentions [25, 26, 27, 28, 29, 30, 31]. However, the tremendous success of DL in the computer science, witnessed in areas of computer vision and image recognition [32], can be mainly attributed to the availability of large-scale labeled data (i.e., “big data”) and the interpolatory nature of their problems. Unfortunately, labeled data for surrogate modeling of fluid systems are often sparse and could be noisy, since they are obtained from either principled CFD simulations or experimental observations, both of which are expensive to obtain. Therefore, in such “small data” regimes [33], the true power of DL cannot be fully exploited by naively using the off-the-shelf DL model in the computer science community as an end-to-end fashion [34] for a data-fit of surrogate modeling.

In conventional ML problems, the mechanism behind the system is usually unknown and thus can only be learned from the labeled data. In contrast, for modeling a physical system, the governing equations are usually known a priori but are difficult to solve efficiently. Instead of learning solely from the labeled data, e.g., solution of the states on certain points in parameter space, the known governing equations can be utilized to constrain (or even drive) the learning to compensate for the insufficiency of the data. Specifically, the training (optimization) of a DNN can be driven by minimizing the residual of the governing equations constructed by the DNN ansatz. This idea of physics-constrained learning is not new and was proposed back in the late ’90s in the context of solving classic differential equations [35, 36, 37]. However, limited by the NN techniques and computational power at that time, this seminal work didn’t have a big impact. Recently, this idea has been revived because of the recent advances in deep learning [32] combined with ever-increasing computational resources. Notably, the physics-informed neural network (PINN) proposed by Raissi et al. was used to solve a number of deterministic one-dimensional (1D) PDEs and two/three-dimensional (2D/3D) PDE-constrained inverse problems with a moderate amount of labeled data [33, 38, 39], e.g., noisy measurements of the velocity field. A similar approach is also applied to learn the constitutive relationship in a Darcy flow [40]. The PINN approach has been recently extended to assimilate multi-fidelity training data [41], and its UQ analyses have been explored based on arbitrary polynomial chaos [42] and adversarial inference [43]. Similar ideas of using physical constraints to regularize the DNN training have also been investigated in [44, 45, 46, 47]. In the aforementioned works, a moderate amount of labeled data either from simulations or experimental measurements are still needed for obtaining an approximation to the solution of the PDEs. In fact, if the initial and boundary conditions are well imposed thus the corresponding PDE problem is well-defined, in principle, the unique solution should be captured by the DNN via PDE-constrained learning without any labeled data. This is one of the main motivations of current work. Note that recently there have been several works on the concept of data-free DNNs, e.g., for solving a handful of computer vision problems [48], deterministic PDEs [49, 50, 33], high-dimensional stochastic partial differential equations (SPDE), and backward stochastic differential equations (BSDE) [51, 52, 53, 54, 55].

In the context of surrogate modeling, Nabian and Meidani [56] and Karumuri et al. [57] applied the PDE-constrained fully-connected neural network (FC-NN) for uncertainty propagation in steady heat equations. Zhu et al. [58] proposed a PDE-constrained, label-free DNN surrogate model for UQ in an elliptic PDE using both the FC-NN and convolutional neural networks (CNN). Moreover, both the deterministic and probabilistic formulations of physics-constrained learning are devised and studied. Their results have shown a significant potential of using the physics-constrained DNN for surrogate modeling, where no labeled data are required during the training. Nonetheless, the success has only been demonstrated in a few model problems governed by linear elliptic PDEs, e.g., diffusion equations. Thus it remains unclear if the physics-constrained learning can be used to handle realistic fluid systems governed by the Navier–Stokes equations in a parametric setting. Significant work is still needed to further explore the real-world problems for broad impacts.

The objective of this paper is to develop a physics-constrained, data-free DNN for surrogate modeling of incompressible flows. A structured FC-NN is devised to approximate the solutions of the parametric Navier–Stokes equations, where the initial/boundary conditions are enforced instead of being penalized together during training in previous works [33]. In addition, contrary to the previous data-driven surrogates, the training of our DNN is solely driven by minimizing the residuals of the governing PDEs (i.e., conservation laws), where no expensive CFD simulation data is used. To the best of our knowledge, the current work is the first attempt to build a parametric DNN surrogate model for fluid simulations without using any simulation data. The effectiveness and merits of the proposed method are demonstrated by investigating a number of internal flows relevant to cardiovascular applications. The rest of the paper is organized as follows. The framework of structured FC-NN surrogate based on the physics-constrained label-free training is introduced in Section 2. Numerical results of surrogate modeling and uncertainty propagation on several vascular flows are presented in Section 3. The performance of soft and hard boundary enforcement approaches, different adaptive activation functions, and data-free/data-driven learning strategies are discussed in Section 4. Finally, conclusion is drawn in Section 5.

2 Methodology

2.1 Overview

Most low-speed flows, e.g., blood flows in large or medium sized vessels, can be described by the incompressible Navier–Stokes equations given as:

[TABLE]

where $t$ and $\mathbf{x}$ are time and space coordinates, respectively; $\Omega_{f,t}\triangleq\Omega_{f}\times[0,T]$; $\boldsymbol{\theta}$ is a $d$-dimensional parameter vector, including input and/or operational parameters such as fluid properties, inlets/outlets, and geometry of the domain; both velocity $\mathbf{u}(t,\mathbf{x},\boldsymbol{\theta})$ and pressure $p(t,\mathbf{x},\boldsymbol{\theta})$ are functions of time, space, variable parameters; $\rho$ and $\nu$ represent density and viscosity of the fluid, respectively; $\mathbf{b}_{f}$ is the body force; $\Omega_{f}\subset\mathbb{R}^{3}$ denotes the fluid domain. The solutions of velocity $\mathbf{u}$ and pressure $p$ can be uniquely determined when suitable initial and boundary conditions are prescribed,

[TABLE]

where both $\mathcal{I}$ and $\mathcal{B}$ are general differential operators that define the initial and boundary conditions, respectively; $\partial\Omega_{f}$ denotes the boundary region. When a set of parameters $\boldsymbol{\theta}$ is given, the flow field, i.e., $\mathbf{u}(t,\mathbf{x})$ and $p(t,\mathbf{x})$, can be solved numerically by discretizing the Eqs. 1 and 2 using FD/FV/FE methods. However, this process involves mesh generation and iteratively solving large linear/nonlinear systems, which is usually time-consuming. Therefore, propagating the parameter uncertainty or inferring the unknown parameters through the FD/FV/FE solver becomes intractable when it comes to parametric problems, e.g., some parameters of $\boldsymbol{\theta}$ are uncertain or unknown. Solving varying-geometry problems is especially challenging since any change of the geometry requires regeneration of the computational meshes.

To enable fast predictions in terms of UQ and optimization applications, a deep neural network (DNN) architecture is built to approximate the solutions of the Navier–Stokes equations in a parametric setting. The DNN-based surrogate is expected to provide a rapid online prediction of the flow field with any given set of parameters $\boldsymbol{\theta}$ after the offline training. A schematic diagram of the proposed framework is shown in Fig. 1. A FC-NN is devised with the input layer composed of time $t$, spatial coordinates $\mathbf{x}$, and variable parameters $\boldsymbol{\theta}$. The raw outputs of the FC-NN (i.e., $\mathbf{u}$ and $p$) are used to construct the state variables (i.e., velocity $\mathbf{\hat{u}}$ and pressure $\hat{p}$) together with contribution from the particular solution that encodes initial/boundary conditions. The FC-NN is trained by minimizing residuals of the Navier–Stokes equations and no data from CFD simulations are needed. Therefore, the DNN predictions are expected to conform to the conservation laws of fluid flows and satisfy the specified initial/boundary conditions. Note that the Navier–Stokes equations will not be solved with any numerical discretization. The details of the physics-constrained training and boundary condition enforcement will be presented in the following subsections.

2.2 Deep Neural Network and Physics-Constrained Training

Neural networks (NN) are a set of algorithms, inspired by the biological neural networks in brains, for classification and regression tasks. There are various types of NNs with different neuron connection forms and architectures, e.g., fully-connected neural networks (FC-NN), convolutional neural networks (CNN), and recurrent neural networks (RNN). In this work, the feedforward FC-NN is considered, where the neurons of adjacent layers are fully connected and outputs of each layer are fed forward as the inputs to the next layer. A FC-NN defines a mapping from the input layer $\mathbf{z}_{0}\in\mathbb{R}^{n_{0}}$ to the output $\mathbf{z}_{L}\in\mathbb{R}^{n_{L}}$. The layers between the input and output layers are called hidden layers $\mathbf{z}_{l}$, where $l=1,\ldots,L-1$. By convention, a neural network with more than one hidden layer is called a “deep” NN. Mathematically, two adjacent layers are connected as,

[TABLE]

where $\mathbf{W}_{l}\in\mathbb{R}^{n_{l-1}\times n_{l}}$ and $\mathbf{b}_{l}\in\mathbb{R}^{n_{l}}$ are the weight matrix and bias vector; the subscript $l$ denotes the index of the layer; $\sigma_{l}(\cdot)$ is an activation function acting element-wise, for which a number of options can be chosen, e.g., sigmoids, rectified linear units (ReLU), and tanh functions. After training, the weights, bias, and activation function at each layer are determined, and the output prediction $\mathbf{z}_{L}$ (i.e., velocity and pressure) can be rapidly computed from any given input vector $\mathbf{z}_{0}$ (i.e., coordinates and parameters) based on the Eq. 3. Since this feedforward algorithm (Eq. 3) only involves a few matrix multiplications, the computational cost for evaluating the trained FC-NN can be neglected compared to that of a CFD simulation.

Traditionally, to build a surrogate model for the CFD simulation of the solution $\mathbf{f}(t,\mathbf{x},\boldsymbol{\theta})$, one can simply consider a black-box surrogate, e.g., FC-NN, or CNN [18], as $\mathbf{z}_{L}(t,\mathbf{x},\boldsymbol{\theta};\mathbf{W},\mathbf{b})$, i.e.,

[TABLE]

where $\mathbf{f}$ is the solution vector as $\mathbf{f}=\begin{bmatrix}\mathbf{u}&p\end{bmatrix}^{\intercal}$, including velocity $\mathbf{u}$ and pressure $p$; $\mathbf{W}$ and $\mathbf{b}$ denote the weights and biases of the entire network. Generally, training of a DNN is purely data-driven, and it consists of finding a set of (sub)optimal DNN parameters ($\mathbf{W},\mathbf{b}$) such that the mismatch between the training data $\mathbf{f}^{d}$ and the DNN predictions $\mathbf{f}^{p}$ is locally minimized. That is, one can formulate an optimization problem as,

[TABLE]

where the loss function $\mathcal{L}_{data}(\mathbf{W},\mathbf{b})$ is named as “data-based loss”, and $\left\lVert\cdot\right\rVert_{\Omega_{f,t}}$ is L2 norm over $\Omega_{f,t}$; $\mathbf{W}^{*},\mathbf{b}^{*}$ denote a set of (sub)optimal NN weights and biases obtained from the optimization.

However, as discussed above, this black-box surrogate modeling requires enormous training data $\mathbf{f}^{d}$, which is too expensive to obtain from a large number of CFD simulations. Instead, following previous PINN framework [33], we consider leveraging the governing PDEs in the loss function by minimizing the violation of the solution $\mathbf{z}_{L}$ in terms of the known governing PDEs for fluid dynamics over a domain of interests without the needs of solving these equations for each parameter with traditional numerical methods. Specifically, only the residuals of the Navier–Stokes equations are computed based on the FC-NN predictions and it corresponds to a constrained optimization as follows,

[TABLE]

[TABLE]

where the loss function $\mathcal{L}_{phy}(\mathbf{W},\mathbf{b})$ here is named as “physics-based loss”; the superscript $p$ indicates that the quantities are predicted by the DNN. To construct the PDE residuals in the loss function, several first and/or second derivative terms of $\mathbf{u}^{p}$ and $p^{p}$ with respect to time $t$ and space coordinates $\mathbf{x}$ are required, which can be computed based on automatic differentiation (AD) [59]. AD is an accurate and efficient way to calculate derivatives in a computational graph, which has started to gain increasing attention in the machine learning community. The general idea of AD is to use the chain rule to back-propagate derivatives from the output layer to the inputs as the connection between each layer of a NN is analytically defined. Compared to numerical differentiation techniques, derivatives calculated from AD are much more accurate since they do not suffer from truncation or round-off errors. Most modern deep learning frameworks such as PyTorch [60], TensorFlow [61], and Theano [62] have the AD implemented. To solve the optimization problem defined in Eq. 6, stochastic gradient descent (SGD) algorithms are used, which are known to be a stochastic approximation of the gradient descent (GD) optimization. In SGD, only a subset of points are randomly sampled from the input space to calculate the direction of the gradient at each iteration. The SGD algorithms are known to work very well to escape bad local minima in the neural network training [63] under one point convexity property. Although the global minimum cannot be guaranteed for a non-convex optimization problem as defined in Eq. 6, an empirically good local minimum is usually found based on the SGD algorithms.

2.3 Boundary Condition Enforcement

If the physics-based loss $\mathcal{L}_{phy}$ becomes identically zero, the DNN predictions of velocity $\mathbf{u}^{p}$ and pressure $p^{p}$ will exactly satisfy the Navier–Stokes equations (Eq. 1). Therefore, penalizing the PDE residuals can regularize the data-driven DNN solutions to be more physical. This idea is known as the physics-informed, weakly-supervised deep learning [33, 38, 45]. To make the problem well-posed, proper initial and boundary conditions (IC/BC) are needed and imposed as constraints (Eq. 7) which are often treated in a “soft” manner by modifying the original loss function with penalty terms [64, 33]. For example, the IC/BC can be imposed in a “soft” way by modifying Eq. 6 as,

[TABLE]

where $\lambda_{i}$ and $\lambda_{b}$ are penalty coefficients. However, the soft IC/BC enforcement methods have several major drawbacks: (1) there is no quantitative guarantee on how accurate the IC/BC being imposed and thus the solution could be unsatisfactory; (2) the optimization performance can depend on the relative importance of each term, but how to assign weight for each term can be difficult. Alternatively, we can impose the IC/BC in a “hard” manner, where a particular solution that solely satisfies the initial/boundary condition is added. Hence, the constraints on IC/BC are automatically fulfilled. A mixed enforcement on IC/BC is proposed in this work, where the Neumann and Dirichlet boundary conditions (BC) are treated separately: the Neumann BC are formulated into the equation loss, i.e., in a soft manner, while the IC and Dirichlet BC are encoded in a hard manner by constructing the DNN ansatz $\hat{\mathbf{u}}^{p}$ and $\hat{p}^{p}$ with a particular solution as follows,

[TABLE]

where $\mathbf{u}_{particular}$ is a particular solution that just satisfies IC and BC: $\mathbf{u}_{particular}(\mathbf{x},0)=\mathbf{u}_{0}(\mathbf{x})$, $\mathbf{p}_{particular}(\mathbf{x},0)=p_{0}(\mathbf{x})$ and $\mathbf{u}_{particular}(\mathbf{x},t)|_{\mathbf{x}\in\partial\Omega_{f}}=\mathbf{u}^{b}(\mathbf{x})$, $\mathbf{p}_{particular}(\mathbf{x},t)|_{\mathbf{x}\in\partial\Omega_{f}}=p^{b}(\mathbf{x})$; $D(t,\mathbf{x},\boldsymbol{\theta})$ is a globally defined smooth function from internal points to the “boundary” in $\Omega_{f,t}$, i.e., a space–time sense. That is, $D$ is zero on the boundary $\partial\Omega_{f}\times[0,T]$ and $\Omega_{f}\times\{0\}$ while increases away from the boundary. For those problems where the IC/BC and the geometry of the domain in $\Omega_{f,t}$ is simple, the function $D$ and particular solution can be written analytically. However, if the geometry is too complex to have an analytic form, e.g., a patient-specific artery, the function $D$ can be pre-trained by a low-capacity NN in the way proposed in [50]. Finally, the constrained optimization problem in Eq. 6 can be reformulated as an unconstrained one, as shown in Eq. 10.

[TABLE]

3 Numerical Results

A number of 2D vascular flows with idealized geometries, including circular pipe flows, stenotic flows, and aneurysmal flows, are studied to evaluate the performance of the proposed data-free DL surrogate model in Eq. 10. Forward propagation of the uncertainties in the fluid properties and domain geometry to the flow QoIs (e.g., velocity, pressure, and wall shear stress) are investigated through the DNN surrogate model. In this study, only steady-state solutions are considered for proof-of-concept, thus the constraint of initial condition can be neglected in Eq. 9.

A composite FC-NN architecture is devised for the surrogate, which is composed of three sub-DNNs with an identical structure of 3 hidden layers with 20 neurons per layer. The Swish activation function [65] with fixed hyperparameters is employed in each layer except the last one, where a linear activation function is used. The three sub-DNNs share the same input layer and separately predict three scalar state variables, i.e., velocity $u,v$, and pressure $p$. All three sub-DNNs are trained simultaneously with a unified physics-based loss function. To solve the unconstrained optimization problem defined in Eq. 10, we used the Adam optimizer [66], a robust variant of the SGD method, where the learning rate is adaptively changed based on the estimates of the moments. The initial learning rate and mini-batch size is set as $1\times 10^{-3}$ and 50, respectively. Because of the adaptivity feature of the Adam optimizer, the hyperparameters in the training are robust to some extent and require little tuning. Properly initializing the DNN parameters is also important. After comparing several widely-used initialization schemes [67, 68], we chose the He’s normal initializer [68], where initial weights are drawn from a truncated normal distribution. In general, a good choice of the DNN architecture, including the number of layers, number of neurons per layer, activation function, and initialization schemes, is important to the learning performance but is still determined by trial and error. A rule of thumb is to achieve a favorable performance using the simplest network structure, which enables fast training and better generalizability. In this work, a handful of DNN architectures with different “depths” are investigated, and the “shallowest” one with a satisfactory learning performance is adopted. To demonstrate the robustness of the physics-constrained learning, the architecture and hyperparameters remain the same for all the cases throughout the paper. Note that a comprehensive parameter study and architecture optimization of the DNN is out of the scope of the current work.

The composite FC-NN is implemented in the PyTorch platform [60]. As discussed in Section 2, only the collocation points are required and they are uniformly sampled in the spatial $\mathbf{x}$ and parameter $\boldsymbol{\theta}$ spaces. Alternatively, one can choose a space-filling Latin hypercube sampling in the $\Omega_{f}$ [33]. In this work, the PDE residuals are extensively evaluated on a large number of collocation points to ensure learning quality. For all test cases, the training of about $10^{6}$ SGD iterations are performed on an NVIDIA GeForce GTX 1080Ti Graphics Processing Unit (GPU) card, and the cost is approximately 3.5 hours. Note that the offline training cost can be potentially reduced by optimizing the DNN architecture. To validate the prediction performance of the trained DNN surrogates, corresponding CFD simulations are also conducted using an open-source FV-based CFD solver, OpenFOAM [69]. Mesh convergence study is performed to ensure the solution accuracy. The code and datasets for this work will become available at https://github.com/Jianxun-Wang/LabelFree-DNN-Surrogate upon publication.

3.1 Circular Pipe Flow

Flow in a pipe/tube is very common in physiological systems, e.g., blood in arteries or air flow in trachea. The pipe flow is often driven by the pressure difference between the two ends of a tube, or by the body force of gravity. In a cardiovascular system, the former one is more dominant since the blood flow is mainly governed by pressure drop due to the heart pumping. In general, simulating the fluid dynamics in a tube requires solving the full Navier–Stokes equations numerically, but if the tube is straight and has a constant circular cross section, analytical solution of the fully-developed steady-state flow is available, which is an ideal benchmark to validate the performance of the proposed method. As a result, we first study the flow in a 2D circular pipe (also known as the Poiseuille flow).

In this case, the pressure inlet and outlet are used to drive the flow since we only focus on the fully-developed regime, and no-slip wall boundary is prescribed on the tube walls. The boundary conditions are encoded into the surrogate model by constructing the DNN ansatz $\hat{u}^{p},\hat{v}^{p}$, and $\hat{p}^{p}$ based on Eq. 9. The no-slip condition of velocity on the wall can be imposed by designing the $\hat{u}^{p},\hat{v}^{p}$ as,

[TABLE]

where $y$ is the radial distance, $d=0.1$ is the diameter of the tube, the raw DNN output is denoted by $u^{p}$. The pressure inlet $p_{in}=0.1$ and outlet $p_{out}=0$ are imposed by designing the $\hat{p}^{p}$ as,

[TABLE]

where $x_{in}$ and $x_{out}$ are coordinates of the two ends of the tube, and the raw DNN output is denoted by $p^{p}$. All three sub-DNNs are trained to capture the spatial flow fields with parameter variation in the fluid viscosity $\nu$. Input (collocation) points in the parameter space of $\nu$ are uniformly sampled in the range $10^{-4}\leq\nu\leq 1.9\times 10^{-3}$, where the corresponding Reynolds number (Re) are moderate (Re $<300$). After training, both velocity and pressure fields can be obtained immediately by evaluating the trained DNN with any given input $\nu$ and a spatial grid on $\mathbf{x}$. Hence, the DNN surrogate can be utilized to propagate the uncertainty in viscosity $\nu$ based on the Monte Carlo (MC) simulation, where a large number of samples are drawn from the $\nu$ distribution and propagated to the QoIs via the DNN surrogate. In the following test cases, 500 MC samples are used to compute desired statistics and distributions.

The DNN surrogate results (shown in Figure 2) are compared against the analytical solution, which is given by

[TABLE]

where $y$ denotes the spanwise coordinate and $\Delta p$ is the pressure difference. It can be observed from Fig. 2a that the DNN-predicted velocity profiles (red dashed lines) of four different $\nu$ samples almost exactly agree with the analytical solutions (blue solid lines), where the Reynolds numbers (Re) of the four cases are 283, 121, 33, and 3, respectively. Actually, the trained DNN is able to accurately predict the pipe flow field with any given viscosity, where the Reynolds number is moderate. Fig. 2b shows the uncertainty of the center velocity $u_{c}$ propagated from a normally distributed $\nu$ with mean of $\overline{\nu}=1\times 10^{-3}$ and variance of $\sigma_{\nu}=2\times 10^{-4}$. It can be seen that the DNN-predicted probability density function (PDF) is almost identical to that of the analytical solutions, demonstrating the excellent performance of the physics-constrained learning for uncertainty propagation.

3.2 Blood Flow with Standardized Geometry

Two types of canonical vascular flows, stenotic flow and aneurysmal flow, with standardized vessel geometries are studied. A stenotic flow refers to the flow through a vessel where there is narrowing and re-expansion of the vessel wall. This local restriction of the vessel is related to many cardiovascular diseases, e.g., arteriosclerosis, stroke, and heart attack [70]. The vascular flow within an aneurysm, which is the expansion of an artery due to the weakness of vessel walls, is called aneurysmal flow. The rupture of an aneurysm may cause life-threatening conditions, e.g., subarachnoid hemorrhage (SAH) due to cerebral aneurysm rupture [71], and investigation of the hemodynamics can improve the diagnosis and fundamental understanding of aneurysm progression and rupture [72].

Whereas realistic vascular geometries are usually irregular and complex, including sites of curvature, bifurcation and junctions, idealized stenosis and aneurysm models are studied here for proof-of-concept. Namely, both the stenotic and aneurysmal vessels are idealized as an axisymmetric tube with a varying cross-section radius, which is parameterized by the following function,

[TABLE]

where $R_{0}$ is the radius at the inlet, which is set as 0.05, and the sign of $A$ determines if the vessel is stenotic or aneurysmal. Namely, a positive and negative sign correspond to the stenosis and aneurysm, respectively. Three control parameters $A$, $\mu$, and $\sigma$ define the shape of the stenotic (aneurysmal) vessel. The scale parameter $A$ controls the curvature along the tube, and a larger $|A|$ leads to narrower stenosis (broader aneurysm). Parameter $\mu$ defines the streamwise location of the minimum (maximum) radius of the stenosis (aneurysm), and $\sigma$ affects the steepness of the geometric variation. In this study, the latter two parameters $\mu$ and $\sigma$ are fixed as $0.5$ and $0.1$, respectively. Only $A$ is considered as a variable parameter to control the degree of the stenosis (aneurysm).

Similar to the pipe flow, the pressure inlet/outlet and no-slip wall boundary conditions are prescribed for both stenosis and aneurysm cases. The boundary-encoded sub-DNNs are constructed to learn the parametric flow solutions, where the wall BC is imposed using the geometric function $D(x,y)=R(x)^{2}-y^{2}$. Contrary to the case studied above, where only the viscosity variation is considered, the DNN surrogate is also trained to capture the varying stenosis (aneurysm) geometry, which is known challenging for mesh-based CFD simulations. Specifically, the solutions of varying viscosity are learned for fixed vascular geometries ($A=5\times 10^{-3}$ for stenosis and $A=-5\times 10^{-3}$ for aneurysm) in the first place, and then the performance for capturing geometry variations is examined at a fixed viscosity ($\nu=1\times 10^{-3}$). Moreover, the uncertainties from the flow viscosity and vessel geometry are propagated to the QoIs through the trained DNN surrogate using MC sampling, and the results are validated by CFD-based MC simulations.

3.2.1 Flow in Idealized Stenosis

The DNN is trained to parameterize the solutions of stenotic flows with varying viscosity, where collocation points are sampled in $\nu$ space within the range of $[5\times 10^{-4},1\times 10^{-2}]$ for physics-based training. \floatsetup[figure]style=plain,subcapbesideposition=top

Figure 3 shows the DNN-predicted flow fields of three different viscosity samples, i.e., $\nu=6.4\times 10^{-4}$, $1.85\times 10^{-3}$, and $2.14\times 10^{-3}$, at moderate Reynolds numbers. The corresponding CFD simulations are performed for comparison. It can be seen that the flow patterns of different $\nu$ are similar, where the fluid is accelerated streamwisely through the converging region and slows down passing the diverging part of the tube. However, the velocity magnitude reduces as the viscosity increases (left to right columns). As shown in Figs 3a and 3b, the DNN-predicted velocity contours of streamwise and spanwise components ($u$ and $v$) agree with the CFD solutions very well, though the magnitude in the case with the smallest viscosity ($\nu=6.4\times 10^{-4}$) is slightly underestimated. Moreover, the nonlinear pressure drops can be accurately captured by the DNN surrogate as the profiles of centerline pressure from the DNN and CFD are almost identical (Fig. 3c).

To learn the flow solutions with varying geometry, the DNN is trained on uniformly sampled points within the range of $0\leq A\leq 1\times 10^{-2}$.

We compare the flow fields predicted by the trained DNN against the CFD benchmarks, and the results of three different samples of stenosis geometries, $A=2\times{10^{-3}}$, $4\times{10^{-3}}$, and $7\times{10^{-3}}$, are shown in Fig. 4. From left to right, the degree of stenosis increases and thus the total flow rate is reduced due to the increased resistance (Fig. 4a and Fig. 4b). The pressure drop becomes more nonlinear as the stenotic vessel turns to be narrower (Fig. 4c). We can see that the DNN predictions can capture these flow features and agree with the CFD benchmarks well. Admittedly, the streamwise velocity of the flow with the narrowest stenotic vessel is slightly underestimated, and the centerline pressure profiles predicted by the DNN and CFD has a small discrepancy. This might be because the increased nonlinearity due to the steep geometric variation poses a challenge on the learning.

After the training, the DNN surrogates are used to rapidly propagate uncertainties in viscosity and vessel geometry, and the effects on QoIs are investigated. Specifically, 500 MC samples are drawn from a normally distributed viscosity $\nu$ and geometric parameter $A$, which are propagated to the center velocities $u_{c}$ (at $x=0.5,y=0.0$) through the DNN surrogates. The probability distributions of $u_{c}$ due to uncertain viscosity and vessel geometry are shown in Figs. 5a and  5b, respectively, where the propagated results through the CFD solver are also plotted for comparison.

It shows that the propagated distributions of $u_{c}$ are non-Gaussian in both cases, which is due to the strong nonlinearity of the Navier–Stokes operator. As expected, the DNN-propagated uncertainties present a good agreement with the CFD-based benchmarks, especially in the case of viscosity uncertainty propagation (Figs. 5a), where the two PDF curves are almost overlapped with each other. As for the geometry uncertainty propagation, although the overall feature of the PDF is captured, the peak density is slightly underpredicted by the DNN. The reason behind this could be that the DNN surrogate tends to underestimate the velocity magnitude in particular for a steep geometry variation.

3.2.2 Flow in Idealized Aneurysm

We first learn the aneurysmal flows with varying viscosity, where the geometry of the aneurysm is fixed ($A=-5\times{10^{-3}}$).

Training is conducted by sampling the viscosity points ranging from $5\times 10^{-4}$ to $1\times 10^{-2}$. Figure 6 shows the DNN-predicted flow fields of three viscosity samples, i.e, $\nu=6.4\times 10^{-4},1.85\times{10^{-3}}$, $2.14\times 10^{-3}$, where the CFD benchmarks are plotted for comparison. In addition to the flow velocity and pressure, here we also investigate the wall shear stress (WSS) $\tau$, which has been demonstrated as a critical factor affecting the aneurysm initialization, progression, and rupture [73]. It can be seen from Fig. 6 that the DNN-predicted flow QoIs are in a good agreement with the CFD solutions. The decrease of the flow velocity and the minimum WSS of the vessel are accurately captured by the DNN surrogate.

Figure 7 shows the performance for learning the geometry-varying solutions of the aneurysmal flow, where the training is conducted by sampling the geometric parameter $A$, ranging from $-2\times 10^{-2}$ to [math]. The predicted flow fields of three different geometry samples are presented, where the size of the aneurysm increases from left to right. The flow decelerates through the expanded region of the vessel and the velocity at the center of the aneurysm is significantly reduced, in particular when the aneurysm becomes larger. It is observed from the contour comparisons (Figs. 7a and 7b), the DNN predictions agree with the CFD solutions pretty well. As for the WSS profile, its shape and magnitude vary as the geometry changes, which can be accurately captured by the DNN surrogate (Figs. 7c).

The uncertainty propagation using the trained DNN surrogates is then conducted.

Uncertainties in viscosity $\nu$ and geometric parameter $A$ with Gaussian distributions are considered, and the QoIs are center velocity $u_{c}$ and the minimum WSS $\tau_{c}$, which are important for aneurysmal flows. As shown in Fig. 6c, the WSS remains invariant due to viscosity perturbation with a fixed vessel geometry, hence the distribution of $u_{c}$ is studied for viscosity uncertainty propagation (shown in Fig. 8a). The PDF obtained by the DNN surrogate almost coincides with the CFD-based benchmark. As for the geometry uncertainty propagation, the center velocity $u_{c}$ is not an interesting quantity since it almost remains the same as geometry changes. Instead, we investigate the propagated uncertainty in the minimum WSS $\tau_{c}$, which is sensitive to geometry variation. It is observed from Fig. 8b that the probabilistic distribution of $\tau_{c}$ propagated by the DNN is in a favorable agreement with the CFD-based benchmark, though the peak of the density is slightly underestimated.

3.2.3 Summary of Training and Prediction Performance

The training and testing performance on the vascular flow cases are summarized in Table 1, where the training loss is defined as the sum of $L_{2}$ norms of momentum and continuity equation residuals (see eq. 10) and test error is defined as the normalized $L_{2}$ difference between DNN and CFD results, shown as,

[TABLE]

where pressure drop $\Delta p=0.1$ is used to normalize the errors.

It can be seen from Table 1 that both the training loss and test errors are reasonably small for all cases after training, and the test error in $u$ is more dominant than the others. Among the four scenarios, the geometric variation of the stenosis is the most challenging to learn since the training loss (i.e., equation residual) remains relatively large and the test errors are one-order bigger than the other cases. This is consistent with the notable discrepancy in the DNN-predicted PDF observed in Fig. 5b.

All the cases presented above are trained on $N_{p}=20$ parameter collocation points. Namely, the equation residuals are minimized on $20$ collocation points uniformly sampled from the parameter space. It is necessary to check if the size $N_{p}$ of parameter collocation points is sufficiently large for the training. Therefore, we conduct a parameter study using different amounts of parameter collocation points for training. The total test errors (sum of test errors in $u,v$, and $p$) against different numbers of training collocation points for all four cases are presented in Fig. 9.

As expected, the errors generally decrease as the number of training collocation points increases. However, the error decreasing rate is quite mild. The test errors remain approximately unchanged when $N_{p}>20$ for all cases while the training cost will increase as more collocation points are used, which justifies the sufficiency of the total collocation points used in this work. Detailed results of training cost and testing errors for all cases are summarized by Tables 3 and 4 in Appendix A.

3.2.4 Computational Cost of Uncertainty Propagation

Finally, we briefly discuss the computational cost of the uncertainty propagation tasks presented above, which is determined by the forward model-evaluation time and number of MC samples. For both stenotic and aneurysmal flows, a structured mesh with $10^{4}$ quadrilateral grids is used for the CFD model. The computational cost of each forward solution of the fully-converged CFD simulation is around 40 CPU seconds, and propagating 500 MC samples in test cases takes about 20,000 CPU seconds on a single CPU core. However, the online evaluation of the trained DNN surrogate is very fast, and the cost of each forward DNN evaluation is less than $2\times 10^{-2}$ CPU seconds. Therefore, the propagation of 500 MC samples by the DNN surrogate only takes about 10 CPU seconds, which gains 2000-times speedup over the CFD-based uncertainty propagation. This advantage can be considerable when a larger number of forward evaluations are needed or more complicated fluid systems are considered. Moreover, the FV/FE-based CFD simulations require mesh generation, which is often a manually-cumbersome and labor-intensive process, in particular for the flows with complex geometries and moving boundaries, e.g, patient-specific cardiovascular simulations. Therefore, the mesh-free feature of the proposed method shows great promises for the many-query analysis in these systems.

4 Discussion

4.1 Pitfall on Using Soft Boundary Enforcement

The initial and boundary conditions (IC/BC) can be imposed in the physics-constrained learning either as a soft or hard constraint. When no labeled data is used in training, a properly enforced IC/BC is crucial to ensure the uniqueness of the learned PDE solutions. Although we have demonstrated the effectiveness of the hard enforcement approach in Sec. 3, it is still interesting to investigate the performance of soft enforcement method [33] in the purely PDE-driven training. Hence, all the test flows are studied again using physics-constrained learning, where BCs are imposed in a soft manner (as Eq. 16). Namely, the BCs are formulated as a boundary loss component $\mathcal{L}_{B}$, which is incorporated into the physics-based loss function $\mathcal{L}_{phy}$ as,

[TABLE]

where $\lambda$ is the penalty coefficient. For the circular pipe flow case, both the hard and soft BC constraints can lead to excellent learning and prediction performance (results not shown). However, when the radius varies along the tube as in the stenosis and aneurysm cases, the DNN with the soft BC enforcement does not perform well as the no-slip BC of the vessel wall is poorly imposed especially near the bottleneck in Fig. 10. Consequently, the solution to the flow field becomes inaccurate. For example, Fig. 10 shows the results for learning a stenotic flow ($A=5\times 10^{-3}$) using the soft constraint with different $\lambda$ values, where the result with hard boundary constraints and CFD benchmark are plotted as well for comparison.

In contrast to the result with hard constraints (Fig. 10b), both the flow patterns and magnitudes predicted with the soft boundary enforcement with different $\lambda$ (Fig. 10c–f) are completely wrong compared to the CFD benchmark (Fig. 10a). Moreover, the nonlinear behavior of the centerline pressure profile cannot be accurately captured (results not shown).

The poor performance reflects the major drawbacks of the soft constraints in the physics-driven training as mentioned above. First, unique PDE solutions are determined by the IC/BC, which cannot be guaranteed by simply penalizing the boundary loss. Moreover, it is difficult to assign the relative weight ($\lambda$) for different components in the loss function, and there could be “competing” effects between the equation loss and boundary loss, which makes the optimization difficult to converge. As shown in Table 2, by assigning a large weight for the boundary loss ($\lambda=1000$), the boundary condition can be well prescribed but the PDE residual remains a large value and cannot be further reduced. On the other hand, when the weight for the boundary is small ($\lambda=1$), the loss of x-momentum equation can be reduced to $O(10^{-4})$ but the BC fails to be imposed. It is important to note that none of these four $\lambda$ leads to a physical stenotic flow pattern as shown in Fig. 10.

4.2 Role of Activation Function in Physics-Constrained Learning

The performance of DNN training is affected by the activation function to a large degree. The widely used activation functions includes ReLU, Sigmoid, Tanh, etc. [74]. However, these activation functions are not guaranteed to be optimized in terms of the convergence rate and accuracy. Recent studies [75, 76] proposed to train an adaptive activation function as well as the neural network weights to achieve better convergence property. Notably, Ramachandran et al. [75] introduced an adaptive activation function called Swish, which is defined as $x\cdot\textrm{Sigmoid}(\beta x)$ and $\beta$ is a trainable parameter. Jagtap and Karniadakis [76] presented a new adaptive function defined as $\textrm{Tanh}(nax)$, where $a$ is an adaptive parameter to be learned and $n$ is a scale factor that potentially speeds up the convergence.

In current work, the training process uses a Swish activation function with fixed $\beta=1$. It’s not clear how the adaptivity of activation function can affect the convergence rate and accuracy. Furthermore, as we discussed and highlighted the necessity of the hard BC enforcement in the physics-constrained data-free learning, it is also interesting to compare the relative importance of boundary condition enforcement and adaptive activation function on model performance. Therefore, we test the effects of different activation functions with/without hard boundary constraints in the stenosis case ($A=5\times 10^{-3}$, $\nu=1\times 10^{-3}$), and the resulting learning curves are shown in Fig. 11, where panels (a) and (b) show the convergence histories of different activation functions with and without a hard BC enforcement, respectively.

The legend “Swish-$\beta$” refers to the Swish function with adaptive $\beta$, while the “Swish” means the Swish function with a fixed $\beta=1$. The results of adaptive Tanh activation function with different $n$ [76] are plotted as well. It can be seen in Fig. 11a that the convergence rate of using adaptive Tanh functions is faster than that of the Swish activation functions within the first 400 epochs, but the training loss of all cases finally converge to the same order. Moreover, the convergence curves almost overlap with each other for the same type of activation function with different hyperparameters. When the BC is imposed in a soft manner (Fig. 11b), the solution can be quickly trapped in a bad local minimum as we discussed above and the final training loss is over one-order larger than that of the cases with hard BC constraints. Using different adaptive activation functions does not help to further decrease the loss and all the convergence curves are overlapped. Moreover, the convergence history of each trainable hyperparameter (i.e., $\beta$ in Swish function and $a$ in Tahn function) is monitored (see Fig. 12 in Appendix B), and we found that these hyperparameters are more likely to converge to the optimized values when the BCs are enforced in a hard manner. To sum up, the way of imposing BCs is found to be more critical than the adaptivity of activation function in our cases in terms of accuracy in the physics-constrained data-free learning.

4.3 Physics-Constrained Data-Free Learning vs. Traditional Data-Driven Learning

We have demonstrated that the solutions of parametric Navier-Stokes equations can be effectively learned by solely minimizing the PDE residuals without using any simulation data. To better evaluate the advantages and limitations of the physics-constrained learning, we also conducted a series of comparison studies between the physics-constrained data-free DNN and traditional data-driven DNN in terms of learning efficiency and prediction accuracy. Namely, all the vascular flow cases discussed above are learned again in a purely data-driven way, where the DNN architecture and hyperparameters remain the same and only the physics-based loss function (i.e., Eq. 10) is replaced with the data-based one (i.e., Eq. 5). A number of CFD simulations with different input parameters are conducted and the simulated velocity and pressure fields are collected as the training data. The detailed comparisons for stenotic and aneurysmal flows cases with varying viscosity and geometry are summarized by Tables 3 and 4 in Appendix A. In general, the prediction results from purely data-driven learning are slightly more accurate than those of physics-constrained learning, and the accuracy of both models improves as the parameter collocation points increase. However, data-driven learning requires additional offline CFD simulations and this computational overhead can quickly grow as more training points are sampled from the parameter space. In this paper, since the CFD cases considered here are not costly to simulate (e.g., each CFD simulation takes about 40 CPU seconds), the computational overhead due to the offline data generation process is not significant. Nonetheless, the advantage of the data-free feature in physics-constrained learning will become more notable when large-scale 3D flow problems are considered, where a single simulation run could be very expensive. It notes that when training data are ready for use, the cost of data-driven training process is approximately similar to that of the physics-constrained training (in the same order), though the latter is slightly slower due to the additional AD calculations for derivatives. To reduce the training cost in data-driven learning, one way is to reduce the spatial dimensionality by projecting the training data (i.e., velocity/pressure fields) onto the POD basis, and learning is performed on POD coefficients instead of spatial collocation points. We also conducted POD-based data-driven learning and found that the learning performance with the current shallow network structure is unsatisfactory (results are not shown here for conciseness).

Admittedly, the current form of physics-constrained DNN has its limitations, for example, the offline training process is still costly, the convergence cannot be guaranteed due to the non-convexity of DNN optimization, and scalability for high-dimensional complex problems is still challenging. The proposed PDE-constrained DNN is not expected to replace the classical CFD (numerical) solvers, which have been developed for decades. However, the development of PDE-constrained DNN for surrogate modeling shows strong promise. Particularly, the proposed method is mesh-free and thus does not require arduous mesh generation labor and intensive domain expertise in numerical modeling, which is suitable for, e.g., rapidly testing ideas in the design phase. We expect that the effectiveness of surrogate model based on physics-constrained DNN will be significantly promoted along with the rapid advances in improvement of DNN training efficiency, e.g., a recent study has suggested that a novel photonic chip has a potential to be used to train deep neural networks 10-million of times more efficiently than current CPUs/GPUs do [77].

5 Conclusion

Surrogate modeling of fluid flows governed by the Navier–Stokes equations is significant for uncertainty quantification, optimization design, and inverse analysis in many engineering systems. As a universal function approximator, DNN is becoming a popular approach for surrogate modeling. However, training of a DNN often requires large number of labeled data, which are usually not available for efficiently developing surrogates since each data point requires an expensive CFD simulation. This paper presented a novel DNN surrogate for fluid simulations without using any labeled data (i.e., CFD simulation data). Specifically, a structured DNN architecture is devised to approximate the solutions of the parametric Navier–Stokes equations, where the initial/boundary conditions are satisfied automatically. Instead of using any simulation data, the DNN is trained by solely minimizing the violation of the mass and momentum conservation laws for fluid flows. Compared to the previous works of physics-constrained learning, this paper focuses on modeling of fluid systems governed by parametric Navier–Stokes equations. The proposed methods were tested on three flow cases relevant to cardiovascular applications, i.e., circular pipe flow, stenotic flow, and aneurysmal flow. The DNNs with equation-based loss were trained to learn the flow fields with parameter variations in, e.g., viscosity and domain geometry. Uncertainties in these parameters are propagated through the trained DNN surrogate and the results are validated against the CFD benchmarks. The comparisons indicate the excellent agreement between the physics-constrained DNN surrogate models and CFD simulations. Without using any labeled data in training, the DNN is able to accurately parameterize the velocity/pressure solutions with varying viscosity and geometries, which can be used to efficiently propagate uncertainties with enormous MC samples. Moreover, the performances of using hard and soft IC/BC enforcement approaches are compared and the issues of soft constraints in physics-constrained are discussed. We also investigated the influence of state-of-art adaptive activation functions and compared the present labeled-data-free learning approach with traditional data-driven learning approach in terms of accuracy and efficiency. In summary, the results have demonstrated the merit of the proposed method and suggest a great promise in developing DNN for surrogate fluid models without the need for CFD simulation data.

Compliance with Ethical Standards

The authors declare that they have no conflict of interest.

Acknowledgement

LS gratefully acknowledge partial funding of graduate fellowship from China Scholarship Council (CSC) in this effort. JXW would acknowledge support from the Defense Advanced Research Projects Agency (DARPA) under the Physics of Artificial Intelligence (PAI) program (contract HR00111890034). The authors also would like to thank Dr. Nicholas Zabaras and Mr. Yinhao Zhu for their helpful discussions during this work.

Appendix A

The physics-constrained, data-free model and purely data-driven model are trained on a single NVIDIA 1080TI GPU card.

Appendix B

Figure 12 shows the convergence histories of trainable parameters of different adaptive activation functions. The comparison is also made between cases where the BCs are imposed in hard and soft manners. It can be found that the parameters of activation function are more likely to converge with hard BC enforcement (panels (a)–(e)) compared to the cases with soft BC enforcement (panels (f)–(j)).