Graph Neural Network Based Surrogate Model of Physics Simulations for
  Geometry Design
Computational Intelligence (CI) techniques have shown great potential as a
surrogate model of expensive physics simulation, with demonstrated ability to
make fast predictions, albeit at the expense of accuracy in some cases. For
many scientific and engineering problems involving geometrical design, it is
desirable for the surrogate models to precisely describe the change in geometry
and predict the consequences. In that context, we develop graph neural networks
(GNNs) as fast surrogate models for physics simulation, which allow us to
directly train the models on 2/3D geometry designs that are represented by an
unstructured mesh or point cloud, without the need for any explicit or
hand-crafted parameterization. We utilize an encoder-processor-decoder-type
architecture which can flexibly make prediction at both node level and graph
level. The performance of our proposed GNN-based surrogate model is
demonstrated on 2 example applications: feature designs in the domain of
additive engineering and airfoil design in the domain of aerodynamics. The
models show good accuracy in their predictions on a separate set of test
geometries after training, with almost instant prediction speeds, as compared
to O(hour) for the high-fidelity simulations required otherwise.
Jian Cheng Wong, Chin Chun Ooi, Joyjit Chattoraj, Lucas Lestandi, Guoying Dong, Umesh Kizhakkinan, David William Rosen, Mark Hyunpong Jhon, My Ha Dao
# Graph Neural Network Based Surrogate Model of Physics Simulation for Geometry Design

###### Abstract

Computational Intelligence (CI) techniques have shown great potential as a surrogate model of expensive physics simulation, with demonstrated ability to make fast predictions, albeit at the expense of accuracy in some cases. For many scientific and engineering problems involving geometrical design, it is desirable for the surrogate models to precisely describe the change in geometry and predict the consequences. In that context, we develop graph neural networks (GNNs) as fast surrogate models for physics simulation, which allow us to directly train the models on 2/3D geometry designs that are represented by an unstructured mesh or point cloud, without the need for any explicit or hand-crafted parameterization. We utilize an encoder-processor-decoder-type architecture which can flexibly make prediction at both node level and graph level. The performance of our proposed GNN-based surrogate model is demonstrated on 2 example applications: feature designs in the domain of additive engineering and airfoil design in the domain of aerodynamics. The models show good accuracy in their predictions on a separate set of test geometries after training, with almost instant prediction speeds, as compared to O(hour) for the high-fidelity simulations required otherwise.

 Graph neural network, fast surrogate model, physics simulation

## I Introduction

Computational simulations of physics have become an essential component of modern science and engineering, and the field is growing rapidly. Physics simulations are frequently used as a model of reality to provide insights for understanding and predicting the behavior of physical systems at various levels of detail and accuracy. Simulation models can be computationally very expensive, due to their complex, multi-physics, and multi-scale nature. Therefore, sole use of high-fidelity simulation models for many tasks--such as material discovery, drug discovery, digital twin, computer-aided design--where one needs to quickly survey a vast number of possible scenarios and/or geometries, can quickly become intractable.

CI techniques have shown great potential to learn the nonlinear input-output relationship from existing simulation results, allowing for fast predictions, albeit at the expense of some accuracy. They are used as fast surrogate models of physics simulation to support discovery, design optimization, and decision-making process. Traditional surrogate models typically only predict a single output or a few outputs. The recent quest is to make predictions of the entire finely resolved physical field using advanced deep learning techniques. Deep convolutional neural networks, such as U-Net, have become a popular choice for the high dimensional surrogate modelling task when the data can be easily represented as 2D pixel or 3D voxel [1], and have been successfully applied to engineering problems such as flow past an airfoil [2, 3, 4, 5], micro fluids channel flow [6], stress distribution in a 3D printed part [7, 8], and topological optimization [9].

