Physics Informed Neural Networks for Modeling of 3D Flow-Thermal
  Problems with Sparse Domain Data
Successfully training Physics Informed Neural Networks (PINNs) for highly
nonlinear PDEs on complex 3D domains remains a challenging task. In this paper,
PINNs are employed to solve the 3D incompressible Navier-Stokes (NS) equations
at moderate to high Reynolds numbers for complex geometries. The presented
method utilizes very sparsely distributed solution data in the domain. A
detailed investigation on the effect of the amount of supplied data and the
PDE-based regularizers is presented. Additionally, a hybrid data-PINNs approach
is used to generate a surrogate model of a realistic flow-thermal electronics
design problem. This surrogate model provides near real-time sampling and was
found to outperform standard data-driven neural networks when tested on unseen
query points. The findings of the paper show how PINNs can be effective when
used in conjunction with sparse data for solving 3D nonlinear PDEs or for
surrogate modeling of design spaces governed by them.
Saakaar Bhatnagar, Andrew Comerford, Araz Banaeizadeh
# Physics Informed Neural Networks for Modeling of 3D Flow-Thermal Problems with Sparse Domain Data

###### Abstract

Successfully training Physics Informed Neural Networks (PINNs) for highly nonlinear PDEs on complex 3D domains remains a challenging task. In this paper, PINNs are employed to solve the 3D incompressible Navier-Stokes (NS) equations at moderate to high Reynolds numbers for complex geometries. The presented method utilizes very sparsely distributed solution data in the domain. A detailed investigation on the effect of the amount of supplied data and the PDE-based regularizers is presented. Additionally, a hybrid data-PINNs approach is used to generate a surrogate model of a realistic flow-thermal electronics design problem. This surrogate model provides near real-time sampling and was found to outperform standard data-driven neural networks when tested on unseen query points. The findings of the paper show how PINNs can be effective when used in conjunction with sparse data for solving 3D nonlinear PDEs or for surrogate modeling of design spaces governed by them.

**Keywords:** Physics Informed Neural Networks; Navier-Stokes Equations; Surrogate Modeling; Design Optimization

## 1 Introduction

Over the last few years, there has been significant growth in the popularity of machine learning algorithms to solve partial differential equations (PDE) or assist PDE solvers, such as computational fluid dynamics (CFD) solvers [1, 2]. A particular application where CFD solvers struggle, due to the computational cost, is iterative design optimization. This is the process of continually updating a design (e.g. an electronics assembly layout) and computing the solution (e.g. flow or thermal fields) to optimize the performance (e.g. constrain the temperatures or reduce the pressure drop). The challenge for CFD is the input-output relationship is one-to-one. Therefore, any changes to the input vector (e.g. geometric variations) need to be re-simulated, leading to high costs when iterating on different design scenarios [3]. Overall high-fidelity iterative design requires a prohibitive level of resources, both computationally and monetarily, and often leads to a sub-optimal outcome. The attraction of Machine Learning (ML) algorithms in these scenarios is the ability to rapidly find solutions for such problem setups that are challenging in conventional CFD, such as large design space explorations [4], turbulence model closure [5] or solving incomplete/ill-posed problems [6].

Conventional ML algorithms usually require large amounts of data to train. This represents a challenge when using ML in engineering applications such as CFD, since experimental data can be difficult and expensive to obtain and may suffer from measurement noise. Furthermore, in many engineering experiments, field data such as temperature and velocity fields can sometimes only be captured at specific locations, and it is difficult to get full field solution results from physical experiments. Research has turned to using simulation data for training ML models, but the computational cost of generating large amounts of data to train models is a major bottleneck.

Physics Informed Neural Networks (PINNs) [7] represent an advance in scientific machine learning that has the potential to solve many of the aforementioned issues. By adding the physics that governs the problem into the loss function, and optimizing the loss, it is possible to have the network learn the
solution of the problem represented by that equation in a data-free manner. PINNs can be used in cases where sporadic experimental field data is available [8, 9] to calculate the rest of the field variable and can be used to solve problems with incomplete or missing physics [10, 11].

Another application area, in which PINNs could be very beneficial is machine learning-based surrogate modeling. Although a relatively new field, several ML architectures and methods have been utilized in the literature. These include: Proper Orthogonal Decomposition (POD) [12], Gappy POD [13] and Manifold Learning [14]. More recently, increased attention has been given to statistical methods like Gaussian processes and neural networks that incorporate Machine Learning (ML) to create surrogate models. Bhatnagar et al. [15] used a CNN architecture to predict aerodynamic flow fields over airfoils and created a surrogate model that generalized between flow conditions and airfoil geometries. Guo et al. [16] also used a Convolutional Neural Network (CNN) architecture to predict steady flows over automotive vehicles. Lee and You [17] used Generative Adversarial Networks (GANs) coupled with physical laws to predict unsteady flow around a cylinder, demonstrating the benefits of using embedded physics. Raissi and Karniadakis [18] use Gaussian processes to model and identify several complex PDEs.

