# Quasi-Monolithic Graph Neural Network for
Fluid-Structure Interaction

Rui Gao, Rajeev K. Jaiman
Department of Mechanical Engineering
University of British Columbia
Vancouver, Canada
{garrygao<EMAIL_ADDRESS>

###### Abstract

Using convolutional neural networks, deep learning-based reduced-order models
have demonstrated great potential in accelerating the simulations of coupled
fluid-structure systems for downstream optimization and control tasks.
However, these networks have to operate on a uniform Cartesian grid due to the
inherent restriction of convolutions, leading to difficulties in extracting
fine physical details along a fluid-structure interface without excessive
computational burden. In this work, we present a quasi-monolithic graph neural
network framework for the reduced-order modelling of fluid-structure
interaction systems. With the aid of an arbitrary Lagrangian-Eulerian
formulation, the mesh and fluid states are evolved temporally with two sub-
networks. The movement of the mesh is reduced to the evolution of several
coefficients via proper orthogonal decomposition, and these coefficients are
propagated through time via a multi-layer perceptron. A graph neural network
is employed to predict the evolution of the fluid state based on the state of
the whole system. The structural state is implicitly modelled by the movement
of the mesh on the fluid-structure boundary; hence it makes the proposed data-
driven methodology quasi-monolithic. The effectiveness of the proposed quasi-
monolithic graph neural network architecture is assessed on a prototypical
fluid-structure system of the flow around an elastically-mounted cylinder. We
use the full-order flow snapshots and displacements as target physical data to
learn and infer coupled fluid-structure dynamics. The proposed framework
tracks the interface description and provides the state predictions during
roll-out with acceptable accuracy. We also directly extract the lift and drag
forces from the predicted fluid and mesh states, in contrast to existing
convolution-based architectures. The proposed reduced-order model via graph
neural network has implications for the development of physics-based digital
twins concerning moving boundaries and fluid-structure interactions.

##  1 Introduction

Accurate and efficient predictions and control of the spatial-temporal
dynamics of fluid-structure systems are essential in various engineering
disciplines. The coupling between fluid and structure can lead to complex
dynamical effects such as vortex-induced vibrations and flutter/galloping [1].
Owing to the highly nonlinear and multiscale characteristics of fluid-
structure coupling, state-of-the-art computational fluid dynamics (CFD) and
finite element analysis (FEA) tools based on the solution of partial
differential equations are considered for engineering analysis and design.
While high-fidelity CFD/FEA simulations can provide accurate prediction, the
high computational cost of these simulations limits their applications for
downstream design optimization and control tasks. High-fidelity CFD/FEA
simulations can serve the role of a data generator for data-driven
predictions, with efficient design iterations or active control of coupled
fluid-structure systems. This work is motivated by the need to make coupled
fluid-structure simulations efficient for the digital twin technology, whereby
multi-query analysis, optimization and control are required.

This limitation of high-fidelity modelling based on partial differential
equations has inspired the development of the so-called data-driven reduced-
order modeling (ROM) techniques. By constructing low-dimensional models, ROM
techniques have the potential to address the limitations of high-fidelity
CFD/FEA models for efficient multi-query analysis, optimization and control
tasks. Such data-driven models can be used in an offline-online manner,
whereby the models can be trained to learn a low-dimensional representation of
the system from the high-dimensional physical data in the offline stage, and
provide efficient predictions during the online stage. Numerous data-driven
techniques have been developed for low-dimensional modelling and predictive
offline-online applications. Popular methods like proper orthogonal
decomposition (POD) [2, 3] and dynamic mode decomposition (DMD) [4], as well
as their variants (e.g., [5, 6, 7]), are usually based on the projection into
a linear subspace. However, these methods encounter difficulty when applied to
scenarios with high Reynolds numbers and convection-dominated problems,
whereas one needs a significantly large number of linear subspaces to achieve
a satisfactory approximation.

In the last decade, deep neural networks have been explored as alternatives to
the aforementioned techniques. Autoencoders, as a non-linear extension of POD,
have been shown to have significant advantages over POD in the compression of
data [8]. Combined with convolution, which supplies the inductive bias of
locality and translational equivariance, the convolutional encoder-propagator-
decoder architecture have been adopted in many applications, including flow
over fixed bodies in both 2D [9, 10, 11, 12] and 3D [13] scenarios, fluid-
structure interactions [14, 15], under-water noise propagation [16], and many
more. Leveraging the computational power of modern graphical processing units
(GPU), the computational times needed for these convolutional neural networks
are usually reduced to the level of milliseconds per step, which are orders of
magnitude faster than the traditional full-order CFD simulations.

While convolutional neural networks have achieved numerous successes, they
have inherent limitations and pose challenges in dealing with fluid-structure
boundaries. As convolutions have to be performed on a uniform Cartesian grid,
the resolution of the grid on different regions within the simulation domain
has to maintain similar resolutions. Considering the flow past a bluff body as
an example, the region near the body’s surface is resolved by a grid of the
same density as the far field. This means that one must either introduce a
significant number of grids in a far-field region that does not feature much
physics, or tolerate a relatively low resolution near the body’s surface. The
former choice of maintaining a dense grid throughout the bulk and interface
requires substantially more computational resources. As a result, the balance
usually aligns with the latter, inevitably leading to the loss of fine
physical details and the difficulty in extracting important physical
statistics like lift and drag coefficients. Despite many efforts to mitigate
this issue, by forming a hybrid model with traditional machine learning
techniques such as proper orthogonal decomposition [17, 18], interpolation &
projection schemes [14], or special mesh design [19], a completely
satisfactory solution is not yet available.

Graph neural networks, recently introduced as a geometric deep learning
framework, have the ability to address this difficulty, although existing
applications are mostly restricted to fluid systems with fixed boundaries. As
its name suggests, graph neural networks operate on graphs, which can be
intuitively converted from any mesh, meaning that one can control the
resolution of the different parts of the simulation domain using the
strategies adopted in CFD. As a result, GNNs can maintain a significantly
better resolution for important regions within the domain with the same mesh
size compared with CNNs. With this clear advantage over CNN, graph neural
networks are recently introduced for modeling fluid flow. Applications include
flow around fixed bodies like cylinder or airfoil [20, 21], reacting flows
[22], flow field super-resolution [23], flow field completion [24], etc.
Additional techniques and designs like on-the-fly graph adaptation [20],
multi-graph with different level of fineness [21, 25], rotational equivariance
[21], quadrature integration-based loss [26], and polynomial processors [22]
are also combined with the graph neural networks to further boost their
performance.

While these works have demonstrated the potential of graph neural networks for
various fluid flow applications, no study exists for the application to fully-
coupled fluid-structure systems. Notably, Pfaff et al. [20] simulated the
dynamics of flags in their work, but the graph is limited to the flag itself,
and the fluid flow surrounding the flag is not simulated. Li et al. [27]
simulated the flow flushing around a free-moving rigid box within a container.
However, both the flow and the box are discretized into large particles, and
therefore losing the fine physical details. Most existing GNN-based works on
fluid flow modelling, to the best of the authors’ knowledge, only focus on
applications with fixed, rigid body and/or domain boundaries.

In this work, we target to fill this gap in the literature by modelling fluid-
structure systems via a graph neural network. Adopting the arbitrary
Lagrangian-Eulerian formulation, we propose a quasi-monolithic GNN framework
for modelling fluid-structure systems. A GNN sub-network is employed to
predict the evolution of the fluid state based on the state of the whole
fluid-structure system. Mesh and solid movements are first projected to a
lower dimension by POD, and then propagated through time by a multi-layer
perceptron sub-network. The solid-state is implicitly modelled by the movement
of the mesh on the fluid-solid boundary. The framework is applied to a
prototypical fluid-structure problem of an elastically-mounted cylinder case,
and is demonstrated to generate accurate roll-out predictions. More
importantly, accurate lift and drag force predictions can be directly
extracted from the predicted system state by simply integrating the Cauchy
stress tensor at the surface of the moving solid body, which is difficult for
existing convolution-based frameworks.

This article is organized as follows. In Section 2, we describe the full-order
FSI system and the individual components such as proper orthogonal
decomposition, multi-layer perceptron and graph neural net for our proposed
reduced-order methodology. Section 3 presents a novel quasi-monolithic GNN
framework by combining the components. Detailed setup and the results of the
experiments on a prototypical fluid-structure interaction problem of an
elastically-mounted cylinder system undergoing vortex-induced vibration are
covered in Section 4. We conclude the work in Section 5.

##  2 Computational methodology: Review

Before presenting our quasi-monolithic methodology, we first provide a review
of full-order and reduced-order representations of the fluid-structure system.
Subsequently, we will introduce the individual components that will be
assembled into a complete framework in Sec. 3.

###  2.1 Full-order system

We briefly summarize the full-order system comprising the Eulerian fluid and
the Lagrangian structure, together with the traction and velocity continuity
conditions at the fluid-structure interface. Under the arbitrary Lagragian-
Eulerian (ALE) framework, the coupled system between isothermal incompressible
viscous fluid flow and an irrotational elastically-mounted rigid body can be
modelled as

$\rho^{f}\frac{\partial\boldsymbol{u}^{f}}{\partial
t}+\rho^{f}(\boldsymbol{u}^{f}-\boldsymbol{w})\cdot\nabla\boldsymbol{u}^{f}=\nabla\cdot\boldsymbol{\sigma}^{f}+\boldsymbol{b}^{f}\quad\text{on}\quad\boldsymbol{\Omega}^{f},$
(1a)
$\nabla\cdot\boldsymbol{u}^{f}=0\quad\text{on}\quad\boldsymbol{\Omega}^{f},$
(1b) $\boldsymbol{m}^{s}\frac{\partial\boldsymbol{u}^{s}}{\partial
t}+\boldsymbol{c}^{s}\boldsymbol{u}^{s}+\boldsymbol{k}^{s}(\boldsymbol{\varphi}^{s}(t)-\boldsymbol{\varphi}^{s}(0))=\boldsymbol{F}^{s}+\boldsymbol{b}^{s}\quad\text{on}\quad\boldsymbol{\Omega}^{s},$
(1c)

with boundary conditions on the surface of the rigid body

$\boldsymbol{u}^{f}(t)=\boldsymbol{u}^{s}(t)\quad\text{on}\quad\boldsymbol{\Gamma}^{fs},$
(2a)
$\boldsymbol{F}^{s}(t)=-\int_{\boldsymbol{\Gamma}^{fs}(t)}{\boldsymbol{\sigma}^{f}\cdot\boldsymbol{n}d\boldsymbol{\Gamma}}.$
(2b)

The superscripts $(\cdot)^{f}$ and $(\cdot)^{s}$ denote the state parameters
for fluid and solid respectively. The fluid, mesh velocity and body force are
denoted by $\boldsymbol{u}^{f}$, $\boldsymbol{w}$ and $\boldsymbol{b}^{f}$
respectively within the fluid domain $\boldsymbol{\Omega}^{f}$, while
$\boldsymbol{\varphi}^{s}$, $\boldsymbol{u}^{s}$, $\boldsymbol{F}^{s}$ and
$\boldsymbol{b}^{s}$ denote the displacement, the velocity, the traction force
and the body force for the solid body $\boldsymbol{\Omega}^{s}$, respectively.
The mass, damping, and stiffness matrices for the rigid body are denoted by
$\boldsymbol{m}^{s}$, $\boldsymbol{c}^{s}$ and $\boldsymbol{k}^{s}$
respectively, while $\rho^{f}$ is the fluid density. At the solid-fluid
interface $\boldsymbol{\Gamma}^{fs}$, $\boldsymbol{n}$ denotes the unit
outward normal vector. Assuming Newtonian fluid, the Cauchy stress tensor is
written as

$\boldsymbol{\sigma}^{f}=-p^{f}\boldsymbol{I}+\mu^{f}(\nabla\boldsymbol{u}^{f}+(\nabla\boldsymbol{u}^{f})^{T}),$
(3)

in which $p^{f}$ denotes the pressure, and $\mu^{f}$ denotes the viscosity of
the fluid. We solve the coupled differential equations using the
Petrov–Galerkin finite element and the semi-discrete time stepping [1]. For
the stable and robust coupling of the fluid–structure interaction, we utilize
the nonlinear partitioned staggered procedure [28].

###  2.2 Reduced-order model

With regard to the reduced-order modelling, the governing Eqs. (1) and (2) can
be re-written in the abstract dynamical form

$\frac{d\boldsymbol{q}}{dt}=\tilde{\boldsymbol{F}}(\boldsymbol{q}),$ (4)

where $\boldsymbol{q}$ denotes the state of the FSI system. In the present
study, the state of the system includes the fluid state, which comprises the
velocity $\boldsymbol{u}^{f}=(u_{x}^{f},u_{y}^{f})$ and the pressure $p^{f}$,
as well as the solid state, which involves two degrees of freedom for the
translational motion of the rigid body. In Eq. (4), the right hand side term
$\tilde{\boldsymbol{F}}$ represents an underlying dynamic model for the fluid-
structure system described in Eqs. 1-3. If we discretize the system with a
certain fixed time step $\delta t$, then Eq. 4 can be rewritten following
forward Euler time integration:

$\boldsymbol{q}(t_{n+1})-\boldsymbol{q}(t_{n})=\delta\boldsymbol{q}=\boldsymbol{F}(\boldsymbol{q}(t_{n})),$
(5)

in which $\boldsymbol{q}(t_{n})$ denotes the system state at time step
$t_{n}$, and $\boldsymbol{F}$ is the full-order system state update function
in the discretized time domain.

With reduced-order modelling, one aims to construct an approximation
$\hat{\boldsymbol{F}}$ to the full-order function $\boldsymbol{F}$

$\boldsymbol{F}(\boldsymbol{q}(t_{n}))=\boldsymbol{q}(t_{n+1})-\boldsymbol{q}(t_{n})\approx\delta\hat{\boldsymbol{q}}=\hat{\boldsymbol{F}}(\boldsymbol{q}(t_{n}),\theta),$
(6)

in which $\theta$ denotes the parameters of the reduced-order model
$\hat{\boldsymbol{F}}$. These parameters are calculated, tuned or trained to
minimize the approximation error between the reduced-order model
$\hat{\boldsymbol{F}}$ and the full-order model $\boldsymbol{F}$

$\theta=\arg_{\theta}\min(E[\boldsymbol{F}-\hat{\boldsymbol{F}}(\theta)]),$
(7)

in which $E$ denotes the expectation operator, which can vary for different
approaches such as projection-based model reduction, system identification or
machine learning. For example, proper orthogonal decomposition as a
projection-based model reduction minimizes the mean-squared error, i.e.,
$E[\boldsymbol{F}-\hat{\boldsymbol{F}}(\theta)]=\|\boldsymbol{F}-\hat{\boldsymbol{F}}(\theta)\|_{2}$,
where $\|\cdot\|_{2}$ denotes an Euclidean norm. In this work, we construct a
fully data-driven reduced-order model $\hat{\boldsymbol{F}}$ with three
components: proper orthogonal decomposition, multi-layer perceptron, and graph
neural network. These building blocks will be introduced in the remaining
parts of this section, and assembled into a complete framework in Sec. 3.

###  2.3 Proper orthogonal decomposition

As will be discussed later in Sec. 3.1, we choose to use the arbitrary
Lagrangian-Eulerian formulation for the proposed quasi-monolithic GNN
framework, which features mesh movements over time. This movement is
inherently low-rank, and therefore suitable for a linear projection-based
reduction. We employ Proper orthogonal decomposition (POD) for this purpose.
Particularly, we use the spatial variant [3] of POD. Given a data matrix
$\boldsymbol{X}$ in which each column contains the mesh displacements at a
certain time step, a singular value decomposition can be performed,

$\boldsymbol{X}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{T}.$ (8)

Each column of the left unitary matrix $\boldsymbol{U}$ is called a mode.
Suppose the singular values on the diagonal of matrix $\boldsymbol{\Sigma}$
are in decreasing order. In that case, the left-most vectors within the matrix
$\boldsymbol{U}$ are the modes capturing the most variance in the data. For a
low-rank system, a reduced-order model can be constructed by projecting the
data onto the first $l$ modes,

$\boldsymbol{C}^{(1:l)}=(\boldsymbol{U}^{(1:l)})^{T}\boldsymbol{X},$ (9a)
$\hat{\boldsymbol{X}}=\boldsymbol{U}^{(1:l)}\boldsymbol{C}^{(1:l)},$ (9b)

where $\boldsymbol{U}^{(1:l)}$ denotes the matrix $\boldsymbol{U}$ truncated
to its left-most $l$ columns (i.e., retaining the first $l$ modes), the
coefficient matrix $\boldsymbol{C}^{(1:l)}$ corresponds to the projected
weight of the data on each of the modes, and $\hat{\boldsymbol{X}}$ denotes
the reduced-order approximation of $\boldsymbol{X}$ using the first few modes.
The temporal evolution of the system is then reduced using the evolution of
the mode $1$ through $l$ coefficients over time.

###  2.4 Multi-layer perceptron

Multi-layer perceptron (MLP), or alternatively called feedforward neural
network, is one of the simplest modern neural networks. For an input vector
$\boldsymbol{z}$, a multi-layer perceptron with $n_{l}-1$ hidden layers can be
written as

$f(\boldsymbol{z})=f_{n_{l}}\circ f_{n_{l}-1}\circ f_{n_{l}-2}\circ\cdots\circ
f_{2}\circ f_{1}(\boldsymbol{z}),$ (10)

in which $\circ$ denotes function composition. Each layer $f_{i}$ is defined
as

$f_{i}(\boldsymbol{z})=\begin{cases}\sigma_{i}(\boldsymbol{W}_{i}\boldsymbol{z}+\boldsymbol{b}_{i})&i=1,2,\ldots,n_{l}-1\\\
\boldsymbol{W}_{i}\boldsymbol{z}+\boldsymbol{b}_{i}&i=n_{l}\\\ \end{cases}$
(11)

where $\boldsymbol{W}_{i}\boldsymbol{z}+\boldsymbol{b}_{i}$ is a linear
transformation of the input vector $\boldsymbol{z}$. The weight matrix
$\boldsymbol{W}_{i}$ does not have to be a square matrix, i.e., the linear
transformation can project the input into a higher or lower dimension matrix.
The number of rows in the weight matrix $\boldsymbol{W}_{i}$ is called the
layer width, and the non-linear function $\sigma_{i}$ is usually called the
activation function. In this work, multi-layer perceptrons are employed to
temporally evolve the POD coefficients of the mesh movements, and are also
used as the encoders, decoders, as well as processors for the graph neural
network, introduced next in Sec. 2.5.

###  2.5 Graph neural network

Graph neural networks, as discussed in the introduction, have demonstrated
their advantage over convolution-based neural networks in modeling fluid
systems, and therefore also utilized in this work. Specifically, we adopt the
variant of the message-passing graph neural network proposed by Battagalia et
al. [29]. This network architecture has been applied to various fluid flow
cases [20, 21, 22] and performed well. The network includes three steps:
Encoding, message-passing, and decoding. Given a graph converted from a mesh
following the illustration in Fig. 1, along with the feature vectors attached
to its nodes $\mathcal{V}$ and edges $\mathcal{E}$, the encoder step can be
written as

$\text{v}_{i}\leftarrow g^{\text{v}}(\text{v}_{i}),$ (12a)
$\text{e}_{ij}\leftarrow g^{\text{e}}(\text{e}_{ij}),$ (12b)

in which $\text{v}_{i}$ is the feature vector attached to node $i$,
$\text{e}_{ij}$ is the feature vector attached to the directed edge pointing
from node $i$ to node $j$, and $g^{\text{v}}$ and $g^{\text{e}}$ are the node
and edge encoders respectively. The left arrow $\leftarrow$ updates the
parameter on its left-hand side by its right-hand side value.

Figure 1: An illustration to transform an unstructured mesh to the
corresponding graph: (a) triangular mesh near the solid-fluid boundary, (b)
equivalent graph representation. Each vertex in the mesh is converted to a
node, while each link between two vertices in the mesh is converted to two
directed edges. The direction of each edge is marked by an arrow. System
states attached to the mesh are also converted and re-attached to the nodes
and edges as their respective feature vectors.

The encoded nodes and edges then proceed through a series of message-passing
layers, with each layer consisting of an edge update step and a node update
step, plotted in Fig. 2. Using the information from each edge itself and the
two nodes it connects, the edge update step can be written as

$\text{e}_{ij}^{\prime}=\phi^{\text{e}}(\text{e}_{ij},\text{v}_{i},\text{v}_{j}).$
(13a) In the subsequent node update step, each node is updated using the
information from the node itself and the aggregated information from all the
edges pointing to it:
$\text{v}_{i}^{\prime}=\phi^{\text{v}}(\text{v}_{i},\text{AGG}_{j}{\text{e}_{ji}^{\prime}}).$
(13b)

The aggregation function AGG is chosen to be the mean function in this work.
Note that the edge and node processors $\phi^{\text{e}}$ and $\phi^{\text{v}}$
are not necessarily shared for different message-passing layers. The outputs
at the end of the message-passing layers are decoded to generate the final
outputs. In this work, only outputs on nodes are needed, and therefore only a
node decoder $h^{\text{v}}$ is required:

$\text{v}_{i}^{\prime}\leftarrow h^{\text{v}}(\text{v}_{i}^{\prime}).$ (14)

With forward Euler time discretization described in Eq. 5, the output nodes
$\text{v}_{i}^{\prime}$ at each time step will be added to the corresponding
entries of the input to serve as the input of the next time step. Further
details will be presented in Sec. 3.2, after the discussion of the formulation
that governs the proposed quasi-monolithic GNN framework.

Figure 2: Schematic of the edge and node update steps within each message-
passing layer of graph neural network: (a) edge update step, and (b) node
update step. In the edge update step, the feature vector $\text{v}_{i}$
attached to node $i$, the feature vector $\text{v}_{j}$ attached to node $j$,
as well the feature vector $\text{e}_{ij}$ attached to the edge connecting the
two nodes, are processed by the edge processor $\phi^{\text{e}}$ to generate
the updated edge feature vector $\text{e}_{ij}^{\prime}$. In the node update
step, the previously updated feature vectors of all the edges connecting to
node $i$ are first aggregated through function AGG, and then processed along
with the node feature vector $\text{v}_{i}$ by the node processor
$\phi^{\text{v}}$ to generate the updated node feature vector
$\text{v}_{i}^{\prime}$. The parts of the graph shown in semi-transparent
colour are not used within the update steps for the edge feature
$\text{e}_{ij}$ in (a) and the node feature $\text{v}_{i}$ in (b).

##  3 Quasi-monolithic graph neural network

In this section, we construct the quasi-monolithic GNN framework for FSI
systems using the components introduced in Sec. 2. Under the arbitrary
Lagrangian-Eulerian formulation, the framework contains two sub-networks: a
GNN for the temporally propagating fluid state, and an MLP to evolve the mesh
state that is projected to a low dimension via POD. A schematic of the data
flow during temporal roll-out is shown in Fig. 3. In the remaining parts of
this section, we will first discuss the choice of formulation and then
describe the setups of the two sub-networks.

Figure 3: Schematic of the data flow within the proposed quasi-monolithic
graph neural network during prediction roll-out. The mesh state at time step
$t_{n}$ is reduced by POD, and then evolved by a multi-layer perceptron to
generate the mesh state at time step $t_{n+1}$. The fluid state is at time
step $t_{n+1}$ is predicted by a graph neural network, using the whole system
state at time step $t_{n}$ as input. The boundary information and the mesh
information are used by the graph neural network in every time step but are
not updated. Residual links are adopted following the forward Euler time
discretization described in Eq. 5.

###  3.1 Choice of formulation, node and edge features

Similar to full-order CFD simulations, the formulation used for the framework
has to be chosen first. Since the Reynolds number of the fluid flow is not
very low, a fully Lagrangian formulation is not preferable. What is left is
the choice between fully Eulerian and arbitrary Lagrangian-Eulerian
formulations. For fully-Eulerian formulation, as the solid body is moving
within the domain, on-the-fly graph adaptation would be preferable to help
maintain a relatively small mesh and graph size. However, numerical
experiments by Pfaff et al. [20] showed that the adaptive step for the graph
could cost significantly more computational time than running the neural
network itself for one step. We therefore choose to use the ALE formulation
for the present framework, which features mesh movement according to the
movement of the solid body, while the expensive graph adaptation is not
necessary.

With the formulation determined, the next step is to convert the mesh to a
graph, and re-attach the flow data on the mesh to the nodes and edges of the
converted graph. As mentioned in Sec. 2.5, we convert the mesh to graph
following the illustration in Fig. 1. The relative location information
between the vertices within the mesh are attached to the edges, following
Pfaff et al. [20],

$\text{e}_{ij}=[x_{i}-x_{j},y_{i}-y_{j},\sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}}],$
(15)