Many scientific and engineering problems involve geometry design. It is desirable for the geometry, e.g., from molecular structure to the shape of airfoil, to be precisely represented by the surrogate models. The voxel-type representation however has its limitation when describing geometry designs with curvature and/or fine features, and hence this may impede the performance and applicability of convolutional neural network-based surrogate models. Usually, data augmentation will be required when training convolutional neural network models, e.g., to incorporate invariances. On the other hand, graph-based methods represent a promising candidate for modeling geometry [10, 11]. For example, in molecular and materials science domain,
GNNs have been explored for learning the topology of atomic and molecular structures and predicting their properties [12, 13]. GNNs have also shown promise in predicting the time evolution of solid and fluid dynamic systems [14].

In that context, we develop GNN-based surrogate models of physics simulation for geometry design. By utilizing graph representation, the model can more precisely describe the change in geometry with complex features and predict the consequences, which is essential for the design and optimization task. We evaluate the GNN-based surrogate models on 2 engineering design problems, i.e., feature design in the additive manufacturing domain and airfoil design in the aerodynamics domain. The models are not only capable of making fast prediction of an overall performance indicator, but also of predicting the entire physical field of a geometry design, such as residual stress distribution of a 3D printed part and surface pressure distribution on an airfoil, with orders of magnitude speed-up relative to typical high-fidelity simulations.

## II Statement of The Problem

We consider the following surrogate modelling task. The goal is to develop a CI model that allows the computationally efficient evaluation of the physical quantities of interest, e.g., the physical field or an aggregated performance indicator, due to the change in geometry, in addition to any possible change in operating condition. The CI model is therefore required to map the input geometry and operating condition onto the output physical quantities of interest.

The data from physics simulations (e.g., finite element, finite volume, or finite difference methods) are used for training the CI model. Typically, simulations numerically compute the physical quantities on a _mesh_, i.e., a network of discretized cells and node points covering the continuous problem domain. In a nutshell, the simulation mesh can be classified into structured and unstructured mesh types. The structured mesh has a regular connectivity on homogenous, voxel-like cells. On the other hand, the unstructured mesh allows for mixed types of mesh cells and irregular connectivity, making it a preferred choice for handling complex geometries. Moreover, the amount of node points and connecting edges may change significantly given different geometry designs of varying shape, feature, and size. The CI model should be generic enough to learn from such unstructured mesh simulation data.

## III Graph Neural Network Based Surrogate Model of Physics Simulation

### _Graph representation of simulation data_

The unstructured simulation mesh with node points \(V\) connected by mesh edges \(E\) can be represented by a graph \(G=\left(V^{G},E^{G}\right).\) Each node \(i\in V^{G}\) is associated with _node feature_\(\mathbf{v}_{i}\). Each edge \(ij\in E^{G}\) connecting node \(i\) to node \(j\) is associated with _edge feature_\(\mathbf{e}_{ij}\). The graph contains bidirectional mesh edges, i.e., for each edge \(ij\in E^{G}\) there exist an edge \(ji\in E^{G}\) connecting from node \(j\) to node \(i\). Each node point is associated with a positional coordinate \(\mathbf{x}_{i}\) in Cartesian coordinate system, which is to be encoded into both node and edge features. To achieve spatial equivariance, we encode relative coordinate \(\mathbf{x}_{i}-\mathbf{x}_{ref(G)}\) to a reference point \(\mathbf{x}_{ref(G)}\) in the given graph \(G\) into the node feature \(\mathbf{v}_{i}\). We encode relative displacement vector \(\mathbf{x}_{li}=\mathbf{x}_{j}-\mathbf{x}_{i}\) and its L2 norm \(\left\|\mathbf{x}_{ji}\right\|_{2}\) into the edge feature \(\mathbf{e}_{ji}\). The graph is a very general representation in the way that we may encode many other problem related parameters, such as the change in operating condition, into the node or edge features.

### _Graph neural networks for learning simulation data_

For this work, we use a GNN to learn the mapping from the input graph \(G=\left(V^{G},E^{G}\right)\) onto the output physical quantities that we want to model. We adopt an encode-process-decode type architecture (Figure 1), in accordance with prior work in literature [15, 14, 10].