Several of the aforementioned studies used purely data-driven models and required the creation of large amounts of training data to generate accurate and generalizable models. PINNs have the capability to greatly reduce these data generation costs, and it has been shown that training surrogates using the physics embedded in the loss function greatly improves predictive accuracy, across a wide range of applications [17, 19, 20, 21].

However, there is currently a lack of research articles applying PINNs to 3-dimensional (3D) problems, particularly for highly nonlinear PDEs like the Navier-Stokes equations. These problems are challenging for PINNs due to a variety of reasons that are discussed later in this paper. Yet, these problems are the most lucrative to solve, as most industrial applications of CFD are done in 3D. This paper provides results that aim to address this gap, by solving several problems with realistic physical parameters, over complex geometries in a data-assisted manner, using very sparse domain data. Further, this paper solves a realistic flow-thermal design optimization problem using a hybrid data-PINN surrogate model and shows how PINN models outperform standard data-driven neural network (NN) surrogates for every test point queried in the Design of experiments (DoE) space for the surrogate modeling problem.

The paper is divided as follows; Section 2 introduces PINNs in more detail and discusses some of the technical challenges with training PINNs. Section 3 outlines some of the important features the authors incorporate in the creation and training of PINNs to enable accurate and fast convergence. Section 4 demonstrates several problems solved using PINNs, and showcases a design optimization problem using PINN-based surrogates. Section 5 discusses how the work shown in this paper can be improved upon.

## 2 Physics Informed Neural Networks (PINNs)

### Setting up a PINN Training

Physics-informed neural networks (PINNs) leverage automatic differentiation to obtain an analytical representation of an output variable and its derivatives, given a parametrization using the trainable weights of the network. By employing the underlying static graph, it is possible to construct the differential equations that govern physical phenomena.

A PDE problem in the general form reads:

\[\mathcal{N}_{\mathbf{x}}[u]=0,\mathbf{x}\in\Omega, \tag{1}\]

\[\Phi(u(\mathbf{x}))=\mathbf{g}(\mathbf{x}),\mathbf{x}\in\partial\Omega \tag{2}\]

where \(\Phi\) can be the identity operator (Dirichlet B.C) or a derivative operator (Neumann/Robin B.C). In order to solve the PDE using the PINN method, the residual of the governing PDE is minimized, which is defined by

\[r_{\theta}(\mathbf{x})=\mathcal{N}_{\mathbf{x}}[f_{\theta}(\mathbf{x})], \tag{3}\]
where \(f_{\theta}\) is the predicted value by the network. The residual value, along with the deviation of the prediction from boundary/initial conditions, is used to construct the loss, which takes the form:

\[L(\theta)=L_{r}(\theta)+\sum_{i=1}^{M}\lambda_{i}L_{i}(\theta), \tag{4}\]

where the index i refers to different components of the loss function, relating to initial conditions, boundary conditions, and measurement/simulation data. \(\lambda_{i}\) refers to the weight coefficient of each loss term. The individual loss terms are constituted as follows:

\[L_{r}=\frac{1}{N_{r}}\sum_{i}^{N_{r}}[r(\mathbf{x}_{r}^{i})]^{2},\ L_{b}=\frac{1 }{N_{b}}\sum_{i}^{N_{b}}[\Phi(\hat{u}(\mathbf{x}_{b}^{i}))-g_{b}^{i}]^{2},\ L_{d}=\frac{1}{N_{d}}\sum_{i}^{N_{d}}[u( \mathbf{x}_{d}^{i})-\hat{u}(x_{d}^{i},t_{d}^{i})]^{2}, \tag{5}\]

where the subscripts r, b, and d refer to collocation, boundary, initial, and data points, respectively. The loss term \(L(\theta)\) can then be minimized to have the network learn the solution to the PDE described by 1,2. A popular method is to use gradient-based optimizers like Adam [22] and L-BFGS to optimize the network weights.

### Current Challenges with PINNs

Although the PINN method shows great promise, it still has a number of unresolved issues. The biggest challenges with PINNs currently lie in the scalability of the algorithms to large 3D problems as well as problems with complex nonlinearities, and unsteady problems. Some of the issues described henceforth are tackled by methods described in Section 3.

#### 2.2.1 Weak imposition of Boundary Conditions

The solution of a PDE problem must obey all initial and boundary conditions imposed on it while minimizing the residual of the governing equation. However, for neural network based solvers it is difficult to impose boundary and initial conditions in an exact manner. This is because the standard way to impose B.C in PINNs is to create a linear combination of loss functions (as described mathematically in the previous section). Each loss either describes the deviation of the network output from a specific boundary condition, or the magnitude of the residual of the governing equations. Therefore, boundary conditions are only satisfied in a weak manner. There has been research demonstrating the utility of exact imposition of boundary conditions [23, 24, 25] or creative multi-network approaches [26], such implementations are mostly problem-specific and do not generalize well.