where $(x_{i},y_{i})$ denotes the spatial location of the vertex corresponding
to node $i$.

The displacement of the mesh from its initial location, however, has to be
attached to the nodes because it would be difficult to recover the correct
location of the nodes during prediction roll-out if this information is
attached to the edges. The fluid state parameters $u_{x},u_{y},p$, as well as
the border information $\gamma$ are also attached to the nodes. We drop the
superscript $(\cdot)^{f}$ in the remaining parts of the paper for brevity. The
full node feature vector is, therefore

$\text{v}_{i}=\left[u_{x,i},u_{y,i},p_{i},\delta x_{i},\delta
y_{i},\gamma_{i}\right],$ (16)

where $u_{x,i}$ and $u_{y,i}$ denote the two velocity components of the flow
corresponding to the node $i$, $p_{i}$ is the pressure, $\delta x_{i}$ and
$\delta y_{i}$ denote the displacement of the vertex corresponding to node $i$
from its location on the input mesh, and the border information $\gamma_{i}$
are encoded in a one-hot manner, which is similar to the work of Pfaff et al.
[20].

###  3.2 GNN for temporal propagation of flow state

We apply the GNN described in Sec. 2.5 for the temporal propagation of the
fluid state. Following Eq. 5, we use the network to predict the difference
between neighboring temporal steps. With the input node and edge features at
time step $t_{n}$ following Eq. 15 and 16, the output from the graph neural
network are