**Encoder.** The encoder embeds each node and edge in the input graph into a latent vector of size \(n_{l}\), using the encoder MLPs \(\varepsilon^{L}\) and \(\varepsilon^{V}\) for edge feature \(\mathbf{e}_{ji}\) and node feature \(\mathbf{v}_{i}\) respectively, i.e.,

\[\mathbf{e}^{\prime}_{ji}\leftarrow\varepsilon^{E}(\mathbf{e}_{ji}), \tag{1a}\]

Figure 1: (a) Comparison of unstructured and structured / voxel meshes on an example geometry from the feature design problem (Section IV-A-1). The voxels are unable to describe features such as slope and hole, despite having a much denser mesh. (b) Example geometry from the airfoil design problem (Section IV-A-2) with nodes and edges. (c) Illustrative schematic of the GNN architecture employed.
\[\mathbf{v^{\prime}}_{i}\leftarrow\mathbf{\varepsilon}^{\nu}(\mathbf{v}_{i}). \tag{1b}\]

**Processor.** The processor recursively updates the embedded latent vectors at each node and edge with \(L\) message passing steps. Each message passing step uses a separate block of processor MLPs \(\rho^{\mathbf{\varepsilon}}\) and \(\rho^{\nu}\) for the latent edge and node features:

\[\mathbf{\dot{e}^{\prime}}_{jl}\leftarrow\rho^{\mathbf{\varepsilon}}(\mathbf{e}^{\prime}_{ jl},\mathbf{v^{\prime}}_{j},\mathbf{v^{\prime}}_{i}), \tag{2a}\] \[\mathbf{\dot{v}^{\prime}}_{i}\leftarrow\rho^{\nu}\big{(}\mathbf{v^{\prime}}_{i },\sum_{j}\mathbf{\dot{e}^{\prime}}_{jl}\big{)}. \tag{2b}\]

The edge processor MLP \(\rho^{\mathbf{\varepsilon}}\) first computes the update for latent edge features by taking the concatenated features of edge and its receiver and sender nodes from the previous step as input. Then, the node processor MLP \(\rho^{\nu}\) takes the concatenated features of node from previous step and all incoming edges (with a sum pooling) at current step as input. Finally, the updates are added back to the latent edge and node features from previous step, creating a residual connection:

\[\mathbf{e^{\prime}}_{jl}\leftarrow\mathbf{e^{\prime}}_{jl}+\mathbf{\dot{e}^{ \prime}}_{jl}, \tag{3a}\] \[\mathbf{v^{\prime}}_{i}\leftarrow\mathbf{v^{\prime}}_{i}+\mathbf{\dot{v}^{ \prime}}_{i}. \tag{3b}\]

**Decoder.** The latent node features are then passed to 2 separate decoders, i.e., a graph decoder MLP \(\delta^{\mathbf{\varepsilon}}\) and a node decoder MLP \(\delta^{\mathbf{\varepsilon}}\):

\[\mathbf{y}_{\mathbf{\varepsilon}}\leftarrow\delta^{\mathbf{\varepsilon}}(\tfrac{1}{n(\mathbf{ \varepsilon}^{\mathbf{\varepsilon}})}\sum_{i\in\mathbf{\varepsilon}^{\mathbf{\varepsilon} }}\mathbf{v^{\prime}}_{i}), \tag{4a}\] \[\mathbf{y}_{i}\leftarrow\delta^{\mathbf{\varepsilon}}(\mathbf{v^{\prime}}_{i},\mathbf{y}_{\mathbf{ \varepsilon}}). \tag{4b}\]

The graph decoder MLP extracts the graph level feature \(\mathbf{y}_{\mathbf{\varepsilon}}\) from the average pooling of all graph nodes. We denote \(n(\mathbf{V}^{\mathbf{\varepsilon}})\) as the total number of nodes in \(\mathbf{\varepsilon}\). For a graph level prediction task, \(\mathbf{y}_{\mathbf{\varepsilon}}\) becomes the output. For node level prediction task, \(\mathbf{y}_{\mathbf{\varepsilon}}\) will then be concatenated with the latent node features \(\mathbf{v}^{\prime}_{i}\) and passed to the node decoder MLP \(\delta^{\mathbf{\varepsilon}}\) to predict the physical quantities \(\mathbf{y}_{i}\) that we want to model.