Weak imposition of boundary conditions also creates another issue, one that is fairly common in multi-task learning and multi-objective optimization: choosing the values of loss term coefficients that make up the linear combination. Choosing these weights is a nontrivial exercise that would require calibration via hyper-parameter search, which is not feasible. Wang et al. [27] introduced a heuristic dynamic weighting algorithm to update and select these weights automatically and continuously during the training, to enable convergence to the correct answer. Additionally, there have been several other algorithms proposed to choose the correct scheme for weighting the losses [28, 29, 30]. This continues to be an active area of research in the PINNs community. Finally, methods have been proposed to impose the boundary conditions in a strong manner by manipulating the output formulations [23] or by utilizing operator networks [31].

#### 2.2.2 Difficult Optimization Problem

A second problem is the nature of the loss landscape itself, in which a reasonable local minimum is required to be found. As seen in Krishnapriyan et al. [32], Gopakumar et al. [33],Subramanian et al. [34] and Basir and Senocak [35], as well as the author's own experiments, different non-dimensional quantities (e.g. Reynolds number) in the governing equations, the number of dimensions of the problem, the point cloud/discretization, the boundary conditions and the complexity of the solution to be predicted can adversely affect the loss landscape of the neural network training. This makes the optimization challenging and can fail to find an adequate local minimum via a gradient descent-based algorithm. Recently, methods borrowing concepts from optimization theory have shown alternate formulations
(e.g. augmented lagrangian method for the loss functions) can aid the convergence properties of the training problem [35, 36]. There have also been efforts towards imposing physical constraints in an integral form [37].

#### 2.2.3 Cost of training

Constructing the PDE loss functions involves several backward passes through the network, which is a costly operation. PINNs on average take longer to train than their data-driven counterparts for exactly this reason; the computation graph of a PINN training is much more complex. Moreover, for the Navier-Stokes equations, it has been seen that although the stream function formulation provides better results (due to exact enforcement of continuity), it is costlier in terms of training time. As seen in NVIDIA's experiments [38], it can take several million iterations for the more complex problems to be solved via PINNs. To reduce the cost of training approaches such as automatic differentiation for finite difference formulations [39], or using first-order formulations [40], have been proposed. However, these solutions tend to be mostly problem-specific and do not necessarily generalize well to increased problem complexity and grid definitions. Meta-learning algorithms [41] have also recently gained significance as an effective way to reduce the cost of training neural networks on new tasks, and some of this work has been extended to PINNs [42] as well.

## 3 Important Features for Creating PINN Models

In this section, the important techniques used to create PINN-based models cost-effectively are outlined. The PINN models in subsequent sections are created by combining these features that have been found to have an effect on the accuracy of the model and the speed of training.

### Hybrid Data-Physics Training

Compared with the original PINNs method proposed by Raissi et al. [7], a plethora of research has been undertaken to improve and expand on the method [43, 44]. From these developments, the PINNs method has been applied to solve PDE-based problems of increasing complexity and dimensionality. However, the PINNs method is currently not suited for solving engineering problems often encountered in industry in a data-free manner. The optimization issues and cost of model training outlined above make the method, presently, unsuitable for use as a forward solver. To get the best of both worlds, the PINNs method can be augmented with data. Figure 1 depicts the tradeoff between using only data or only physics, and that the sweet spot lies in using both. In addition to the discussed benefit of hybrid data-physics training reducing the cost of generating data, there have been several examples showing that the inclusion of sparse solution data in the training loss function significantly improves the convergence capabilities of the PINNs method [33, 43, 45].

In this paper, we take inspiration from this and use very sparse solution data to solve 3D flow-thermal problems and inform our surrogate models with physics while creating them.

Figure 1: The spectrum of data-driven versus physics-informed models. Incorporating governing physics information into the models during creation serves as an effective form of regularization and often helps reduce the amount of data required to achieve the same accuracy levels.
### Modified Learning Rate Annealing

As described in Section 2.2.1, the learning rate annealing algorithm has proved to be very effective in mitigating the stiffness of the PINN training problem. However, utilizing this method over a broader spectrum of problems highlighted an issue with stability. The following outlines this issue:

As shown in Equation 4 the PINN loss function being optimized takes the form:

\[L(\theta)=L_{r}(\theta)+\sum_{i=1}^{M}\lambda_{i}L_{i}(\theta) \tag{6}\]

At any training step, the update to the loss coefficient is calculated [27] as

\[\hat{\lambda}_{i}=\frac{max_{\theta}|\nabla_{\theta}L_{r}(\theta)|}{|\nabla_{ \theta}L_{i}(\theta)|},i=1,....,M\]