$\hat{\text{v}}_{i,out}\approx\left[u_{x,i},u_{y,i},p_{i}\right]_{t_{n+1}}-\left[u_{x,i},u_{y,i},p_{i}\right]_{t_{n}},$
(17)

where $\hat{\text{v}}_{i,out}$ denotes the output for node $i$. The fluid
state at time step $t_{n+1}$ is recovered during the temporal roll-out by
adding the neural network output to the corresponding entries of the node
feature at time step $t_{n}$.

###  3.3 POD-MLP for temporal propagation of mesh

We use a hybrid POD-MLP sub-network to evolve the mesh state. The mesh
displacements are first projected to a low dimension using POD, effectively
compressing the displacements of all the nodes over time into the evolution of
$l$ coefficients. Given the displacement of a mesh containing $n_{v}$ vertices
(nodes) over $n_{t}$ time steps, the displacements in the $x$ and $y$
directions are stacked to create a matrix $X$ of size $2n_{v}\times n_{t}$,

$\boldsymbol{X}=\begin{bmatrix}\delta x_{1,1}&\delta x_{1,2}&\cdots&\delta
x_{1,n_{t}}\\\ \delta x_{2,1}&\delta x_{2,2}&\cdots&\delta x_{2,n_{t}}\\\
\vdots&\vdots&\ddots&\vdots\\\ \delta x_{n_{v},1}&\delta
x_{n_{v},2}&\cdots&\delta x_{n_{v},n_{t}}\\\ \delta y_{1,1}&\delta
y_{1,2}&\cdots&\delta y_{1,n_{t}}\\\ \delta y_{2,1}&\delta
y_{2,2}&\cdots&\delta y_{2,n_{t}}\\\ \vdots&\vdots&\ddots&\vdots\\\ \delta
y_{n_{v},1}&\delta y_{n_{v},2}&\cdots&\delta y_{n_{v},n_{t}}\end{bmatrix},$
(18)

and then decomposed and truncated following Eqs. 8 and 9. An MLP is then used
for the temporal prediction of these coefficients.

######  Remark.

It is worth mentioning that this idea of a hybrid POD-neural network model has
also been used in our previous work [14], where an LSTM [30] network was used
along with POD to propagate the mesh location. In this work, we take a further
step from that idea, coupling the mesh and fluid prediction during roll-out.
We also simplified the network used from an LSTM network to a simple multi-
layer perceptron.

The inputs to the MLP are the history of the coefficients corresponding to the
first $l$ modes over the last $n_{h}$ time steps, while the output is the
difference between the coefficients at time step $t_{n+1}$ and time step
$t_{n}$:

$\text{MLP}\left(c_{t_{n}}^{(1:l)},c_{t_{n-1}}^{(1:l)},\ldots,c_{t_{n-n_{h}+1}}^{(1:l)}\right)=\hat{\delta
c}^{(1:l)}\approx c_{t_{n+1}}^{(1:l)}-c_{t_{n}}^{(1:l)},$ (19)

in which $\hat{\delta c}^{(1:l)}$ denotes the predicted change of coefficients
between time step $t_{n}$ and time step $t_{n+1}$. During temporal roll-out,
the POD coefficients at time step $t_{n+1}$ are recovered by adding the neural
network output to the POD coefficients at time step $t_{n}$. The coefficients
are then multiplied to the corresponding modes and summed to generate the mesh
displacements $(\boldsymbol{\delta{x}},\boldsymbol{\delta{y}})_{t_{n+1}}$ for
all the nodes. The solid movement is implicitly modelled by the movement of
the fluid-solid boundary nodes.

##  4 Experiments with elastically-mounted cylinder in flowing fluid

In this section, we evaluate the framework constructed in Sec. 3, focusing on
its capability to learn the temporal evolution of the FSI system and make
accurate temporal roll-out predictions. We consider a representative FSI
system of an elastically-mounted cylinder in a laminar flow. In the remaining
parts of this section, we will first introduce the generation of the ground
truth train and testing data sets, followed by the setup of the model and how
it is trained. We then present the metrics used for the evaluation, as well as
the results.

###  4.1 Ground truth data generation

The ground truth data used in this work is the flow around an elastically-
mounted cylinder simulated at laminar flow condition. The vibration of a
cylinder is governed by the four key non-dimensional parameters, namely mass-
ratio $\left(m^{*}\right)$, Reynolds number $\left(Re\right)$, reduced
velocity $\left(U_{r}\right)$ and critical damping ratio $(\zeta)$ defined as

$\displaystyle m^{*}=\frac{m^{s}}{m^{f}},\ \ \ \
Re=\frac{\rho^{f}U_{\infty}D}{\mu^{\mathrm{f}}},\ \ \ \
U_{r}=\frac{U_{\infty}}{f_{n}D},\ \ \ \
\zeta=\frac{c^{s}}{2\sqrt{k^{s}m^{s}}},$ (20)

where $m^{s}$ is the mass of the body, $c^{s}$ and $k^{s}$ are the damping and
stiffness coefficients, respectively for an equivalent mass-spring-damper
system of a vibrating structure, $U_{\infty}$ and $D$ denote the free-stream
speed and the diameter of cylinder, respectively. The natural frequency of the
body is given by $f_{n}=(1/2\pi)\sqrt{k^{s}/m^{s}}$, the mass of displaced
fluid by the structure is $m^{f}=\rho^{f}\pi D^{2}/4$ for a cylinder cross-
section. For the current work, we consider the mass ratio $m^{*}=10$, the
Reynolds number $Re=200$, the reduced velocity $U_{r}=5$ and the damping
coefficient $c^{s}=0$.

We solve the full-order system (Eqs. 1 and 2) using Petrov-Galerkin finite
elements and a semi-discrete time-stepping scheme [28]. The domain and mesh
for the ground-truth simulation are plotted in Fig. 4 and 5a. The inlet has a
uniform flow $u_{x}=U_{\infty}=1,u_{y}=0$, the outlet boundary $\Gamma_{out}$
is traction free, and a slip-wall boundary condition is implemented at top and
bottom boundaries $\Gamma_{top}$ and $\Gamma_{bottom}$. The distance between
the side walls is $30D$, and the distances to the upstream and the downstream
boundaries are $10D$ and $20D$ respectively. The simulation results are
validated with the reference values from Li et al. [31]. A total of 2499
continuous steps are sampled from a stable condition with a time step $\Delta
t=0.04$, which is equivalent to the ground truth time step of the full-order
simulation.

Figure 4: Schematic of an elastically-mounted cylinder for full-order
simulation. Figure 5: Mesh used for the ground truth simulation and the neural
network-based simulation: (a) CFD mesh, (b) graph neural network mesh. Both
meshes contain only quadrilateral elements with refinement near the surface of
the cylinder and in the wake region. The blue square in (b) marks out the part
of the domain being plotted in Fig. 7.

For the neural network, we use a graph transformed from a coarser mesh (shown
in Fig. 5b), which has 2204 nodes and 8728 directed edges. The ground truth
ALE mesh movement of the coarser mesh at each time step is calculated from the
ground truth displacement $\boldsymbol{\varphi}^{m}=(\boldsymbol{\delta
x},\boldsymbol{\delta y})$ of the cylinder by modeling the mesh as an elastic
material in equilibrium, which is given by

$\nabla\cdot\left((1+k_{m})[\nabla\boldsymbol{\varphi}^{m}+(\nabla\boldsymbol{\varphi}^{m})^{T}+(\nabla\cdot\boldsymbol{\varphi}^{m})\boldsymbol{I}]\right)=0,$
(21)