## IV Results and Discussion

We evaluate the GNN-based surrogate models on 2 engineering problems, namely the feature design in the additive manufacturing domain and airfoil design in the aerodynamics domain.

### _Description of problem domain and input encoding_

#### Iv-A1 Feature design in additive manufacturing

The laser powder bed fusion process (LPBF) in additive manufacturing can fabricate parts with intricate geometries, providing a high level of design freedom for designers to improve the functional performance of products. While most existing studies in additive manufacturing domain have focused on simple geometry, e.g., regular, bulky shapes, such surrogate models may not be sufficient for the co-design of product and process. A more general surrogate model that can quickly evaluate the functional performance due to any change in shapes and features is needed. With the aim of predicting a wide range of geometries that can arise during the design process, [7] trained a U-Net model on a large set of 3D geometries comprising different combinations of basic geometric features, e.g., circular struts, square struts, and walls. The U-Net model however can have poor voxelization when representing certain features such as circular and thin struts.

In this work, we train a GNN on a similar set of geometries. The objective of the surrogate model is to predict the physical field, e.g., residual stress (measured through the scalar von Mises stress) on every node point, for features with different shape, sizes, orientations, and intersections. We utilize graph (unstructured mesh) representation which can more precisely describe these features (Figure 1a). Moreover, we do not need to perform the data augmentation as described in [7], since our encoding already accounts for spatial equivariance. 12 geometries (test set 1) are reserved for testing. The remaining 458 geometries (train set) are used for model training. Furthermore, 7 challenging geometries (test set 2) are designed to evaluate the out-of-distribution performance of the model: 3 contained more complex combinations of the circular struts, square struts, and walls; 2 are the original and optimized designs of an engineering bracket; the other 2 are the original and optimized designs of

Figure 2: In blue: selected feature designs used for model training (train set). In red: 7 challenging geometries (test set 2) are designed to evaluate the out-of-distribution performance of the model.
a bike stem. Selected feature designs are displayed in Figure 2. A high fidelity LPBF process full order simulation [16] was used to generate the ground truth for all these geometries.

When computing relative coordinates given a graph \(G\), the reference point \(\mathbf{x}_{ref(G)}\) is chosen to be \(\left(x_{median}\cdot y_{median},0\right)\) which is the median node points in \(x\) and \(y\) axes. This is because the geometry is fixed at the printing direction which is along \(z\) axis. In addition to the relative coordinate, we encode its L1 norm \(\left\|\mathbf{x}_{i}-\mathbf{x}_{ref(G)}\right\|_{1}\) into the node feature \(\mathbf{v}_{l}\). A node \(i\) may belong to multiple cells which are from a different cell type. We further encode cell types into \(\mathbf{v}_{l}\) using one hot vectors \(\mathbf{I}_{cell,type_{l}}\). The \(\mathbf{v}_{l}\) also includes the number of connecting nodes, \(n_{e_{l}}\). The additional information may help the GNN to better identify abstract features across different graphs. The input features and target output are normalized to an appropriate range for GNN training.

#### Iii-A2 Airfoil design in aerodynamic