It can be seen that if the loss \(L_{i}\) decreases much faster than \(L_{r}\) during the training, the value of \(\hat{\lambda}_{i}\) increases. This then leads to a larger coefficient for that loss term and an associated faster decay of the loss.

This instability has the unintended consequence of the optimizer getting stuck in minima where it minimizes the loss \(L_{i}\) very well but is unable to optimize for the loss of the other constraints. The proposed updated algorithm to mitigate this issue is shown in Algorithm 1. The values of thresholds are hyper-parameters, but if the inputs and outputs of the network have been normalized (using standard score normalization, for example), then selecting values between \(10^{-3}\) and \(10^{-5}\) works well in practice.

```
for update step = 1 to \(N\)do if\(L_{i}(\theta)\leq(threshold)_{i}\)then \[\hat{\lambda}_{i}=0\] else  Compute \(\hat{\lambda}_{i}\) by \[\hat{\lambda}_{i}=\frac{max_{\theta}|\nabla_{\theta}L_{r}(\theta)|}{|\nabla_{ \theta}L_{i}(\theta)|},i=1,....,M\] endif  Update weights \(\lambda_{i}\) as \[\lambda_{i}=(1-\alpha)\lambda_{i}+\alpha\hat{\lambda}_{i}\]  Update network parameters via gradient descent: \[\theta_{n+1}=\theta_{n}-\eta\nabla_{\theta}L_{r}(\theta)-\eta\sum_{i=1}^{M} \lambda_{i}\nabla_{\theta}L_{i}(\theta)\] endfor  We set the hyper-parameter \(\alpha=0.1\) and \(\eta=10^{-3}\). Threshold values are chosen somewhere between \(10^{-3}\) and \(10^{-5}\).
```

**Algorithm 1** Modified Learning Rate Annealing

For a problem with the loss function

\[L(\theta)=L_{r}(\theta)+\lambda_{neu}L_{neu}(\theta)+\lambda_{dir}L_{dir}(\theta) \tag{7}\]

where \(L_{r}(\theta)\), \(L_{neu}(\theta)\) and \(L_{dir}(\theta)\) correspond to the PDE, Neumann and Dirichlet loss respectively, Figure 2 shows the training curves for the individual losses, and the value of the adaptive coefficients when they are calculated using Algorithm 1. It can be seen that when the boundary loss terms in Figures 1(c) and 1(d) go below their thresholds (set to \(10^{-5}\)), the associated coefficients shown in Figures 1(a) and 1(b) start decaying. Following this, the PDE loss starts improving much faster. If the term \(L_{i}(\theta)\) goes above its threshold, it leads to a spike in the adaptive constant \(\lambda_{i}\) which brings it down again.
### Fourier Feature Embeddings

As described in Tancik et al. [46], Artificial Neural Networks suffer from a spectral bias problem. To overcome this, they introduced a Fourier feature embedding that allows models to capture high-frequency components of the solution effectively. This has the effect of markedly improving the ability of the networks to capture sharp gradients in the solutions, which requires the network to be able to learn high-frequency components of the solution quickly.

Following the implemenation in Tancik et al. [46], for an input vector

Figure 2: Adaptive coefficients and loss terms from Equation 7 during training. (a) Evolution of the Dirichlet loss adaptive constant during training.(b) Evolution of the Neumann loss adaptive constant during training. (c) Dirichlet B.C loss term \(L_{dir}(\theta)\) (d) Neumann B.C loss term \(L_{neu}(\theta)\) (e) The PDE loss during training. Once the values of both the adaptive constants start dropping, the PDE loss improves much more rapidly.
\[\mathbf{v}=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]\]

instead of using \(\mathbf{v}\) as the input we compute the Fourier feature mapping:

\[\gamma(\mathbf{v})=[\cos(2\pi\mathbf{b}_{1}^{T}\mathbf{v}),\sin(2\pi\mathbf{b}_ {1}^{T}\mathbf{v}),.....,\cos(2\pi\mathbf{b}_{m}^{T}\mathbf{v}),\sin(2\pi \mathbf{b}_{m}^{T}\mathbf{v})] \tag{8}\]

where m is a hyper parameter and the frequencies \(\mathbf{b}_{j}^{T}\) are selected randomly from an isotropic distribution. Then \(\gamma(\mathbf{v})\) is passed into the network.

The Fourier feature embedding was shown to be highly effective in training PINNs models by Wang et al. [47], and several results were shown for 1D and 2D problems. We extend this implementation to solve 3D flow problems via PINNs and use it to create our hybrid data-PINN surrogate for flow thermal problems.

In addition, there have been other proposed solutions for the spectral bias problem for applications to PDE problems, such as the Siren activation [48], Fourier Neural Operators [49], and weighting schemes derived from the theory of Neural Tangent Kernels (NTK) [28].

## 4 Experiments and Results

In this section, some example problems are solved using PINNs. Sections 4.1 and 4.2 solve the 3D incompressible Navier-Stokes equations through a data-assisted approach, where very sparse solution data is provided in the domain.