with boundary conditions

$\boldsymbol{\varphi}^{m}(t)=\varphi^{s}(t)-\varphi^{s}(0)\quad\text{on}\quad\boldsymbol{\Gamma}^{fs},$
(22a)
$\boldsymbol{\varphi}^{m}=0\quad\text{on}\quad\partial\boldsymbol{\Omega}^{f}\backslash\boldsymbol{\Gamma}^{fs},$
(22b)

where $\partial\boldsymbol{\Omega}^{f}(0)\backslash\boldsymbol{\Gamma}^{fs}$
denotes the boundary of the fluid domain excluding fluid-solid interface, and
$k_{m}$ is the local element level mesh stiffness parameter [1, 32] which we
choose as 0 here for simplicity. The fluid state parameters $u_{x},u_{y},p$
are then interpolated from the CFD mesh (Fig. 5a) to the coarser mesh (Fig.
5b) via the Clough-Tocher 2D interpolation scheme available in SciPy package
[33].

###  4.2 Train-test split and pre-processing

The interpolated fluid flow data are then attached to the nodes and edges
converted from the coarse mesh following Eqs. 15 and 16, and subsequently
split into training and testing data sets. Within all 2499 time steps, the
first 2061 time steps are used for training purposes, and the rest are used
for testing purposes.

For the fluid state prediction, we construct two sets of training data. In the
one-step prediction training set, the input nodes $\text{v}_{i}|_{t_{n}}$ and
edges $\text{e}_{ij}$ corresponds to the output nodes
$[u_{x,i},u_{y,i},p_{i}]_{t_{n+1}}-[u_{x,i},u_{y,i},p_{i}]_{t_{n}}$. In the
five-step roll-out prediction training set, the input nodes
$\text{v}_{i}|_{t_{n}}$ and edges $\text{e}_{ij}$ corresponds to output nodes
$[u_{x,i},u_{y,i},p_{i}]_{t_{n+5}}-[u_{x,i},u_{y,i},p_{i}]_{t_{n+4}}$. The
input and output nodes and edges are then normalized to zero mean and unit
variance, respectively.

For the mesh displacement prediction, we only use the one-step prediction data
for training, in which each input
$\left[c_{t_{n}}^{(1:l)},c_{t_{n-1}}^{(1:l)},\ldots,c_{t_{n-n_{h}+1}}^{(1:l)}\right]$
corresponds to an output $\left[c_{t_{n+1}}^{(1:l)}-c_{t_{n}}^{(1:l)}\right]$,
with the number of history steps chosen as $n_{h}=3$. The proper orthogonal
decomposition is performed on the centred mesh displacement data set over all
the 2061 time steps used for training. The coefficients are normalized to zero
mean and unit variance before constructing the inputs and outputs, which are
normalized again afterwards. We only use $l=2$ modes for the reduced-order
model, as the mesh displacement matrix is inherently rank-2.

The network is expected to generate accurate roll-out predictions, and
therefore the test set is designed such that the neural network roll-out can
start from any time step within the test set. At each time step within the
data set, the data available are the inputs to the fluid state prediction
$\text{v}_{i}|_{t_{n}}$ and $\text{e}_{ij}$, as well as the inputs to the mesh
displacement prediction
$\left[c_{t_{n}}^{(1:l)},c_{t_{n-1}}^{(1:l)},\ldots,c_{t_{n-n_{h}+1}}^{(1:l)}\right]$.
These data are normalized using the mean and variance statistics calculated
from the training data sets.

Table 1 summarizes the data sets.

Table 1: Training and testing data sets | Name | Inputs | Output | Total samples
---|---|---|---|---
Training sets | One step, fluid | $\text{v}_{i}|_{t_{n}},\text{e}_{ij}$ | $[u_{x,i},u_{y,i},p_{i}]_{t_{n+1}}-[u_{x,i},u_{y,i},p_{i}]_{t_{n}}$ | 2048
5 steps, fluid | $[u_{x,i},u_{y,i},p_{i}]_{t_{n+5}}-[u_{x,i},u_{y,i},p_{i}]_{t_{n+4}}$
One step, mesh | $[c^{(1:2)}_{t_{n}},c^{(1:2)}_{t_{n-1}},c^{(1:2)}_{t_{n-2}}]$ | $[c^{(1:2)}_{t_{n+1}}]-[c^{(1:2)}_{t_{n}}]$
Testing sets | Roll-out, fluid | $\text{v}_{i}|_{t_{n}},\text{e}_{ij}$ | N/A | 435
Roll-out, mesh | $[c^{(1:2)}_{t_{n}},c^{(1:2)}_{t_{n-1}},c^{(1:2)}_{t_{n-2}}]$ | N/A
Table 2: Neural network training details. | Sub-network
---
(total params; optimizer)
Stage | Training set | Batch size | Epochs |  | Starting
---
learning rate
| Learning rate
---
decay per epoch
Fluid prediction (2.33 M; Adam) | 1 | One step, fluid | 4 | 6 | $10^{-6}$ | 2.152
2 | 50 | $10^{-4}$ | 1
3 | 100 | $10^{-4}$ | 0.955
4 | 50 | $10^{-6}$ | 1
5 | 5 steps, fluid | 2 | 50 | $10^{-4}$ | 0.912
| Mesh prediction
---
(267 K; AdamW)
1 | One step, mesh | 512 | 200 | $10^{-4}$ | 0.98
Table 3: Time required per prediction roll-out step Component | Time per step
---|---
Total |  $\approx 21.4$ ms
Fluid prediction sub-network |  $\approx 12.4$ ms
Mesh prediction sub-network |  $\approx 0.2$ ms
Other overhead |  $\approx 8.8$ ms

###  4.3 Model and training details

With the training and testing data sets available, we proceed to set up and
train the model. The implementation is completed with the PyTorch package
[34]. For the graph neural network used in the fluid state prediction, the
encoders, decoder, and the node and edge processors for each message-passing
layer are all selected to be MLPs with two hidden layers and Gaussian Error
Linear Unit activation function [35]. The layer widths of all these MLPs are
fixed at 128 following Pfaff et al. [20]. The weights and biases are
initialized with PyTorch default setting. We tested over different numbers of
passage-passing layers since we use a quadrilateral mesh rather than the
triangular mesh in most existing works, and it seems that 15 message-passing
layers are still a reasonable choice to achieve the balance between accuracy
and computational cost. Layer normalization [36] is utilized after all
encoders as well as all node and edge processors to facilitate convergence.
The graph message-passing is implemented using the gather-scatter scheme
following the practice in PyTorch Geometric package [37].

For the MLP used for the mesh displacement prediction, we use two hidden
layers (i.e., $n_{l}=3$) with a layer width of 512. A sinusoidal activation
function [38] and the corresponding initialization scheme is adopted. With
history size $n_{h}=3$, the input size is six, and the output size is 2.

The two sub-networks are trained separately on a single Nvidia RTX 3090 GPU
with a central processing unit (CPU) being AMD Ryzen 9 5900 @ 3 GHz $\times$
12 cores. The fluid state prediction sub-network is trained with an Adam
optimizer [39] for a total of 206 epochs, divided into five stages listed in
table 2. The mesh displacement prediction sub-network is trained for 200
epochs with batch size 512. An AdamW [40] optimizer is used with a learning
rate starting at $10^{-4}$ and an exponential learning rate decay of $0.98$
per epoch. For both networks, we use a mean-square error loss normalized on
each batch as given by:

$\mathcal{L}=\frac{\|\psi-\hat{\psi}\|_{2}^{2}}{\|\psi\|_{2}^{2}},$ (23)

where $\psi$ denotes the ground truth output for the batch, and $\hat{\psi}$
denotes the predicted output from the neural network. Note that this loss is
not the normalized mean square loss, as we sum over the numerator and
denominator before the division under the consideration of numerical
stability.

Figure 6: Predicted vs. ground truth mode coefficients for the mesh
displacement during prediction roll-out over 400 steps, starting from the
first time step within the test data set.

###  4.4 Evaluation metrics

After training, the framework is evaluated on the testing data set to evaluate
its performance in producing temporal roll-out predictions. We consider
several metrics to facilitate the evaluation. Qualitatively, we plot out the
predicted pressure fields at different time steps during the prediction roll-
out, and compare them with the ground truth. Quantitatively, we use the
coefficient of determination $R^{2}$ to evaluate the accuracy of predictions
for the system states, which can be calculated given the system state
parameter $\boldsymbol{s}$ and its prediction $\hat{\boldsymbol{s}}$

$R^{2}=1-\frac{\|\boldsymbol{s}-\hat{\boldsymbol{s}}\|_{2}^{2}}{\|\boldsymbol{s}-\bar{\boldsymbol{s}}\|_{2}^{2}},$
(24)

where $\bar{\cdot}$ denotes the calculation of mean. A coefficient of
determination close to 1 means that the error of the prediction is low. In
addition, we also extract the lift coefficient $C_{l}$ and drag coefficient
$C_{d}$ on the cylinder from the predicted system state. The calculation is
performed by integrating the Cauchy stress tensor $\boldsymbol{\sigma}^{f}$
(defined in Eq. 3) from the first layer of element away from the cylinder,

$C_{l}=\frac{1}{\frac{1}{2}\rho^{f}(U_{\infty})^{2}D}\int_{\boldsymbol{\Gamma}^{fs}}(\boldsymbol{\sigma}^{f}\cdot\boldsymbol{n})\cdot\boldsymbol{n}_{y}d\boldsymbol{\Gamma},$
(25a)
$C_{d}=\frac{1}{\frac{1}{2}\rho^{f}(U_{\infty})^{2}D}\int_{\boldsymbol{\Gamma}^{fs}}(\boldsymbol{\sigma}^{f}\cdot\boldsymbol{n})\cdot\boldsymbol{n}_{x}d\boldsymbol{\Gamma},$
(25b)

where $U_{\infty}=1$ is the inlet velocity, and the characteristic length $D$
is the diameter of the cylinder in this work. It should be noted that this way
of calculation is different from our previous work [14], as only pressure lift
and drag coefficients were considered in that work, but the contribution from
viscous forces were not calculated. It is also different from works like [15],
as they did not calculate the lift and drag coefficients from the predicted
flow and mesh states, but rather model and evolve them separately.

###  4.5 Results

In this subsection, we present the results when evaluating the framework on
the test data set, using the metrics listed in Sec. 4.4. Table 3 summarizes
the approximate computational time needed for each prediction roll-out step
when running on an Nvidia RTX 2060 GPU with AMD Ryzen 7 4800H CPU @ 2.90 GHz
$\times$ 8 cores. Approximately 2.1 GB GPU memory is required for the roll-out
during the tests, including the model parameters, the buffers and other
overhead.

Figure 7: Predicted vs. ground truth normalized pressure field
$p^{*}=2p/\rho^{f}U_{\infty}^{2}$ near the cylinder, at 100, 200, and 300
steps during prediction roll-out. The region plotted is marked out by a blue
square in Fig. 5b.