The airfoil design is an important problem in aerodynamic [5]. Variables of interest of the design problem include the pressure distribution on the airfoil surface, the total lift and drag coefficients. As an airfoil is a thin object with a sharp and long tail, precise representation of its shape is crucial for accurate prediction of its aerodynamic performance. In this work, we train the GNNs to predict the surface pressure on airfoil, as well as 2 derived performance indicators: drag and lift coefficients \(C_{D}\), \(C_{L}\). We leverage the huge database of Reynolds-averaged Navier-Stokes (RANS) simulations generated by [3], comprising of 1,505 2D airfoil shapes obtained from the UIUC Airfoil Coordinates Database [17]. Each airfoil shape was simulated with several freestream conditions \(\left(u_{o},v_{o}\right)\). We extracted and post-processed a total of 6,372 simulation solutions (train set) with 1,448 airfoils for model training. Specifically, we construct a graph for the airfoil shape by using the surface node points. Each node is only connected to its nearest left and right neighbor nodes along the airfoil surface (Figure 0(b)). We also extracted and post-processed 90 simulation solutions (test set 1) with 30 airfoils that are not used for training, to evaluate the generalization of the model.

When building the graph \(G\) for a given airfoil, we use \(\mathbf{x}_{ref(G)}\)\(=\) (0,0) as reference point. This is because all the airfoil shapes are already scaled to \(x\in\) [0,1], with the leading edge positioned at (0,0). We also encode the given freestream condition \(\left(u_{o},v_{o}\right)\) into the node features \(\mathbf{v}_{l}\), for all the nodes \(i\in V^{G}\). In addition, we use one hot vectors \(\mathbf{I}_{u/l_{1}}\) to indicate whether a given node \(i\) belongs to the upper or lower surface. The input features are normalized to an appropriate range for GNN training. Like [3], we normalize the target pressure by dividing it with the freestream velocity \(vel\)\(=\)\(u_{o}\)\({}^{2}+\)\(v_{o}\)\({}^{2}\), followed by a de-mean. We do not apply normalization to the target \(C_{D}\) and \(C_{L}\).

### _Model architecture, training, and evaluation_

We build a node level target model for residual stress from the feature design dataset, and a node level target model for pressure from the airfoil design dataset. For the same dataset, we build separate graph level target model for \(C_{D}\) and \(C_{L}\). The architecture and training configurations used for these 4 models are summarized in Table I. For each model, we tested several architecture variants and training parameters. It was found that they are robust to many choices, such as the depth \(n_{d}\) and width \(n_{w}\) of each MLP, the size of latent features \(n_{l}\), the number of message passing steps \(L\), and the batch size, as long as they are within a certain range.

All the MLPs use _sine_ activation, which is found to be a better choice compared to other activations such as _tanh_ and ReLU, in terms of solution quality. However, a _linear_ activation is used at the final output layer for pressure, \(C_{D}\) and \(C_{L}\) models since their values span from negative to positive. The residual stress is strictly nonnegative hence we use a ReLU activation at the final output layer of the model.

We initialize the models using _He_ method and train them with ADAM optimizer for a fixed number of epochs. We use the mean absolute error (MAE) loss with L1 regularization, which tends to give better training and generalization