Section 4.3 uses a hybrid data-PINN approach to generate a surrogate model for a given design space of a heat sink with a chip underneath it, undergoing cooling via forced convection. Then, given certain constraints on the running metrics of the chip-sink setup (like max temperature in the chip), the optimal set of parameters in the Design of Experiments (DoE) space that satisfy the constraints while maximizing an objective are obtained via rapid design optimization using the created surrogate.

Details on hyper-parameters used in the model training for each experiment that follows can be found in Appendix Section A.1.

### Forward Solve of 3D Stenosis Problem

Flow through an idealized 3D stenosis geometry at a physiologically relevant Reynolds number is demonstrated, see Figure 3 for details about the geometry. To the author's best knowledge, flow through a stenosis has been solved using PINNs only at a low Reynolds number of approximately 6 (based on inlet diameter) [23]. Flow through irregular geometries has been solved at a higher Re (500), but in 2D [50]. In this paper, the stenosis problem is solved at Re 150, and in 3 dimensions.

As discussed in Section 2.2, at higher Reynolds numbers the standard PINN implementation struggles to achieve a good local minimum. This was confirmed using a standard PINN implementation. To alleviate this issue a data-assisted approach where sporadic solution data can be added throughout the domain of interest (depicted on a slice in Figure 4). The data was given in the form of concentric rings at the radii depicted on the cut plane.

#### 4.1.1 Problem Setup

The flow problem through the stenosis is solved by solving the steady-state incompressible Navier-Stokes equations:

\[\nabla\cdot\mathbf{u}=0, \tag{9}\]

\[(\mathbf{u}\cdot\nabla)\mathbf{u}=-\frac{1}{\rho}\nabla\mathbf{p}+\nu\nabla \cdot(\nabla\mathbf{u}), \tag{10}\]

subject to

\[\mathbf{u}(x_{b1})=g(x_{b1}),\,x_{b1}\in\partial\Omega_{1},\]
\[\mathbf{u}(x_{b2})=0,\,x_{b2}\in\partial\Omega_{2},\]

\[\nabla u_{i}(x_{b3})\cdot\mathbf{n}=0,\,x_{b3}\in\partial\Omega_{3},i=1,2,3\]

\[p(x_{b3})=0,\,x_{b3}\in\partial\Omega_{3}\]

where \(g(x_{b3})\) represents a profiled input to the stenosis. \(\rho\) and \(\nu\) are the density and kinematic viscosity of the fluid(air) respectively, and \(\mathbf{u}\) and p are the velocity vector and pressure respectively.

In the present problem, a parabolic profiled input is provided with a peak value inflow of 0.15 m/s. The ratio of areas of the throat to the inlet is 0.36.

The output of the network is approximated as \(G_{\theta}\), which is a 4-component output:

\[G_{\theta}=\left[\begin{array}{c}u\\ v\\ w\\ p\end{array}\right]\]

#### 4.1.2 Results

Figure 5 compares the velocity magnitude returned by the trained PINN model and Altair AcuSolve(r) through a 2D slice of the stenosis. As can be seen, the essential features of the flow are captured. Figure 5(a) and 5(b) compare the velocity and pressure profile through the center of the stenosis. The differences between the line plots are attributed to differences in mesh density between the two cases. The CFD mesh was an unstructured mesh of around 63,000 nodes with a boundary layer, while the

Figure 4: Stenosis diagram (not to scale) showing planes where solution data is provided randomly.

Figure 3: Visual description of stenosis problem
point cloud used with the PINN was around 87,000 nodes randomly distributed points except near the boundary where the sampling was finer.

Another approach that was investigated to solve the 3D stenosis problem was that of using "continuity planes" as defined by Hennigh et al. [38] in their experiments solving 3D flow problems using PINNs. In this approach, the authors added constraints on the mass flow through a plane and added these constraints to the loss function. While this approach was found to aid the convergence of the PINN model to the correct solution, there were several issues found to exist with this method:

1. It is difficult to generate continuity planes for complex geometries such as those shown in Sections 4.2 and 4.3.
2. The quality of the solution from the PINN depends heavily on the integration scheme used to calculate the mass flow rate, and the fineness of the points on the continuity plane.

Figure 5: Solution Comparison. (a) Altair AcuSolve® Solution to stenosis problem (b) PINN forward solve to stenosis problem.

Figure 6: Centerline solution comparisons: PINN versus Altair AcuSolve® (a) Total Velocity Comparison (b) Pressure Comparison
Hence, in the next section, random and sparsely distributed data was used in the domain to aid convergence.

### Flow over a Printed Circuit Board (PCB)

#### 4.2.1 Problem Setup

Flow over a PCB consisting of a heat sink, chip, and capacitor is solved at a Reynolds number of approximately 1500, based on the length of the PCB board and air as the fluid. The geometry and flow orientation are shown in Figure 7. This represents a forced convection problem common in electronics design and is a challenging problem for PINNs because it is in 3D, with a complex geometry and large gradients involved.