The roll-out prediction of mesh POD coefficients over 400 time steps using the
mesh prediction sub-network, starting from the first time step within the test
data set (step 2064 within all the 2499 time steps), is plotted in Fig. 6. It
is clear that the mesh prediction sub-network produces stable and accurate
predictions of mesh displacements over time. We then assemble it with the
fluid prediction sub-network for the roll-out predictions of the whole system
state, following the data flow shown in Fig. 3.

Figure 8: Coefficient of determinations of the system state and its components
over 400 prediction roll-out steps.

#### System state prediction.

Starting from the first time step within the test data set, we predict the
system states over the next 300 time steps. Figure 7 shows the predicted
versus the ground truth normalized pressure $p^{*}=2p/\rho^{f}U_{\infty}^{2}$
near the cylinder at steps 100, 200 and 300. The plotted area is marked by a
blue square in Fig. 5b. For the purpose of plotting, the predicted field is
projected onto a $512\times 512$ uniform Cartesian grid using a linear
interpolation scheme available in the SciPy package [33]. The cylinder is
plotted by linking the predicted locations of its boundary nodes. The model is
able to produce visually accurate prediction results at time step 100. At 200
and 300 steps, the error starts to become more pronounced visually in the
regions below and above the cylinder, where the mesh (cf. Fig. 5) is
relatively sparse. This is a result of accumulated errors over the roll-out
time steps. Nevertheless, the predictions in the wake region, where the mesh
resolution is relatively high, remain visually accurate over all 300 predicted
steps. We also plot the coefficient of determination for the predicted system
state, as well as its components, over 400 time steps in Fig. 8.
Quantitatively, the coefficients of determination for the system state, as
well as its components, are very close to 1 in the first 150 steps, indicating
accurate predictions of the temporal evolution of the system. After that, the
error starts to accumulate slowly for the fluid state predictions, while mesh
state predictions remain accurate.

Figure 9: Predicted vs. ground truth lift and drag coefficients during
prediction roll-out over 300 steps, starting from different time steps within
the test set. From top to bottom, time steps include 1, 26, 51, 76, 101 within
the test set.

#### Lift and drag coefficients prediction.

In addition to calculating the coefficient of determinations of the fluid and
mesh state parameters, we also quantify the error of the predicted system
state by calculating the lift and drag coefficients over time following Eq.
25. It should be emphasized again that we are directly integrating the lift
and drag from the predicted system states, without any interpolation or
smoothing procedures. Starting from different time steps within the testing
data set, we plot out the predicted and ground truth lift and drag
coefficients over 300 time steps in Fig. 9. Accurate predictions of lift and
drag coefficients can be extracted from the predicted flow field accurate over
the first 200 steps for all these cases. After that, error starts to occur for
the predicted lift coefficient, while the predicted drag coefficient stays
accurate.

##  5 Conclusions

We presented a quasi-monolithic graph neural network-based framework for data-
driven reduced-order modeling of fluid-structure interaction systems. Adopting
the arbitrary Lagrangian-Eulerian formulation, the proposed framework relies
on a multi-layer perceptron for predicting POD coefficients to evolve mesh
displacements, and a graph neural network for predicting fluid states from the
current system state. The framework is applied to the flow around an
elastically mounted cylinder at Reynolds number $Re=200$. Starting from a
given system state, the model can produce accurate predictions of the system
state over 150 time steps, and the error only starts to accumulate slowly
after that in poorly-resolved regions. Accurate fluid and mesh state
predictions on a non-uniform mesh, refined near the surface of the cylinder,
enable the direct and accurate extraction of lift and drag forces over 200
time steps via the integration of the Cauchy stress tensor on the cylinder
surface, which is difficult with existing convolution-based works. It should
be noted that the framework proposed in this work applies to fluid-structure
interaction systems with rigid solid bodies. In our future work, we plan to
extend the application to fluid-structure interaction systems with deformable
elastic structures.

## Acknowledgments

The authors would like to acknowledge the Natural Sciences and Engineering
Research Council of Canada (NSERC) and Seaspan Shipyards for the funding. This
research was supported in part through computational resources and services
provided by Advanced Research Computing at the University of British Columbia.

## References

  * [1] Rajeev K Jaiman and Vaibhav Joshi.  Computational Mechanics of Fluid-Structure Interaction.  Springer, 2022.
  * [2] John L Lumley.  The structure of inhomogeneous turbulent flows.  In Atmospheric Turbulence & Wave Propagation, 1967.
  * [3] Lawrence Sirovich.  Turbulence and the dynamics of coherent structures. part I. coherent structures.  Quarterly of Applied Mathematics, 45(3):561–571, 1987.
  * [4] Peter J Schmid.  Dynamic mode decomposition of numerical and experimental data.  Journal of fluid mechanics, 656:5–28, 2010.
  * [5] Aaron Towne, Oliver T Schmidt, and Tim Colonius.  Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis.  Journal of Fluid Mechanics, 847:821–867, 2018.
  * [6] Oliver T Schmidt and Peter J Schmid.  A conditional space–time pod formalism for intermittent and rare events: example of acoustic bursts in turbulent jets.  Journal of Fluid Mechanics, 867, 2019.
  * [7] Hao Zhang, Clarence W Rowley, Eric A Deem, and Louis N Cattafesta.  Online dynamic mode decomposition for time-varying systems.  SIAM Journal on Applied Dynamical Systems, 18(3):1586–1609, 2019\.
  * [8] Geoffrey E Hinton and Ruslan R Salakhutdinov.  Reducing the dimensionality of data with neural networks.  science, 313(5786):504–507, 2006.
  * [9] Francisco J Gonzalez and Maciej Balajewicz.  Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems.  arXiv preprint arXiv:1808.01346, 2018.
  * [10] Nils Thuerey, Konstantin Weißenow, Lukas Prantl, and Xiangyu Hu.  Deep learning methods for reynolds-averaged navier–stokes simulations of airfoil flows.  AIAA Journal, 58(1):25–36, 2020.
  * [11] Sandeep R Bukka, Rachit Gupta, Allan R Magee, and Rajeev K Jaiman.  Assessment of unsteady flow predictions using hybrid deep learning based reduced-order models.  Physics of Fluids, 33(1):013601, 2021.
  * [12] Pranshu Pant, Ruchit Doshi, Pranav Bahl, and Amir Barati Farimani.  Deep learning for reduced order modelling and efficient temporal evolution of fluid simulations.  Physics of Fluids, 33(10):107101, 2021.
  * [13] Rachit Gupta and Rajeev K Jaiman.  Three-dimensional deep learning-based reduced order model for unsteady flow dynamics with variable reynolds number.  Physics of Fluids, 34(3):033612, 2022.
  * [14] Rachit Gupta and Rajeev K Jaiman.  A hybrid partitioned deep learning methodology for moving interface and fluid–structure interaction.  Computers & Fluids, 233:105239, 2022.
  * [15] Xinshuai Zhang, Tingwei Ji, Fangfang Xie, Changdong Zheng, and Yao Zheng.  Data-driven nonlinear reduced-order modeling of unsteady fluid–structure interactions.  Physics of Fluids, 34(5):053608, 2022.
  * [16] Wrik Mallik, Rajeev K Jaiman, and Jasmin Jelovica.  Predicting transmission loss in underwater acoustics using convolutional recurrent autoencoder network.  The Journal of the Acoustical Society of America, 152(3):1627–1638, 2022.
  * [17] Zheng Wang, Dunhui Xiao, Fangxin Fang, Rajesh Govindan, Christopher C Pain, and Yike Guo.  Model identification of reduced order fluid dynamics systems using deep learning.  International Journal for Numerical Methods in Fluids, 86(4):255–268, 2018.
  * [18] Sandeep B Reddy, Allan R Magee, Rajeev K Jaiman, J Liu, W Xu, A Choudhary, and A A Hussain.  Reduced order model for unsteady fluid flows via recurrent neural networks.  In International Conference on Offshore Mechanics and Arctic Engineering, volume 58776, page V002T08A007. American Society of Mechanical Engineers, 2019.
  * [19] Li-Wei Chen and Nils Thuerey.  Towards high-accuracy deep learning inference of compressible turbulent flows over aerofoils.  arXiv preprint arXiv:2109.02183, 2021.
  * [20] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter W Battaglia.  Learning mesh-based simulation with graph networks.  arXiv preprint arXiv:2010.03409, 2020.
  * [21] Mario Lino, Stathi Fotiadis, Anil A Bharath, and Chris D Cantwell.  Multi-scale rotation-equivariant graph neural networks for unsteady eulerian fluid dynamics.  Physics of Fluids, 34(8):087110, 2022.
  * [22] Jiayang Xu, Aniruddhe Pradhan, and Karthikeyan Duraisamy.  Conditionally parameterized, discretization-aware neural networks for mesh-based modeling of physical systems.  Advances in Neural Information Processing Systems, 34:1634–1645, 2021.
  * [23] Filipe De Avila Belbute-Peres, Thomas Economon, and Zico Kolter.  Combining differentiable pde solvers and graph neural networks for fluid flow prediction.  In international conference on machine learning, pages 2402–2411. PMLR, 2020.
  * [24] Xiaodong He, Yinan Wang, and Juan Li.  Flow completion network: Inferring the fluid dynamics from incomplete flow information using graph neural networks.  arXiv preprint arXiv:2205.04739, 2022.
  * [25] Leon Migus, Yuan Yin, Jocelyn Ahmed Mazari, and Patrick Gallinari.  Multi-scale physical representations for approximating pde solutions with graph neural operators.  arXiv preprint arXiv:2206.14687, 2022.
  * [26] Han Gao, Matthew J Zahr, and Jian-Xun Wang.  Physics-informed graph neural galerkin networks: A unified framework for solving pde-governed forward and inverse problems.  Computer Methods in Applied Mechanics and Engineering, 390:114502, 2022.
  * [27] Yunzhu Li, Jiajun Wu, Russ Tedrake, Joshua B Tenenbaum, and Antonio Torralba.  Learning particle dynamics for manipulating rigid bodies, deformable objects, and fluids.  arXiv preprint arXiv:1810.01566, 2018.
  * [28] Rajeev K Jaiman, M Z Guan, and Tharindu P Miyanawala.  Partitioned iterative and dynamic subgrid-scale methods for freely vibrating square-section structures at subcritical reynolds number.  Computers & Fluids, 133:68–89, 2016.
  * [29] Peter W Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanchez-Gonzalez, Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, et al.  Relational inductive biases, deep learning, and graph networks.  arXiv preprint arXiv:1806.01261, 2018.
  * [30] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [31] Zhong Li, Weigang Yao, Kun Yang, Rajeev K Jaiman, and Boo Cheong Khoo.  On the vortex-induced oscillations of a freely vibrating cylinder in the vicinity of a stationary plane wall.  Journal of Fluids and Structures, 65:495–526, 2016.
  * [32] Arif Masud and Thomas JR Hughes.  A space-time galerkin/least-squares finite element formulation of the navier-stokes equations for moving domain problems.  Computer methods in applied mechanics and engineering, 146(1-2):91–126, 1997.
  * [33] Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J van der Walt, Matthew Brett, Joshua Wilson, K Jarrod Millman, Nikolay Mayorov, Andrew R J Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E A Quintero, Charles R Harris, Anne M Archibald, Antônio H Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [34] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  Advances in neural information processing systems, 32, 2019.
  * [35] Dan Hendrycks and Kevin Gimpel.  Gaussian error linear units (gelus).  arXiv preprint arXiv:1606.08415, 2016.
  * [36] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton.  Layer normalization.  arXiv preprint arXiv:1607.06450, 2016.
  * [37] Matthias Fey and Jan E Lenssen.  Fast graph representation learning with PyTorch Geometric.  In ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019.
  * [38] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  Advances in Neural Information Processing Systems, 33:7462–7473, 2020.
  * [39] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [40] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  arXiv preprint arXiv:1711.05101, 2017.