\begin{table}
\begin{tabular}{c c c c c c c c c} \hline \hline Problem & No. training & Output & Input graph: & Input graph: & MLP architecture (depth \(n_{d}\), & No. & Training \\ domain & training & output & type & node \(\mathbf{v}_{l}\) & edge \(\mathbf{e}_{l}\) & width \(n_{w}\) output size \(n_{l}\) for & message & (loss, no. epoch, \\  & sample & (graph) & & & & & encoder \(\pi^{x}\), \(\varepsilon^{y}\), processor \(\varepsilon^{y}\), & passing & batch size, initial \\  & (graph) & & & & & & \(\rho^{y}\), and decoder \(\delta^{c}\), \(\delta^{y}\) & step, \(L\) & learning rate) \\ \hline
1. Feature design & 458 & Residual & Node & \(\mathbf{x}_{l}-\mathbf{x}_{ref(G)}\), & \(\mathbf{x}_{l_{1}}\), & \(\left\|\mathbf{x}_{i}\right\|_{2}\) & \(\varepsilon^{x}\), \(\varepsilon^{y}\), \(\rho^{x}\), \(\varepsilon^{y}\): & 6 & MAE \(+\) L1 reg. \\ design & \(\uparrow\) & stress & level & \(\left\|\mathbf{x}_{l}-\mathbf{x}_{ref(G)}\right\|_{1}\) & & \(n_{d}\)\(=\)\(4,n_{w}\)\(=\)\(64\), & & epoch=\(2000\), \\ (3D) & 46 – 3492 & & & & \(\mathbf{I}_{cell,type_{l}}\), \(n_{e_{l}}\) & & \(\delta^{c}\): \(n_{l}\)\(=\)\(4,n_{w}\)\(=\)\(64\), \(n_{i}\)\(=\)\(4\) & & batch size=\(16\), \\  & nodes & & & & & \(\delta^{c}\): \(n_{u}\)\(=\)\(4,n_{w}\)\(=\)\(64\), \(n_{i}\)\(=\)\(1\) & & initial\_r=\(5\)e-4 \\ \hline Airfoil design & 6372 & Pressure & Node & \(\mathbf{x}_{l}-\mathbf{x}_{ref(G)}\), & \(\mathbf{x}_{l_{1}}\), & \(\left\|\mathbf{x}_{i}\right\|_{2}\) & \(\varepsilon^{x}\), \(\varepsilon^{y}\), \(\rho^{x}\), \(\varepsilon^{y}\): & 5 & MAE \(+\) L1 reg. \\ design & \(\uparrow\) & level & & & \(n_{d}\)\(=\)\(5,n_{w}\)\(=\)\(64\), \(n_{i}\)\(=\)\(64\) & & epoch=\(3000\), \\ (2D) & 194 – 660 & & & & & \(\delta^{c}\): \(n_{u}\)\(=\)\(5,n_{w}\)\(=\)\(64\), \(n_{i}\)\(=\)\(4\) & & batch size=\(32\), \\  & nodes & & & & & \(\delta^{c}\): \(n_{u}\)\(=\)\(5,n_{w}\)\(=\)\(64\), \(n_{i}\)\(=\)\(1\) & & initial\_r=\(5\)e-4 \\ \hline Airfoil design & 6372 & Drag & Graph & \(\mathbf{x}_{l}-\mathbf{x}_{ref(G)}\), & \(\mathbf{x}_{l_{1}}\), & \(\left\|\mathbf{x}_{i}\right\|_{2}\) & \(\varepsilon^{x}\), \(\varepsilon^{y}\), \(\rho^{x}\), \(\varepsilon^{y}\): & 5 & MAE \(+\) L1 reg. \\ design & \(\uparrow\) & coefficients, & level & \(\mathbf{I}_{u/l_{1}}\), & \(\mathbf{I}_{u/l_{1}}\), & \(\mathbf{I}_{u/l_{1}}\), & \(n_{u}\)\(=\)\(5,n_{w}\)\(=\)\(64\), \(n_{i}\)\(=\)\(64\) & & epoch=\(500\), \\ (2D) & 194 – 660 & \(C_{D}\) & & & & \(\delta^{c}\): \(n_{u}\)\(=\)\(5,n_{w}\)\(=\)\(64\), \(n_{i}\)\(=\)\(1\) & & batch size=\(64\), \\  & nodes & & & & & \(\delta^{c}\): \(n_{u}\)\(=\)\(5,n_{w}\)\(=\)\(64\), \(n_{i}\)\(=\)\(1\) & & initial\_r=\(5\)e-4 \\ \hline Airfoil design & 6372 & Lift & Graph & \(\mathbf{x}_{l}-\mathbf{x}_{ref(G)}\), & \(\mathbf{x}_{l}\), & \(\left\|\mathbf{x}_{i}\right\|_{2}\) & \(\varepsilon^{x}\), \(\varepsilon^{y}\), \(\rho^{x}\), \(\varepsilon^{y}\): & 5 & MAE \(+\) L1 reg. \\ design & \(\uparrow\) & coefficients, & level & \(\mathbf{I}_{u/l_{1}}\), & \(\mathbf{I}_{u/l_{1}}\), & \(n_{u}\)\(=\)\(5,n_{w}\)\(=\)\(32\), \(n_{i}\)\(=\)\(32\) & & epoch=\(500\), \\ (2D) & 194 – 660 & \(C_{L}\
outcomes. The models are sensitive to the learning rate. We set an initial learning rate at 5e-4 or 1e-3, then reducing it by half on plateauing. The training batch size is set to be 16, 32, or 64. During the mini-batch training, all samples (graphs) within the batch are merged into a single disconnected global graph for performing forward and backward passes.

Finally, the trained models are evaluated by the relative L2 norm error,

\[\epsilon_{R}=\frac{\left\|\mathbf{v}_{target}-\mathbf{v}_{predict}\right\|_{2}}{\left\| \mathbf{v}_{target}\right\|_{2}}\times\mathbf{100\%}. \tag{5}\]

For the node level target models, we first compute their \(\epsilon_{R}\) based on all node-wise predictions in each individual graph. We then report the statistics, i.e., median (min., max.), of all the graphs' \(\epsilon_{R}\). For the graph level target models, we report a single \(\epsilon_{R}\) value computed from all the predictions. The performance of models on respective data sets are reported in Table II.

### _Modelling results on feature design problem_

The residual stress model achieves a performance of about 9% median \(\epsilon_{R}\) on both train and test set 1. From the visualization results presented in Figure 2(a), we observe a good agreement between the model prediction and the ground truth solution for the geometries in test set 1. This shows that the GNN-based surrogate model can adequately learn from the training geometries and predict the quantitative outcome on the test geometries with a similar accuracy. Although not previously seen during the training process, the geometries in test set 1 and train are homogenous, since they are designed from the same set of basic features with different sizes, orientations, and intersections.

Fig. 3: For both (a) & (b), the simulation solutions used as ground truth are displayed at the top rows. The model predictions and their absolute errors are displayed at the middle and bottom rows, respectively.
The prediction results of more challenging feature designs in test set 2 are visualized in Figure 2(b). For the 3 geometries comprising more complex combinations of the basic features, their \(e_{R}\)'s are 14.55%, 20.07%, and 28.49%--almost 2-3 times higher than the median \(\epsilon_{R}\) seen in train and test set 1. Note that the model has difficulty predicting the high residual stress around the intersection and joint areas, which could have great impact on the quality of product--the designer usually needs to pay attention to the high stress area. The original and optimized engineering bracket designs contain new features, such as triangle, hole, and arch, that are unseen by the model during the training. Their \(\epsilon_{R}\)'s are 17.14% and 34.82% respectively. Large prediction errors are found near the bottom holes. Moreover, the \(\epsilon_{R}\) for original and optimized like stem designs are 43.12% and 48.96%, respectively. These 2 designs are bulky, and they contain a lot of curvatures, which are significantly different from the training features in both size and shape.

Overall, the model performance is worse when being applied to out-of-distribution designs. Although it is possible to further reduce the prediction error on train and test set 1 by using a different set of model and training hyperparameters, these models usually do not demonstrate any accompanying performance improvement on test set 2. The performance of our GNN models is consistent to the prior study on a similar dataset while using the state-of-the-art U-Net model [7]. The U-Net uses a different (voxel) representation and is trained on a slightly larger set of geometries with data augmentation. The model gives prediction error between 14-29% on 4 out-of-distribution designs (which is the 1\({}^{\text{st}}\), 3\({}^{\text{st}}\)-5\({}^{\text{th}}\) geometry in our test set 2).

Although it is difficult to directly compare the computation time for different methods due to differences in implementation, our GNN model leads to a significant reduction in runtime as compared to the full order simulation, which requires approximately 1-4hours for simulating the residual stress distribution for each geometry [7], as compared to the instant prediction (\(<\) 1second) by the GNN model.

### _Modelling results on airfoil design problem_

The pressure model achieves a performance within the range of \(\epsilon_{R}\)=3.47% and 28.03% on the test set, which consists of 30 airfoil designs that are not seen by the model during the training, and there are 3 random freestream conditions for each design. Their median \(\epsilon_{R}\) is 8.71%. Figure 3(a) compares the prediction with the simulation ground truth, for 4 different airfoil designs and freestream conditions. The predicted pressure values are very accurate judging by the visual

Figure 4: (a) The pressure prediction on upper (U) and lower (L) surfaces of the airfoil is plotted against the simulation ground truth. The nested plot shows the airfoil design and the freestream condition that are used as the model input. (b-c) The model prediction for \(C_{D}\) and \(C_{L}\) are plotted against the simulation ground truth. The diagonal dotted line indicates a perfect fit to the simulation ground truth.

Figure 5: The relative L2 norm error, \(\epsilon_{R}\) for pressure, \(C_{D}\) and \(C_{L}\) between model prediction and simulation ground truth on test data. Predictions with high \(\epsilon_{R}\), as indicated by reds and dark rods, are more commonly found near to the edge of the distribution in the 2D freestream condition input space.
inspection, even though the quantitative \(\epsilon_{R}\) is around 9%. Note that the high-pressure profile may occur at either the upper or lower surface, depending on the freestream condition \(\nu_{o}\). Figure 4b-c present the prediction results for \(\mathcal{C}_{D}\) and \(\mathcal{C}_{L}\). Their \(\epsilon_{R}\)'s are 14.23% and 8.41% respectively on the same test set.

It is observed that there are some designs, e.g., a very thick airfoil, and certain freestream conditions where the model prediction has a large \(\epsilon_{R}\). These samples generally have a more complex pressure profile. Even so, the pressure model can still capture the general trend, as shown in the bottom panel of Figure 4a. Although it is difficult to verify whether those with large \(\epsilon_{R}\) are indeed entirely due to out-of-distribution designs, we note that the worse predictions are usually clustered around the edge of sample distribution in the freestream condition space, as shown in Figure 5.

Overall, we obtain satisfying results on surface pressure, \(\mathcal{C}_{D}\) and \(\mathcal{C}_{L}\) models. We believe that the issue of out-of-distribution designs will not be too severe, given that these models are trained on a large set of airfoil designs and freestream conditions. These models have the potential to replace simulations for the design of new airfoil on a range of operating conditions. The OpenFOAM simulation time for computing the flow solution on a new airfoil in current dataset takes around 40-70seconds.

## V Conclusion

In this work, we apply GNN model for the surrogate modelling task of predicting the functional performance of different geometrical designs in engineering applications. The GNN allows us to directly train the model on unstructured mesh simulations using graph representation, which can more precisely describe the change in geometry with complex features and predict the consequences.

We present the modelling results on 2 applications: feature design in the domain of additive engineering and airfoil design in the domain of aerodynamics. For the feature design task, the model is designed to predict the residual stress--which has great impact on the quality of the printed product--for features with different shape, sizes, orientations, and intersections. For the airfoil design task, we develop models for predicting the pressure on airfoil surface, as well as key performance indicators in drag and lift coefficients. These models have been shown to provide good accuracy on a separate set of test geometries which share similar features and characteristics to the training geometries. However, we have found that generalization remains a key issue for the residual stress model where the number of training geometries is smaller. The prediction error on out-of-distribution designs is almost 2-5 times higher than the geometries used for training. Further tests on different scenarios and application domains, ideally with an extended dataset, to confirm the potential of GNN-based surrogate model can be beneficial.

The strength of the surrogate models is their prediction speed. They can make predictions for new geometry designs almost instantly, as compared to O(hour) for the high-fidelity simulations typically used in most engineering applications. Our future work includes integrating the GNN-based surrogate models into generative design process in conjunction with other CI techniques such as generative adversarial network.

## Acknowledgment

This research is supported by A*STAR under the AME Enzymatic Programme: "Explainable Physics-based AI for Engineering Modelling and Design (ePAI)" Award No. A20H5b0142 and the IAF-PP project: "Industrial Digital Design and Additive Manufacturing Workflows (IDDAMW)" Award No. A19E10097.