Let \(D\) represent the set of all nodes in the domain. To train the PINN model, the CFD solution was first computed. Next, 1% of the nodes in the solution domain were randomly selected (call this set \(D_{1}\subset D\))). This is a selection of roughly 2,300 node points (from a mesh of roughly 230,000 nodes). The experiment was then divided into three parts:

1. **Case A**: A network was trained on the CFD solution at all points in \(D_{1}\) (i.e \(\forall\mathbf{x}\in D_{1}\)) following which the physics of the problem was **enforced at every node location** in \(D\) (i.e \(\forall\mathbf{x}\in D\)), by including the physics-based loss in the training, and then the network was asked to predict the solution in the entire domain \(D\).
2. **Case B**: A network was trained on the CFD solution at the points contained in \(D_{1}\) (i.e, \(\forall\mathbf{x}\in D_{1}\)) **without any physics enforcement** and then asked to predict the solution in the entire domain (i.e \(\forall\mathbf{x}\in D\)).
3. **Case C**: Finally, the same experiment as Case A was repeated but with a new set \(D_{2}\) consisting of only 0.2% of the nodes in \(D\), which were again randomly selected.

The governing equations for this problem are the Reynolds Averaged Navier-Stokes Equations:

\[\nabla\cdot\mathbf{u}=0, \tag{11}\]

\[(\mathbf{u}\cdot\nabla)\mathbf{u}=-\frac{1}{\rho}\nabla\mathbf{p}+(\nu+\nu_{t })\nabla\cdot(\nabla\mathbf{u}), \tag{12}\]

\(\rho\), \(\nu\), and \(\nu_{t}\) represent the density, kinematic viscosity and eddy viscosity of the system. The inflow is set to a constant velocity of 0.15 m/s and the outflow is set to the stress-free condition. It should be noted that in the current study, eddy viscosity is obtained directly from the CFD solver using the Spalart-Allmaras turbulence model. Turbulence modeling in PINNs is a field of active research with a few articles investigating it [38, 51, 52], and it is a future work to effectively incorporate turbulence models into PINN-based models.

Figure 7: Geometry of a PCB with a chip, sink, and capacitor assembly.
#### 4.2.2 Results

Figure 8 shows the ANN predictions for the different cases. It is evident that by using sparse data, the network is able to better converge toward the CFD solution (shown in Figure 8d) using the physics-based regularizer. However, as evident in Figure 8c, the network failed to converge to a physical solution when the amount of data provided was insufficient, highlighting the importance of a certain amount and fineness of the required data. Table 1 shows the Mean Squared Errors (MSE) for each experiment, for the velocity and the pressure, taking the CFD solution as the ground truth. The MSE is calculated as

\[\text{MSE}=\sqrt{\frac{\sum_{i=1}^{N_{nodes}}(x_{i,pred}-x_{i,truth})^{2}}{N_{ nodes}}} \tag{13}\]

Figure 9 shows the fraction of node points for each case that are above a certain Mean Absolute Error (MAE) value. The lower the fraction, the better the solution. We note from Figure 9 that even for Case A, there are outliers to the solution where the MAE is relatively high, indicating poor convergence to the solution at those nodes. The convergence of PINNs toward the correct solution for highly nonlinear systems is an open and challenging problem, especially in 3 dimensions. Nonetheless, these results open exciting possibilities about using physics-based regularizers in the future and represent a step forward for solving the 3D Navier-Stokes Equations at high Reynolds Numbers using PINNs. Furthermore, data generation costs to create surrogate models using PINNs can be greatly reduced by providing solution data on a coarser grid and solving the physics on a finer grid.

### Surrogate Modeling and Design Optimization of a Heat Sink

In this section, the PINNs surrogate modeling technique is demonstrated for rapid design optimization of a heat sink assembly. The assembly utilizes a chip that generates heat and a fin-type heatsink on top to dissipate heat into the surrounding fluid. The chip-heatsink assembly is cooled by forced convection of air. The geometry and setup are shown in Figure 10.

The goal is to optimize the heat sink design and the running conditions of the assembly, subject to feasibility constraints placed on chip temperature and channel pressure drop. This represents a common design optimization problem in electronics cooling. More specifically, if \(\dot{Q}_{src}\) is the total power being generated by the chip, the optimization problem can be framed as follows:

\[\text{Maximize }\dot{Q}_{src}\text{ s.t} \tag{14}\]

\[\text{Pressure drop across the heat sink channel ($\Delta$P)}\leq 11\text{ Pa} \tag{15}\]

\[\text{Maximum temperature anywhere on the chip }\leq 350\text{ K} \tag{16}\]

The pressure at the outflow is fixed to 0 Pa, and the pressure drop across the heat sink channel is hence calculated as the average pressure over the inflow of the channel:

\[\Delta\text{P}=\overline{\text{P}_{\text{inlet}}} \tag{17}\]

The term to be maximized \(\dot{Q}_{src}\) is also one of the design axes and an input parameter(P3) to the network.

The design variables that can be altered for this present optimization are:

\begin{table}
\begin{tabular}{|c|c|c|c|} \hline Case & Description & MSE (Velocity) & MSE (Pressure) \\ \hline Case A & 1\% domain data + physics & **0.0135** & **0.0037** \\ \hline Case B & 1\% domain data only & 0.0222 & 0.00472 \\ \hline Case C & 0.2\% domain data + physics & 0.0245 & 0.00545 \\ \hline \end{tabular}
\end{table}
Table 1: Mean Squared Errors (MSE) for velocity and pressure for cases A,B, and C
* Inflow Velocity
* Fin height
* Source term in the chip (has to be maximized)

The upper and lower limits of each of the design variables mentioned above are summarized in Table 2. The inlet velocity is set based on typical values found in literature [53] and corresponds to a Reynolds number range of Re 10,300 to Re 24,000.

The governing equations solved for this conjugate heat transfer problem are the same as in Section 4.2 for the flow problem, subject to no-slip boundary conditions on the chip-heatsink assembly with a

Figure 8: Neural Network (NN) prediction with and without physics, for very coarse data supplied on a plane through the domain. (a) **Case A:** Trained on 1% data and physics (b) **Case B:** Trained on 1% solution data only (c) **Case C:** Trained on 0.2% data and physics (d) True Solution from CFD solver

Figure 9: Node fractions of points above a certain MAE value, for each case. (a) MAE of Velocity (b) MAE of Pressure
variable freestream inflow velocity, causing forced convection. As in Section 4.2, the eddy viscosities are taken from the CFD solutions.

The energy equation in both fluid and solid reads:

\[k\nabla^{2}T+\dot{q}_{src}-\rho s\mathbf{u}\cdot\nabla T=0, \tag{18}\]

where T represents the temperature, \(\dot{q}_{src}\) represents the volumetric source term, and \(k\) and \(s\) are the conductivity and specific heat of the material respectively. At the interface between the fluid and solid domain (fluid-sink, sink-chip, and fluid-chip) the interface condition is applied by minimizing the following loss terms as shown in [54];

\[L_{flux}=\frac{1}{N_{int}}\sum_{i=1}^{N_{int}}(f_{d_{1}}(\mathbf{u}(x_{i})) \cdot\mathbf{n}_{d_{1}}+f_{d_{2}}(\mathbf{u}(x_{i}))\cdot\mathbf{n}_{d_{2}})^ {2}, \tag{19}\]

\[L_{val}=\frac{1}{N_{int}}\sum_{i=1}^{N_{int}}(\mathbf{u}_{d_{j}}(x_{i})- \overline{\mathbf{u}_{d_{j}}(x_{i})})^{2}, \tag{20}\]

where \(\mathbf{n}_{d1}=-\mathbf{n}_{d2}\) and j=1,2. The average is taken over j. \(d_{1}\) and \(d_{2}\) refer to the domains on both sides of the interface, and \(N_{int}\) is the number of node points on the interface.

#### Model Creation and Evaluation

The sampling of the above Design of Experiments (DoE) space is done via an efficient space-sampling method to optimally fill the DoE space [55]. The sampled DoE space for training is shown in Figure 11, along with the location of points at which the surrogate model is tested. The reader is referred to Section A.3.1 for a complete tabular description of the DoE space. Note that for this

\begin{table}
\begin{tabular}{|c|c|c|c|} \hline Parameter No. & Parameter Name & Lower Value & Upper Value \\ \hline P1 & Inflow Velocity (\(m/s\)) & 3 & 7 \\ \hline P2 & Fin Height (\(mm\)) & 15 & 23 \\ \hline P3 & Source Term (\(W\)) & 30 & 60 \\ \hline \end{tabular}
\end{table}
Table 2: Design of Experiments space axes ranges for the heat sink design optimization

Figure 10: Basic problem geometry and flow depiction
example, we use full field data at each DoE point to train the surrogate as opposed to a small fraction of it (like in Section 4.1 and 4.2), as the objective is to get a surrogate that is as accurate as possible. Table 3 shows the MSE for the predictions by the hybrid data-PINN model at the test points, calculated w.r.t the CFD solution at the same mesh resolution. Also shown is the MSE for predictions by a standard data-driven NN without leveraging key features described in Section 3, which are used extensively in industry for surrogate modeling applications. The hybrid data-PINN model outperforms the standard data-driven NN for all predictions. Section A.4 shows some more qualitative comparisons between test point results from the PINNs model versus standard data-driven NNs.

#### Solving the Design Optimization Problem