# Newtonian Gravity on an N-Sphere

T. Curtright1 and H. Alshal2
${}^{1}\ $Department of Physics, University of Miami, Coral Gables, FL 33124
${}^{2}\ $Department of Physics, Santa Clara University, Santa Clara, CA 95053
curtright<EMAIL_ADDRESS>
###### Abstract

We consider some elementary features of Newtonian gravity, or electrostatics,
as defined on an N-sphere. In particular, we present and discuss “the shell
theorem” for this system.

## Introduction

The analogue of Newtonian gravity on a closed Riemannian manifold has been an
interesting subject for quite a long time, especially among mathematicians
(e.g. see [2] and the literature cited therein). The N-sphere presents a
tractable example that can be analyzed without approximations (e.g. see [4]
and the literature cited therein).

The physics of Newtonian gravity on the N-sphere is emphasized here.
Alternatively, one may think of the following as electrostatics on the
N-sphere. Either way, perhaps some insight and intuition can be gleaned from
this simple model.

Taking this point of view, we follow a path not very far from that trodden by
G. Green so long ago [5]. We identify the point-particle gravitational or
electrostatic potentials with a Green function of the system, and we define
the force on such particles in terms of the gradient of that potential, as is
usual in physics. We then consider the analogue of Newton’s shell theorem for
this simple system. We find a generalization of the shell theorem which is
simple to state as well as pleasing to our taste.

## N-Sphere Laplacian Green Function

On an N-sphere of radius $R$, acting on functions $f\left(\theta\right)$ with
only polar angle dependence:

$\nabla^{2}=\frac{1}{R^{2}}\left(\frac{1}{\left(\sin\theta\right)^{N-1}}\frac{\partial}{\partial\theta}\left(\left(\sin\theta\right)^{N-1}\frac{\partial}{\partial\theta}\right)\right)$
(1)

For small angles,
$\nabla^{2}\sim\frac{1}{R^{2}}\left(\frac{1}{\theta^{N-1}}\frac{\partial}{\partial\theta}\left(\theta^{N-1}\frac{\partial}{\partial\theta}\right)\right)=\frac{1}{R^{2}}\left(\frac{\partial^{2}}{\partial\theta^{2}}+\frac{N-1}{\theta}\frac{\partial}{\partial\theta}\right)$.
Also recall the “surface area” of the unit $N$-sphere $S_{N}$ is
$A_{N}=2\pi^{\left(N+1\right)/2}/\Gamma\left(\frac{N+1}{2}\right)$. This is
the “total solid angle” for Euclidean space $\mathbb{E}_{N+1}$.

A singular solution of the homogeneous equation

$0=\frac{1}{\sin^{N-1}\theta}\frac{\partial}{\partial\theta}\left(\left(\sin^{N-1}\theta\right)\frac{\partial}{\partial\theta}I\left(\theta,N\right)\right)$
(2)

is obviously given by the integral

$I\left(\theta,N\right)=\int_{\pi/2}^{\theta}\frac{1}{\left(\sin\phi\right)^{N-1}}d\phi$
(3)

_almost everywhere_ , except at the singularities $\theta=0$ or $\theta=\pi$.
All $I\left(N\right)$ for even $N$ involve logarithms, but all
$I\left(N\right)$ for odd $N$ do not.

Similarly, a singular solution of the inhomogeneous equation

$1=\frac{1}{\sin^{N-1}\theta}\frac{\partial}{\partial\theta}\left(\left(\sin^{N-1}\theta\right)\frac{\partial}{\partial\theta}J\left(\theta,N\right)\right)$
(4)

is just as obviously given almost everywhere by the double integral

$J\left(\theta,N\right)=\int_{\pi/2}^{\theta}\frac{1}{\left(\sin\vartheta\right)^{N-1}}\left(\int_{\pi/2}^{\vartheta}\left(\sin\phi\right)^{N-1}d\phi\right)d\vartheta$
(5)

By combining these two solutions we arrive at explicit expressions for a
Laplacian Green function on the $N$-sphere, in a form that exhibits the
physics of the model in simple terms. The result is

$\displaystyle
G\left(\widehat{r},\widehat{s},N\right)=\frac{1}{2A_{N-1}}~{}I\left(\arccos\left(\widehat{r}\cdot\widehat{s}\right),N\right)-\frac{1}{A_{N}}~{}J\left(\arccos\left(\widehat{r}\cdot\widehat{s}\right),N\right)$
(6) $\displaystyle
R^{2}\nabla^{2}G\left(\widehat{r},\widehat{s},N\right)=\delta^{N}\left(\widehat{r}-\widehat{s}\right)-\frac{1}{A_{N}}$
(7)

As usual in physics, $G\left(\widehat{r},\widehat{s},N\right)$ is interpreted
here as the potential at position $\widehat{r}$ produced by a unit point
charge at position $\widehat{s}$.

It is straightforward to express the results (6) in terms of hypergeometric
functions [4], which _Mathematica_ will do without much coaxing, but the
hypergeometric functions in question always reduce to combinations of
elementary functions. (See the Appendix for $N\leq 10$.) Alternatively, the
Green function can be written as a sum of bilinears in a complete set of
hyperspherical harmonics [1] divided by their Laplacian eigenvalues, but
excluding the zero mode solution [2], hence the $-1/A_{N}$ term in (7)

The Dirac delta produced by acting with $\nabla^{2}$ at the “north pole” of
the N-sphere is most easily exhibited by small angle expansions. For example:

$\displaystyle G\left(\widehat{r},\widehat{s}=\widehat{z},N=2\right)$
$\displaystyle=\frac{1}{4\pi}\ln\left(1-\cos\theta\right)=\frac{\ln\theta}{2\pi}-\frac{1}{4\pi}\ln
2-\frac{\theta^{2}}{48\pi}+O\left(\theta^{4}\right)$ (8) $\displaystyle
G\left(\widehat{r},\widehat{s}=\widehat{z},N=3\right)$
$\displaystyle=\frac{1}{4\pi^{2}}\left(\theta-\pi\right)\cot\theta=\allowbreak-\frac{1}{4\pi\theta}+\frac{1}{4\pi^{2}}+\frac{\theta}{12\pi}+O\left(\theta^{2}\right)$
(9)

etc. For $N>2$ the leading term is always
$-1/\left(\left(N-2\right)A_{N}~{}\theta^{N-2}\right)$, thereby revealing the
Dirac delta when the Laplacian acts as
$\frac{1}{R^{2}}\left(\frac{\partial^{2}}{\partial\theta^{2}}+\frac{N-1}{\theta}\frac{\partial}{\partial\theta}\right)$
for small $\theta$.

With $G\left(\widehat{r},\widehat{s},N\right)$ interpreted as a potential,
then the gradient of $G$ is to be interpreted as a repulsive/attractive force
exerted on a unit point charge/mass at position $\widehat{r}$ produced by an
identical point charge/mass at position $\widehat{s}$. Placing a source point
charge at the north pole of the N-sphere will produce a repulsive force on
another, same sign point charge at position $\left(\theta,\cdots\right)$. In
this particular situation, the repulsive force depends only on the polar angle
$\theta$, is independent of all the other “azimuthal” angles, and has only a
$\theta$ component as given by

$F\left(\theta,N\right)=\frac{d}{d\theta}G\left(\theta,N\right)=\frac{1}{4\pi^{N/2}}\left(\Gamma\left(\frac{N}{2}\right)-\frac{2}{\sqrt{\pi}}~{}\Gamma\left(\frac{N+1}{2}\right)K\left(\theta,N\right)\right)\frac{1}{\left(\sin\theta\right)^{N-1}}$
(10)

where
$K\left(\theta,N\right)=\int_{\pi/2}^{\theta}\left(\sin\phi\right)^{N-1}d\phi=\left(\sin\theta\right)^{N-1}\frac{d}{d\theta}J\left(\theta,N\right)$.
For example:

The result for $F\left(\theta,N\right)$ evinces Gauss’ law on the N-sphere, as
Faraday would surely have realized. The field line flux through any $S_{N-1}$
subsphere at polar angle $\theta$, of radius $R\sin\theta$, is $\Phi=\int
F\left(\theta,N\right)R^{N-1}\sin^{N-1}\theta
d^{N-1}\Omega=F\left(\theta,N\right)R^{N-1}A_{N-1}\sin^{N-1}\theta$. This flux
is proportional to the total charge on the hyperspherical cap
$\left.S_{N}\right|_{\leq\theta}$. Thus $\Phi$ decreases monotonically as
$\theta$ increases for $0<\theta<\pi$, and goes to zero as
$\theta\rightarrow\pi$ corresponding to zero total charge on the entire
hypersphere. For example:

$\displaystyle\int F\left(\theta,2\right)\sin\theta d\phi$
$\displaystyle=\frac{1}{2}\left(1+\cos\theta\right)=1-\int_{0}^{\theta}\frac{1}{4\pi}\sin\vartheta
d\vartheta\int_{0}^{2\pi}d\phi$ (11) $\displaystyle\int
F\left(\theta,3\right)\sin^{2}\theta d^{2}\Omega$
$\displaystyle=\frac{1}{\pi}\left(\pi-\theta+\cos\theta\sin\theta\right)=1-\int_{0}^{\theta}\frac{1}{2\pi^{2}}\sin^{2}\vartheta
d\vartheta\int d^{2}\Omega$ (12)

etc. The “$1$” on the RHS is the Dirac delta contribution, while the integral
is due to the uniform negative charge density on the hypersphere.

## The Shell Theorem

An interesting physical feature of this N-sphere model is embodied in the
“shell theorem” — or perhaps the “new shell theorem” would be a more
appropriate description.

In Euclidean space the well-known shell theorem is the statement that an ideal
uniform spherical shell of charge or mass will produce no force on a test
particle placed within the shell. Actually, the theorem is true for general
closed, charged, equipotential shells (somewhat obviously from the point of
view of the uniqueness theorem) as was first proven for spheroidal shells by
Isaac Newton, in the Principia, and much later for general ellipsoidal shells
by James Ivory, remarkably before potential theory was invented.

However, a uniformly charged $S_{N-1}$ sub-sphere located on the N-sphere
_will_ exert force on a test particle located at almost all points on
$S_{N}$.111More generally, on a generically curved, closed Riemannian
manifold, the interior of a charged, closed, equipotential submanifold is
_not_ necessarily at a _constant_ potential, despite what one might naively
expect from harmonic function theory. Recall the potential, identified as a
Green function, is _not_ necessarily harmonic in the interior of the
submanifold due to terms analogous to the $-1/A_{N}$ in (7). The only points
on $S_{N}$ where the force on the test particle will be zero are the “two
antipodal Faraday points” of the system, as well as all points on the
uniformly charged $S_{N-1}$ itself. Note that the antipodal Faraday points (or
“Far points” for brevity) may be defined as the two points on $S_{N}$ for
which all points on the embedded $S_{N-1}$ are equidistant, with each Far
point as far as possible from the sub-sphere as well as from each other. These
two Faraday points maximize, locally, their respective distances from the
$S_{N-1}$ sub-sphere. (Of course, to specify these Far points for the
N-sphere, distances may be computed either intrinsically or else
extrinsically, if in the latter case both the $S_{N-1}$ and the $S_{N}$
manifolds are canonically embedded in $\mathbb{E}_{N+1}$.)