The surrogate model is used to solve the design optimization problem described in Equations 14-16. The goal is to show that the surrogate model can accurately predict the solution in the entire DoE space by returning a solution that satisfies all applied constraints while maximizing the objective. The created surrogate models are interfaced with an optimizer that solves a generic constrained optimization problem via an iterative process, described thoroughly in Appendix Section A.2.

Each snapshot in Figure 12 represents a design iteration, and each particle represents a point in the DoE space. Each axis of a plot represents a parameter axis.

Figure 11: Training and testing points in the 3D DoE space

Figure 12: Design optimization iterations of the heat sink problem (a) Iteration 0 (b) Iteration 5 (c) Iteration 10
For the given constraints, the particles converge to a much smaller region of the DoE space. The design point returned by the optimizer in this case is:

**Inflow Velocity**: 6 m/s

**Chip Power**: 50W

**Fin Height**: 17mm

To test that the result satisfies the constraints, the returned design point is solved by the Altair AcuSolve(r), at the same mesh fineness and another mesh with 10x fineness, keeping all essential mesh features such as boundary layers and refinement zones. As shown in Figures 13 and 14, not only does the given design point satisfy the design constraints, but the finer mesh solution is very close to the coarser solution, and a little tweaking of the design point using CFD with the higher resolution mesh will yield a highly optimized solution to the designer. This optimization is done several orders of magnitude faster than if using traditional CFD, and the reader is referred to Appendix Section A.3.2 for a quantitative description of the same.

## 5 Conclusions and Future Work

In this paper, Physics Informed Neural Networks were used to solve the 3D Navier-Stokes equations in a data-assisted setting, for complex geometries with realistic physical parameters. It was shown that even for problems being solved at high Reynolds Numbers in 3D, PINNs can be trained to produce a

Figure 13: Temperature plot through a slice at the rear of the sink (from bottom to top).

The comparison between the high-fidelity solution on the fine mesh and the PINN prediction on a coarser mesh shows good agreement.

\begin{table}
\begin{tabular}{l c c c} \hline \hline  & **Velocity MSE** & **Pressure MSE** & **Temperature MSE** \\ \hline
**Test Point 1** & & & \\ Hybrid data-PINN Model & **0.65** & **2.62** & **1.81** \\ Standard Data-driven NN & 0.93 & 2.83 & 2.05 \\ \hline
**Test Point 2** & & & \\ Hybrid data-PINN Model & **0.39** & **1.19** & **2.67** \\ Standard Data-driven NN & 0.58 & 1.42 & 2.97 \\ \hline
**Test Point 3** & & & \\ Hybrid data-PINN Model & **0.76** & **3.31** & **1.86** \\ Standard Data-driven NN & 1.10 & 3.51 & 2.18 \\ \hline
**Test Point 4** & & & \\ Hybrid data-PINN Model & **0.33** & **0.99** & **2.87** \\ Standard Data-driven NN & 0.52 & 1.19 & 3.15 \\ \hline \hline \end{tabular}
\end{table}
Table 3: MSE for 4 test points shown in Table 5. The PINN-based model consistently outperforms the standard data-driven NN on all test points.
good solution in the presence of very sparse solution data randomly scattered in the solution domain. However, using too little solution data causes the model to converge to an unphysical solution. PINNs were also demonstrated for 3D flow-thermal surrogate modeling and the PINN-based surrogates consistently outperformed standard data-driven NN on test case examples. The PINN surrogates were also interfaced with a design optimization algorithm to solve a constrained optimization problem. This optimization returned a design point that when solved with high-fidelity CFD was consistent with the requirements of the design constraints, highlighting the suitability of the method to produce surrogates for 3D flow-thermal surrogate modeling problems.

There are multiple avenues through which the work shown in this paper can be improved. Research has to be done to improve convergence and offer guarantees of PINNs training toward the local minimums that represent physical solutions, in a data-free manner. This will further reduce data requirements for the creation of physically consistent PINN models which can greatly improve their surrogate modeling capabilities, by reducing the cost of training and improving predictive accuracy. Further work needs to be done to investigate turbulence modeling in PINNs so that high Reynolds number problems can be solved in a data-free manner. There are also many themes like uncertainty quantification [56, 57, 58] of surrogates and effective surrogate modeling of different geometries [59, 60, 61] that are active fields of research in PINNs, which could be included in future works that build on these results.

## Acknowledgements

This research did not receive any specific grant from funding agencies in the public or not-for-profit sectors, or from any external commercial entities. The authors gratefully acknowledge the use of Altair Engineering Inc.'s computing facilities for running experiments.

## CRediT authorship contribution statement

**Saakaar Bhatnagar:** Formal Analysis, Investigation, Methodology, Software, Validation, Writing-original draft. **Andrew Comerford:** Conceptualization, Investigation, Project Administration, Supervision, Writing- review and editing **Araz Banaeizadeh:** Conceptualization, Project Administration, Supervision, Writing- review and editing