For such a uniformly charged or massive $S_{N-1}\subset S_{N}$, the shell
theorem is supplanted by the following (new shell theorem):

###### Proposition 1

The force on any test particle placed on $S_{N}\smallsetminus S_{N-1}$ is the
same as would result if all the uniformly distributed charge or mass on
$S_{N-1}$ were moved away from the test particle and relocated entirely at the
opposing Faraday point.

To clarify the last _point_ , note that the sub-sphere partitions
$S_{N}\smallsetminus S_{N-1}$ into two disjoint regions, with one region
containing the test particle and the other containing the “opposing” Faraday
point. Thus the position of the test particle determines which of the Far
points is the opposing one.

So stated, the spherical shell theorem in Euclidean space becomes a special
case of the Proposition. If a test particle is placed outside a uniformly
charged $S_{N-1}$ spherical shell embedded in $\mathbb{E}_{N}$, the force on
that test particle is the same as though all the charge were concentrated at
the center of the shell, as is well-known. On the other hand, if the test
particle is placed inside the shell, there is no force, but that is exactly
the same null result that would be obtained if all the charge on the shell
were moved out to the “point at infinity”. In this case, the center of the
charged shell and the point at infinity play the roles of the antipodal
Faraday points for the system. That is to say, the Euclidean space shell
theorem for $\mathbb{E}_{N}$ results from taking the radius of $S_{N}$ to
infinity while keeping fixed the radius of the charged $S_{N-1}$ sub-sphere.

All this is most easily visualized and verified mathematically for the
2-sphere.

## The Simplicity of $N=2$

For a 2-sphere of radius $R$, let $\overrightarrow{r}$ and
$\overrightarrow{s}$ be two points on the sphere, so $r=R=s$. Then the
standard Green function for the Helmholtz equation is

$\displaystyle
G_{\text{Helmholtz}}\left(\overrightarrow{r},\overrightarrow{s}\right)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\frac{Y_{lm}\left(\widehat{r}\right)Y_{lm}^{\ast}\left(\widehat{s}\right)}{k^{2}-l\left(l+1\right)/R^{2}}=\sum_{l=0}^{\infty}\frac{2l+1}{4\pi}\frac{P_{l}\left(\widehat{r}\cdot\widehat{s}\right)}{k^{2}-l\left(l+1\right)/R^{2}}$
(13)
$\displaystyle\left(\nabla^{2}+k^{2}\right)G_{\text{Helmholtz}}\left(\overrightarrow{r},\overrightarrow{s}\right)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}Y_{lm}\left(\widehat{r}\right)Y_{lm}^{\ast}\left(\widehat{s}\right)=\delta^{2}\left(\widehat{r}-\widehat{s}\right)$
(14)

where $\widehat{r}$ has direction $\left(\theta,\phi\right)$, $\widehat{s}$
has direction $\left(\vartheta,\varphi\right)$, and
$\widehat{r}\cdot\widehat{s}=\cos\Theta$, with

$\cos\Theta=\cos\theta\cos\vartheta+\sin\theta\cos\phi\sin\vartheta\cos\varphi+\sin\theta\sin\phi\sin\vartheta\sin\varphi$
(15)

and
$\delta^{2}\left(\widehat{r}-\widehat{s}\right)=\frac{1}{\left|\sin\Theta\right|}~{}\delta\left(\Theta\right)\delta\left(\Phi\right)$.
Note that $\sqrt{1-\cos\Theta}$ is proportional to the length of the chord
(_not_ the arc length) that connects the two points on the sphere.
Unfortunately, _Mathematica_ cannot readily carry out the sum over $l$ to
obtain a simple closed form for $G_{\text{Helmholtz}}$, even after
regularizing the sum through use of an $i\epsilon$ prescription.

This is in contrast to the Laplacian Green function on the 2-sphere. In this
case _Mathematica_ can perform the relevant sum.

$\displaystyle
G_{\text{Laplace}}\left(\overrightarrow{r},\overrightarrow{s}\right)=\text{constant}+\sum_{\text{NB
}l=1}^{\infty}\sum_{m=-l}^{l}\frac{Y_{lm}\left(\widehat{r}\right)Y_{lm}^{\ast}\left(\widehat{s}\right)}{-l\left(l+1\right)/R^{2}}$
$\displaystyle=\text{constant}+\sum_{l=1}^{\infty}\frac{2l+1}{4\pi}\frac{P_{l}\left(\cos\Theta\right)}{-l\left(l+1\right)/R^{2}}$
$\displaystyle=\text{constant}+\frac{R^{2}}{4\pi}\left(\ln\left(1-\cos\Theta\right)+1-\ln
2\right)$ (16)

Note the issue with the $l=0$ mode (for example, see
https://en.wikipedia.org/wiki/Laplace-Beltrami_operator).

It is convenient to choose the constant to be $\frac{R^{2}}{4\pi}\left(-1+\ln
2\right)$ such that

$G_{\text{Laplace}}\left(\overrightarrow{r},\overrightarrow{s}\right)=\frac{R^{2}}{4\pi}\ln\left(1-\cos\Theta\right)$
(17)

Also note that interchange of the Laplacian with the sum gives

$\nabla^{2}G_{\text{Laplace}}\left(\overrightarrow{r},\overrightarrow{s}\right)=\sum_{\text{NB
}l=1}^{\infty}\sum_{m=-l}^{l}Y_{lm}\left(\widehat{r}\right)Y_{lm}^{\ast}\left(\widehat{s}\right)=\delta^{2}\left(\widehat{r}-\widehat{s}\right)-\frac{1}{4\pi}$
(18)

That is to say,
$\int_{S_{2}}\nabla^{2}G_{\text{Laplace}}\left(\overrightarrow{r},\overrightarrow{s}\right)d^{2}r=0$,
no matter what the location of $\overrightarrow{s}\in S_{2}$. With the
interpretation that $\nabla^{2}G_{\text{Laplace}}\propto\sigma$, the surface
charge (or mass) density, the fact that $\nabla^{2}G_{\text{Laplace}}$
integrates to zero means the total charge (or mass) on the sphere is always
zero (cf. Faraday’s field line interpretation). This _fact_ is quite curious
and may be surprising when first encountered. And while it is certainly not a
big deal for electric charge, it is rather more provocative for mass. It
remains to be seen if this fact requires modification of some long-established
ways of thinking about cosmology in a Newtonian framework [3].

In any case, we will _define_ the position dependent Newtonian gravitational
potential energy on the sphere, between point charges at
$\overrightarrow{r_{1}}$ and $\overrightarrow{r_{2}}$, with masses $m_{1}$ and
$m_{2}$, to be

$V\left(\overrightarrow{r_{1}},\overrightarrow{r_{2}}\right)=\kappa~{}m_{1}m_{2}\ln\left(1-\widehat{r_{1}}\cdot\widehat{r_{2}}\right)$
(19)

where $\kappa=G/R$ is the “Newtonian gravitational constant on the sphere”
with units for $G$ chosen to agree with those in 3D Euclidean space.

### 2-Sphere Shell theorem

Now it should be more or less obvious that the proof of the 2D Euclidean space
shell theorem, as given by Newton, does _not_ work on the 2-sphere.

For example, put a uniform ring of total mass $M$ but negligible thickness on
a circle of latitude specified by $\vartheta$. The Faraday points of this
system are obviously the north and south poles of the sphere. Next, consider
the force on an ideal point test particle, of mass $m$, located at
$\left(\theta,\phi\right)$. The uniform mass density of the ring is
$\lambda=M/\left(2\pi R\sin\vartheta\right)$ since its radius is
$R\sin\vartheta$, thus $dM=\frac{1}{2\pi}Md\phi$ for azimuthal angular
segments of the ring. Therefore the net force on the test particle is

$\overrightarrow{F}\left(\theta,\phi\right)=F_{\theta}\widehat{\theta}+F_{\phi}\widehat{\phi}=\frac{\kappa
mM}{2\pi
R}\int_{0}^{2\pi}\left(-\widehat{\theta}\frac{\partial}{\partial\theta}\ln\left(1-\widehat{r_{1}}\cdot\widehat{r_{2}}\right)-\widehat{\phi}\frac{1}{\sin\theta}\frac{\partial}{\partial\phi}\ln\left(1-\widehat{r_{1}}\cdot\widehat{r_{2}}\right)\right)d\varphi$
(20)

Clearly the $F_{\phi}$ component will integrate to zero, by symmetry, but
perhaps less obviously, for $0\leq\theta\leq\pi$,

$\displaystyle f\left(\theta,\vartheta\right)$ $\displaystyle\equiv 2\pi
RF_{\theta}/\left(\kappa mM\right)$ (21)
$\displaystyle=2\pi\left(\frac{1}{1+\cos\theta}\operatorname{Heaviside}\left(\vartheta-\theta\right)-\frac{1}{1-\cos\theta}\operatorname{Heaviside}\left(\theta-\vartheta\right)\right)\sin\theta$

In other words, there is in general an attractive force on the test particle,
due to this ring of total mass $M$, and that force is the same as the force
which would be produced by a _point_ charge of mass $M$, located at the
Faraday point in the spherical region on the other side of the ring, opposite
the test particle. That is to say,

$f\left(\theta,\vartheta\right)=f\left(\theta,\pi\right)\operatorname{Heaviside}\left(\vartheta-\theta\right)-f\left(\theta,0\right)\operatorname{Heaviside}\left(\theta-\vartheta\right)$
(22)

To help visualize this new shell theorem, we plot
$f\left(\theta,\vartheta\right)$ for various ring latitudes, along with the
force curves produced by point particles at either north (i.e.
$f\left(\theta,0\right)$) or south (i.e. $f\left(\theta,\pi\right)$) poles of
the sphere. First, an equatorial ring.

Next, a ring placed at the Tropic of Cancer, latitude $23^{\circ}\ 26\mrq\
10.8"\ =\ 23.43632^{\circ}$N.

Note the discontinuity in $f$ at the ring’s latitude is $4\pi/\sin\vartheta$.

Generalizing all this to higher N-spheres is now too obvious to warrant
further discussion here. But again, it remains to be seen if this requires
modification of some long-established ways of thinking about cosmology in a
Newtonian framework [3].

### Miscellaneous Additional Notes for the 2-Sphere

The perspicacious reader will have noticed that there is an obvious solution
to the inhomogeneous equation

$\frac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\left(\left(\sin\theta\right)\dfrac{\partial}{\partial\theta}h\left(\theta\right)\right)=-\dfrac{1}{4\pi}$
(23)

Namely,

$h\left(\theta\right)=\dfrac{1}{4\pi}\ln\left(\sin\theta\right)+const$ (24)

Combining with the previous result,

$\dfrac{1}{4\pi}\ln\left(1-\cos\theta\right)-\dfrac{1}{4\pi}\ln\left(\sin\theta\right)=\dfrac{1}{4\pi}\ln\left(\tan\frac{\theta}{2}\right)$
(25)

Check for $\theta\neq 0$:
$\dfrac{\partial}{\partial\theta}\left(\left(\sin\theta\right)\dfrac{\partial}{\partial\theta}\left(\dfrac{1}{4\pi}\left(\ln\left(\tan\left(\theta/2\right)\right)\right)\right)\right)=0$,
or at least this is true so long as $\theta\neq 0$. Also check the Dirac
delta:

$\int_{\text{cap at north
pole}}\nabla^{2}\ln\left(\tan\frac{\theta}{2}\right)d^{2}\Omega=\frac{1}{R^{2}}\int_{\partial\text{cap}}\left(\frac{\partial}{\partial\theta}\ln\left(\tan\frac{\theta}{2}\right)\right)\sin\theta
d\phi=\frac{2\pi}{R^{2}}$ (26)

upon using
$\left(\frac{\partial}{\partial\theta}\ln\left(\tan\frac{\theta}{2}\right)\right)\sin\theta=\allowbreak
1$. Therefore,
$\nabla^{2}\ln\left(\tan\frac{\theta}{2}\right)=\frac{2\pi}{R^{2}}~{}\delta^{2}\left(\widehat{z}\right)$,
or at least that is the case near the north pole of the sphere. What about the
south pole?

$\int_{\text{cap at south
pole}}\nabla^{2}\ln\left(\tan\frac{\theta}{2}\right)d^{2}\Omega=-\frac{1}{R^{2}}\int_{\partial\text{cap}}\left(\frac{\partial}{\partial\theta}\ln\left(\tan\frac{\theta}{2}\right)\right)\sin\theta
d\phi=-\frac{2\pi}{R^{2}}$ (27)

So there is another point source at the south pole. That is to say,

$\nabla^{2}\ln\left(\tan\frac{\theta}{2}\right)=\frac{2\pi}{R^{2}}~{}\left(\delta^{2}\left(\widehat{z}\right)-\delta^{2}\left(-\widehat{z}\right)\right)=\frac{4\pi}{R^{2}}~{}\delta^{2}\left(\widehat{z}\right)-\frac{2\pi}{R^{2}}~{}\left(\delta^{2}\left(\widehat{z}\right)+\delta^{2}\left(-\widehat{z}\right)\right)$
(28)

Subtracting the $\ln\left(\sin\theta\right)$ term cancelled the uniform
density on the sphere due to the $\ln\left(1-\cos\theta\right)$ term, but in
addition, the $\ln\left(\sin\theta\right)$ term contributed equal strength
Dirac deltas at both the north and the south poles, such that the overall
contribution of $\ln\left(\sin\theta\right)$ to the total charge was still
zero. Hence, the modification

$G\left(\overrightarrow{r},\overrightarrow{s}\right)=\frac{R^{2}}{4\pi}\ln\left(\tan\left(\frac{1}{2}\arccos\left(\widehat{r}\cdot\widehat{s}\right)\right)\right)$
(29)

yields a Green function whose source is only Dirac deltas, but again with
total charge zero.

$\nabla^{2}G\left(\overrightarrow{r},\overrightarrow{s}\right)=\frac{1}{2}~{}\delta^{2}\left(\widehat{r}-\widehat{s}\right)-\frac{1}{2}~{}\delta^{2}\left(\widehat{r}+\widehat{s}\right)\propto\sigma_{\text{total}}\text{
\ \ with \ \ }\int_{S_{2}}\sigma_{\text{total}}d^{2}r=0$ (30)

Yet another way to look at the $\ln\left(\sin\theta\right)$ term is to write

$\ln\left(\sin\theta\right)=\frac{1}{2}\ln\left(\sin^{2}\theta\right)=\frac{1}{2}\ln\left(1-\cos^{2}\theta\right)=\frac{1}{2}\left(\ln\left(1-\cos\theta\right)+\ln\left(1+\cos\theta\right)\right)$
(31)

which reveals the two Dirac deltas more immediately.

Alternatively, one could just as well _subtract_ from the
$\ln\left(1-\widehat{r}\cdot\widehat{s}\right)$ term a contribution due to a
single, equal strength Dirac delta at _any_ _other location_. For example, if
the original point source is at the north pole of the sphere, put an opposite
sign point source at the south pole, to obtain

$\ln\left(\frac{1-\cos\theta}{1+\cos\theta}\right)=\ln\left(\frac{\sin^{2}\theta/2}{\cos^{2}\theta/2}\right)=2\ln\left(\tan\frac{\theta}{2}\right)$
(32)

Thus, if $\theta\neq 0$ and $\theta\neq\pi$, once again check:
$\frac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\left(\left(\sin\theta\right)\dfrac{\partial}{\partial\theta}\ln\left(\tan\frac{\theta}{2}\right)\right)=\allowbreak
0$. Moreover,

$G\left(\overrightarrow{r},\overrightarrow{s}\right)=\frac{R^{2}}{2\pi}\ln\left(\tan\left(\frac{1}{2}\arccos\left(\widehat{r}\cdot\widehat{s}\right)\right)\right)\
,\ \ \
\nabla^{2}G\left(\overrightarrow{r},\overrightarrow{s}\right)=\delta^{2}\left(\widehat{r}-\widehat{s}\right)-\delta^{2}\left(\widehat{r}+\widehat{s}\right)$
(33)

In fact, an arbitrary number of Dirac deltas, representing point particles
with various strengths at various locations on $S_{2}$, will eliminate the
constant charge density provided all the Dirac delta charges sum to zero.
Similar remarks apply for other $S_{N}$.

## Conclusion

We identified Newtonian and Coulombic potentials, and forces, with Green
functions on $S_{N}$, and their gradients. We then considered the analogue of
Newton’s shell theorem for this simple system. We found a tasteful
generalization of the shell theorem which is simple to state. We believe a
similar shell theorem can be established for spheroidal and ellipsoidal if not
more general closed manifolds, as we intend to show in subsequent work.

Acknowledgement: This work resulted from participation in BASIC 2022. One of
us (TLC) thanks the participants of that conference for their comments on the
$N=2$ example.

## Appendix

The general integral expressions for the Green functions on $S_{N}$, as given
by (6) in the text, are not appropriate for the circle, $S_{1}$, when $\theta$
is allowed to be negative with $-\pi\leq\theta\leq\pi$. In that case $G$ is
best obtained from first principles, with the result

$G\left(\widehat{r},\widehat{s}=\widehat{z},N=1\right)\equiv
G\left(\theta\right)=\frac{1}{2}\left|\theta\right|-\frac{1}{4\pi}\theta^{2}\
,\ \ \ \text{for}\ \ \ -\pi\leq\theta\leq\pi$

The $\theta^{2}$ term eliminates a Dirac delta at $\theta=\pi$, implicit in
the second derivative of the $\left|\theta\right|$ term, by ensuring the slope
of $G\left(\theta\right)$ is continuous at the south pole of the circle, i.e.
$\lim_{\theta\rightarrow\pm\pi}dG\left(\theta\right)/d\theta=0$. Thus

$\frac{d^{2}}{d\theta^{2}}G\left(\theta\right)=\delta\left(\theta\right)-\frac{1}{2\pi}$

Otherwise, the general integral expressions in the text are easily evaluated
for specific integer $N\geq 2$. Here we list
$G\left(\widehat{r},\widehat{s}=\widehat{z},N\right)$ for $2\leq N\leq 10$, as
computed by Maple directly from those integral expressions.

$\displaystyle G\left(\widehat{r},\widehat{s}=\widehat{z},N=2\right)$
$\displaystyle=\frac{1}{4\pi}\ln\left(1-\cos\theta\right)$ $\displaystyle
G\left(\widehat{r},\widehat{s}=\widehat{z},N=3\right)$
$\displaystyle=\frac{\cos\theta}{4\pi^{2}\sin\theta}\left(\theta-\pi\right)$
$\displaystyle G\left(\widehat{r},\widehat{s}=\widehat{z},N=4\right)$
$\displaystyle=\frac{1}{8\pi^{2}}\left(\ln\left(1-\cos\theta\right)-\frac{\cos\theta}{1-\cos\theta}\right)$
$\displaystyle G\left(\widehat{r},\widehat{s}=\widehat{z},N=5\right)$
$\displaystyle=\frac{\cos\theta}{8\pi^{3}\sin^{3}\theta}\left(3\left(\theta-\pi\right)+2\left(\pi-\theta\right)\cos^{2}\theta-\sin\theta\cos\theta\right)$
$\displaystyle G\left(\widehat{r},\widehat{s}=\widehat{z},N=6\right)$
$\displaystyle=\frac{1}{16\pi^{3}}\left(3\ln\left(1-\cos\theta\right)+\frac{\left(4\cos\theta-5\right)\cos\theta}{\left(1-\cos\theta\right)^{2}}\right)$
$\displaystyle G\left(\widehat{r},\widehat{s}=\widehat{z},N=7\right)$
$\displaystyle=\frac{\cos\theta}{16\pi^{4}\sin^{5}\theta}\left(\left(\theta-\pi\right)\left(3+4\sin^{2}\theta+8\sin^{4}\theta\right)-3\left(\sin\theta\cos\theta\right)\left(1+2\sin^{2}\theta\right)\right)$
$\displaystyle G\left(\widehat{r},\widehat{s}=\widehat{z},N=8\right)$
$\displaystyle=\frac{1}{32\pi^{4}}\left(15\ln\left(1-\cos\theta\right)-\frac{\cos\theta}{\left(1-\cos\theta\right)^{3}}\left(33-54\cos\theta+23\cos^{2}\theta\right)\right)$
$G\left(\widehat{r},\widehat{s}=\widehat{z},N=9\right)=\\\
\frac{\cos\theta}{32\pi^{5}\sin^{7}\theta}\left(3\left(\theta-\pi\right)\left(5+6\sin^{2}\theta+8\sin^{4}\theta+16\sin^{6}\theta\right)-\left(\sin\theta\cos\theta\right)\left(15+28\sin^{2}\theta+44\sin^{4}\theta\right)\right)$
$G\left(\widehat{r},\widehat{s}=\widehat{z},N=10\right)=\\\
\frac{1}{64\pi^{5}}\left(105\ln\left(1-\cos\theta\right)\allowbreak+\frac{\cos\theta}{\left(1-\cos\theta\right)^{4}}\left(176\cos^{3}\theta-599\cos^{2}\theta+696\cos\theta-279\right)\right)$

## References

  * [1] H. Alshal and T. Curtright, “Grounded hyperspheres as squashed wormholes” J.Math.Phys. 60 (2019) 032901. https://arxiv.org/abs/1806.03762
  * [2] T. Aubin, Some Nonlinear Problems in Riemannian Geometry, Springer (1998), especially Chapter 4.
  * [3] C. Callan, R.H. Dicke, and P.J.E. Peebles, “Cosmology and Newtonian Mechanics” Am.J.Phys. 33 (1965) 105-108.
  * [4] R. Chapling, “A Hypergeometric Integral with Applications to the Fundamental Solution of Laplace’s Equation on Hyperspheres” SIGMA 12 (2016) 079.
  * [5] G. Green, “Mathematical investigations concerning the laws of the equilibrium of fluids analogous to the electric fluid, with other similar researches” Trans.Cambridge Phil.Soc. 5 (1835) 1–63.