The Port-Hamiltonian Structure of Continuum Mechanics

Abstract

In this paper, we present a novel approach to the geometric formulation of solid and fluid mechanics within the port-Hamiltonian framework, which extends the standard Hamiltonian formulation to non-conservative and open dynamical systems. Leveraging Dirac structures, instead of symplectic or Poisson structures, this formalism allows the incorporation of energy exchange within the spatial domain or through its boundary, which allows for a more comprehensive description of continuum mechanics. Building upon our recent work in describing nonlinear elasticity using exterior calculus and bundle-valued differential forms, this paper focuses on the systematic derivation of port-Hamiltonian models for solid and fluid mechanics in the material, spatial, and convective representations using Hamiltonian reduction theory. This paper also discusses constitutive relations for stress within this framework including hyper-elasticity, for both finite and infinitesimal strains, as well as viscous fluid flow governed by the Navier–Stokes equations.

1 Introduction

Due to the invariance of physical laws and Noether’s theorem, energy is at the heart of many physical theories and is fundamental to the Lagrangian and Hamiltonian approaches to continuum mechanics. In the Hamiltonian formulation, the governing equations of continuum mechanics are expressed using symplectic or Poisson structures which are geometric objects that embody the invariance of physical laws in the state space and describe the energetic structure of the equations of motion. This holds in all classical fields and clearly also for both fluid and solid mechanics in the spatial, convective, or material representation (Arnold 1965; Marsden and Weinstein 1974; Simo et al. 1988).

Due to the inherent skew-symmetric nature of symplectic and Poisson structures, the standard Hamiltonian formalism is limited to conservative closed systems that do not exchange energy with their surrounding. In fluid mechanics, this is usually imposed by assuming the velocity field to be tangent to the spatial domain’s boundary (Arnold 1965), whereas in solid mechanics one is limited to fixed boundary conditions (Simo et al. 1988).

On the contrary, the port-Hamiltonian formalism (Van Der Schaft and Maschke 2002) is applicable to non-conservative and open systems and relies on Dirac structures, instead of symplectic or Poisson structures, for characterizing the energetic structure of the equations of motion. This framework is suitable for describing fluid and solid mechanics with energy exchange within the spatial domain or through its boundary: it can incorporate state constraints (e.g., incompressibility) as well as viscous effects unlike the canonical Hamiltonian formalism. Being based on energy, which is the “lingua franca” of physics, the port-Hamiltonian formalism has been widely applied to numerous multi-physics domains such as fluid mechanics (Califano et al. 2021), solid mechanics (Brugnoli et al. 2021), fluid structure interaction (Califano et al. 2022), thermodynamics (Califano et al. 2022), and magneto-hydrodynamics (Vu et al. 2016).

This paper is a sequel to our recent work (Rashad et al. 2023) dealing with the geometric modeling of nonlinear elasticity using exterior calculus. In the current paper, we focus on the formulation of nonlinear elasticity in the port-Hamiltonian framework using bundle-valued differential forms. The main goal is to express the governing equations of elasticity as a network of energetic modules interconnected using power ports and Dirac structures, graphically represented in Fig. 1. Such identification of the energetic structure underlying the equations has numerous benefits. On the one hand, it provides conceptual insights to better understand the theory of continuum mechanics in a coordinate-free manner. On the other hand, it can be exploited for structure-preserving discretization, model order reduction, energy-based control, and analysis. The interested reader may refer to Rashad et al. (2020) for an extensive survey of the potential of the port-Hamiltonian framework.

The contributions of this paper are as follows:

• 1. We formulate the governing equations of motion of solid and fluid mechanics in the port-Hamiltonian framework.

• 2. In contrast with other efforts in the literature that usually rewrite the equations in a port-Hamiltonian form, we use Hamiltonian reduction theory to derive the governing equations from first principles. Demonstrating this port-Hamiltonian modeling process in itself is a contribution because it can be applied to other physical systems.

• 3. We treat the material, convective, and spatial representations of continuum mechanics and highlight the similarities between them explicated by our exterior calculus formulation.

• 4. We also provide a number of constitutive relations for stress showing how they fit within the port-Hamiltonian framework and highlighting the importance of the convective representation for constitutive modeling.

The outline of the paper is as follows: In Sect. 2, we provide a summary of the geometric formulation of continuum mechanics in tensor calculus and exterior calculus detailed in Rashad et al. (2023). In Sects. 3 and 4, we present the port-Hamiltonian model for nonlinear elasticity and fluid mechanics, separately to highlight the similarities and differences between them. We present the constitutive relations for hyper-elasticity and viscous fluid flow in Sect. 5 and conclude the paper in Sect. 6.

2 Geometric Formulation of Nonlinear Elasticity

In this section, we recall the geometric approach to continuum mechanics both in tensor calculus and exterior calculus. We refer the reader to Marsden and Hughes (1994) for an introduction to the tensor calculus formulation and to Rashad et al. (2023) for an intrinsic treatment of the topic in exterior calculus using bundle-valued forms.

2.1 Notation

Throughout this paper, we will distinguish material quantities by a $\sim$ on top, convective quantities by a $\wedge$ on top, and spatial quantities with none. For any manifold M, we shall denote the space of (real-valued) smooth functions on M by $C^\infty (M)$ and sections of any vector bundle $\mathbb {E}$ over M by $\varGamma (\mathbb {E})$. For instance, the space of p-contravariant and q-covariant smooth tensor fields is denoted by $\varGamma (T^p_q M)$. The space of (scalar-valued) differential k-forms, i.e., totally asymmetric k-covariant tensor fields, will be denoted by $\varOmega ^{k}(M)$. The space of generic two-point tensor fields on the manifolds M and N related by ${\varphi }:{M}\rightarrow {N}$ will be denoted by $\varGamma (T^p_q M \otimes \varphi ^* T^r_s N)$. If M is endowed with a Riemannian metric g, we shall interchangeably denote the index lowering operation of a vector field $u\in \varGamma (TM)$ to a one form in $\varOmega ^{1}(M)$ either by $u^\flat$ or $g\cdot u$. The latter explicit notation will be useful since we will use different metrics for the spatial, convective, and material representations.

2.2 Tensor Calculus

The configuration of an elastic body is described by a smooth embedding ${\varphi }:{\mathcal {B}}\rightarrow {\mathscr {A}}$ of the body manifold $\mathcal {B}$ into the ambient space $\mathscr {A}$. A motion of the elastic body is represented by a smooth curve ${c_\varphi }:{\mathbb {R}}\rightarrow {\mathscr {C}}$ with $\mathscr {C}$ denoting the configuration space of smooth embeddings of $\mathcal {B}$ in $\mathscr {A}$, i.e., $\mathscr {C}:= \text {Emb}^\infty (\mathcal {B},\mathscr {A})$. We denote the image of the whole body by $\mathcal {S}:= \varphi (\mathcal {B}) \subset \mathscr {A}$ and consider only the case $\dim (\mathcal {B}) = \dim (\mathcal {S}) = 3$. We assume $\mathscr {A}$ has a Riemannian structure with $g$ denoting its metric, while $\mathcal {B}$ does not have an intrinsic Riemannian structure and only inherits a configuration-dependent metric $\hat{g}:=\varphi _t^*(g)$ through the embedding $\varphi$ via pullback. Instead, $\mathcal {B}$ is equipped intrinsically with the mass form $\hat{\mu }\in \varOmega ^{n}(\mathcal {B})$.

The material (Lagrangian) velocity of the body is denoted by ${\tilde{v}_t}:{\mathcal {B}}\rightarrow {T\mathcal {S}}$ while the spatial (Eulerian) and convective velocities are denoted by $v_t\in \varGamma (T\mathcal {S})$ and ${\hat{v}}_t\in \varGamma (T\mathcal {B})$, respectively. The three representations of the velocity are related by

$$\begin{aligned} v_t:= \tilde{v}_t \circ \varphi {\scriptstyle ^{-1}_t}, \qquad \qquad {\hat{v}}_t:= T\varphi {\scriptstyle ^{-1}_t}\circ \tilde{v}_t = T\varphi {\scriptstyle ^{-1}_t}\circ v_t \circ \varphi _t, \end{aligned}$$

(1)

where ${T\varphi {\scriptstyle ^{-1}_t}}:{T\mathcal {S}}\rightarrow {T\mathcal {B}}$ denotes the tangent map of $\varphi {\scriptstyle ^{-1}_t}$. The convective, material, and spatial mass density functions will be denoted by $\hat{\rho }, \tilde{\rho }\in C^\infty (\mathcal {B}),\rho \in C^\infty (\mathcal {S})$ which are related by

$$\begin{aligned} \tilde{\rho }= J_{\varphi _t} \hat{\rho }_t = J_{\varphi _t}(\rho _t \circ \varphi _t), \end{aligned}$$

with $J_{\varphi _t}\in C^\infty (\mathcal {B})$ denoting the Jacobian of $\varphi _t$.

The material equations of motion governing the evolution of $(\varphi ,\tilde{v})\in T\mathscr {C}$ on the tangent bundle of $\mathscr {C}$ are given by

$$\begin{aligned} \partial _t \varphi =&\tilde{v}\end{aligned}$$

(2)

$$\begin{aligned} D_t \tilde{v}=&\frac{1}{\tilde{\rho }} \widetilde{\textrm{div}}(\tilde{\sigma })\end{aligned}$$

(3)

$$\begin{aligned} D_t F =&\tilde{\nabla }\tilde{v}\end{aligned}$$

(4)

$$\begin{aligned} \tilde{\sigma }=&\tilde{\rho }\ \left( \frac{\partial \tilde{e}}{\partial F}\right) ^\sharp (F), \end{aligned}$$

(5)

where $\tilde{\sigma }\in \varGamma (T\mathcal {B}\otimes \varphi ^*T\mathcal {S})$ denotes the 1^st Piola–Kirchhoff stress, $F:= T\varphi$ denotes the deformation gradient, $D_t$ denotes the material derivative, and $\tilde{e}(F) \in C^\infty (\mathcal {B})$ denotes the material internal energy function. Furthermore, $\tilde{\nabla }$ denotes the material covariant derivative, $\widetilde{\textrm{div}}:=\textrm{tr}\circ \tilde{\nabla }$ denotes its corresponding divergence operator, and the index raising (i.e., $\sharp$ map) in (5) is with respect to $g^{ij} \circ \varphi _t\in C^\infty (\mathcal {B})$.

The spatial equations of motion governing the evolution of $(\rho ,v) \in C^\infty (\mathcal {S})\times \varGamma (T\mathcal {S})$ are given by

$$\begin{aligned} \partial _t \rho =&-\textrm{div}(\rho v) \end{aligned}$$

(6)

$$\begin{aligned} \partial _t v=&- \nabla _{v}v+ \frac{1}{\rho } \textrm{div}(\sigma ) \end{aligned}$$

(7)

$$\begin{aligned} \sigma =&2 \rho \frac{\partial e}{\partial g}(F,g). \end{aligned}$$

(8)

where $\sigma \in \varGamma (T_0^2\mathcal {S})$ denotes the symmetric Cauchy stress tensor field and $e(F,g) \in C^\infty (\mathcal {S})$ denotes the spatial internal energy function. We denote by $\nabla$ the spatial covariant derivative and by $\textrm{div}:= \textrm{tr}\circ \nabla$ its corresponding divergence operator.

Finally, the convective counterparts of (6-8) that govern the evolution of $(\hat{\rho },{\hat{v}})\in C^\infty (\mathcal {B})\times \varGamma (T\mathcal {B})$ are given by

$$\begin{aligned} \partial _t \hat{\rho }=&- \hat{\rho }\widehat{\textrm{div}}({\hat{v}}) \end{aligned}$$

(9)

$$\begin{aligned} \partial _t {\hat{v}}=&- \hat{\nabla }_{{\hat{v}}}{\hat{v}}+ \frac{1}{\hat{\rho }} \widehat{\textrm{div}}(\hat{\sigma }) \end{aligned}$$

(10)

$$\begin{aligned} \partial _t \hat{g}=&\mathcal {L}_{{\hat{v}}}{\hat{g}} \end{aligned}$$

(11)

$$\begin{aligned} \hat{\sigma }=&2 \hat{\rho }\frac{\partial \hat{e}}{\partial \hat{g}} (\hat{g}) \end{aligned}$$

(12)

where $\hat{\sigma }\in \varGamma (T_0^2\mathcal {B})$ denotes the symmetric convective stress tensor field, $\hat{e}(\hat{g})\in C^\infty (\mathcal {B})$ denotes the convective internal energy function, and $\hat{g}:= \varphi _t^*(g) \in \mathcal {M}(\mathcal {B})$ denotes the convective metric. Furthermore, $\hat{\nabla }$ denotes the convective covariant derivative, and $\widehat{\textrm{div}}:=\textrm{tr}\circ \hat{\nabla }$ denotes its corresponding divergence operator. We denote by $\mathcal {M}(\mathcal {B})$ the space of Riemannian metrics on $\mathcal {B}$, which represents the space of deformations. As described in Rashad et al. (2023); Fiala (2016); Kolev and Desmorat (2021) and discussed later, strain is geometrically represented as a geodesic on this space. Thus, $\mathcal {M}(\mathcal {B})$ will play an important role for constitutive modeling in Sect. 5.

2.3 Exterior Calculus

The geometric formulation presented above relies on the representation of the physical variables appearing in the equations of motion as thermodynamically intensive (i.e., volume independent) variables. The mathematical objects used in such formulation consist of scalar fields, vector fields, and second-rank tensor fields. In the spatial, convective, and material representations, these, respectively, are $(\rho ,v,\sigma )$, $(\hat{\rho },{\hat{v}},\hat{\sigma })$, and $(\varphi ,\tilde{v},\tilde{\sigma })$, in addition to the variables that characterize elastic deformation.

A key limitation associated with the tensor calculus operations used in such geometric formulation is that the underlying geometric and topological structure of the equations of motion is entangled together. Exterior calculus on the other hand provides an elegant machinery for distinguishing between topology and geometry that can highlight the rich structures hidden in the partial differential equations above. The physical variables of continuum mechanics are represented in this exterior calculus formulation using differential forms, both scalar valued and bundle valued.

The space of scalar-valued differential k-forms on an n-dimensional manifold M is denoted by $\varOmega ^{k}(M)$. A k-form $\phi ^k\in \varOmega ^{k}(M)$ is locally a totally asymmetric (0, k) tensor with values in $\mathbb {R}$ at the point $p\in M$. One has that 0-forms are isomorphic to scalar fields and 1-forms isomorphic to covector fields, i.e., $\varOmega ^{0}(M)\cong C^\infty (M)$ and $\varOmega ^{1}(M)\cong \varGamma (T^*M)$. We denote the wedge product, exterior derivative, and Hodge star operators, respectively, by

$$\begin{aligned} & {\wedge }:{\varOmega ^{k}(M)\times \varOmega ^{l}(M)}\rightarrow {\varOmega ^{k+l}(M)} \\ & {\textrm{d}}:{\varOmega ^{k}(M)}\rightarrow {\varOmega ^{k+1}(M)} \qquad \qquad {\star }:{\varOmega ^{k}(M)}\rightarrow {\varOmega ^{n-k}(M)}. \end{aligned}$$

Other operations will be introduced later in the sequel.

The space of bundle-valued differential k-forms on M will be denoted by $\varOmega ^{k}(M;\mathbb {E})$ where $\mathbb {E}$ denotes any vector bundle over M. A k-form $\varPhi ^k\in \varOmega ^{k}(M;\mathbb {E})$ is locally a totally asymmetric (0, k) tensor that takes values in the fiber $\mathbb {E}$ at the point $p\in M$. One has that $\varOmega ^{0}(M;\mathbb {E}) \cong \varGamma (\mathbb {E})$; thus, $\varOmega ^{0}(M;TM)$ and $\varOmega ^{0}(M;T^*M)$ are isomorphic to vector fields and covector fields, respectively. Throughout this paper, we shall refer to elements of $\varOmega ^{0}(M;TM)$ and $\varOmega ^{0}(M;T^*M)$ as vector-valued and covector-valued forms, respectively. A trivial vector-(or covector-) valued k-form is one which is equivalent to the tensor product of a vector field (or covector field) and an ordinary k-form. For example, we say that $\zeta \in \varOmega ^{k}(M;TM)$ is trivial if it is composed of a vector field $u \in \varGamma (TM)$ and a k-form $\alpha ^k \in \varOmega ^{k}(M)$ such that $\zeta = \alpha ^k\otimes u.$

Bundle-valued forms can also elegantly incorporate two-point tensor fields that appear in the material representation of continuum mechanics. We shall denote the space of such forms by $\varOmega ^k_\varphi (M;\mathbb {F})$ where $\mathbb {F}$ is any vector bundle over the manifold N related to M via the map ${\varphi }:{M}\rightarrow {N}$. We denote the wedge dot product, exterior covariant derivative, and complementary Hodge star operators, respectively, by

$$\begin{aligned} & {\ \dot{\wedge }\ }:{\varOmega ^{k}(M;\mathbb {E})\times \varOmega ^{l}(M;\mathbb {E}^*)}\rightarrow {\varOmega ^{k+l}(M)}\\ & {\textrm{d}_{\nabla }}:{\varOmega ^{k}(M;\mathbb {E})}\rightarrow {\varOmega ^{k+1}(M;\mathbb {E})} \qquad \qquad {\star _c}:{\varOmega ^{k}(M;\mathbb {E})}\rightarrow {\varOmega ^{n-k}(M;\mathbb {E}^*)}. \end{aligned}$$

The exterior covariant derivative $\textrm{d}_{\nabla }$ generalizes the standard covariant derivative $\nabla$, and for the case $k=0$ one has that $\textrm{d}_{\nabla }= \nabla$. The exact definitions and coordinate-based expressions of the above operators¹ can be found in Rashad et al. (2023).

In our exterior calculus formulation of continuum mechanics we shall represent physical variables as thermodynamically extensive variables using a combination of scalar-valued and bundle-valued forms. In particular, we shall represent the spatial mass, momentum, and stress by

$$\begin{aligned} (\mu ,\mathcal {M},\mathcal {T}) \in \varOmega ^{n}(\mathcal {S})\times \varOmega ^{n}(\mathcal {S};T^*\mathcal {S})\times \varOmega ^{n-1}(\mathcal {S};T^*\mathcal {S}). \end{aligned}$$

Furthermore, we shall represent the convective metric, momentum, and stress variables by

$$\begin{aligned} (\hat{g},\hat{\mathcal {M}},\hat{\mathcal {T}}) \in \mathcal {M}(\mathcal {B})\times \varOmega ^{n}(\mathcal {B};T^*\mathcal {B})\times \varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B}), \end{aligned}$$

with the tangent and cotangent bundles of $\mathcal {M}(\mathcal {B})$ identified with $\varOmega _{\textrm{sym}}^{1}(\mathcal {B};T^*\mathcal {B})\subset \varOmega ^{1}(\mathcal {B};T^*\mathcal {B})$ and $\varOmega _{\textrm{sym}}^{n-1}(\mathcal {B};T\mathcal {B})\subset \varOmega ^{n-1}(\mathcal {B};T\mathcal {B})$, respectively. Working with symmetric subsets of the above spaces follows from the geometric structure of $\mathcal {M}(\mathcal {B})$ that comprises symmetric positive definite tensor fields (Rashad et al. 2023). Finally, we shall represent the configuration, material momentum, and stress variables by

$$\begin{aligned} (\varphi ,\tilde{\mathcal {M}},\tilde{\mathcal {T}}) \in \mathscr {C}\times \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})\times \varOmega _\varphi ^{n-1}(\mathcal {B};T^*\mathcal {S}), \end{aligned}$$

with the tangent and cotangent bundles of $\mathscr {C}$ identified with $\varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})$ and $\varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})$, respectively.

2.4 Related Works

Our geometric formulation based on exterior calculus is a sequel to Rashad et al. (2023) which is founded on the work of Frankel (2019); Kanso et al. (2007); Gilbert and Vanneste (2023) where stress and momentum are treated as bundle-valued differential forms, which is the natural mathematical representation of these physical variables from a topological perspective.

The majority of geometric formulations of continuum mechanics are based on tensor calculus. In the Hamiltonian mechanics literature, solid and fluid mechanics have been studied extensively using Poisson reduction in Simo et al. (1988); Marsden and Weinstein (1982, 1983); Marsden et al. (1984); Holm et al. (1986); Lewis et al. (1986); Mazer and Ratiu (1989). On the Lagrangian side, Euler–Poincare reduction has been also applied, e.g., in Holm et al. (1998); Gay-Balmaz et al. (2012). In contrast with these works, our approach to deriving the governing equations of motion in the port-Hamiltonian framework is based on Dirac structures and uses the bottom-up philosophy of port-based modeling to represent the dynamics decomposed into a number of energetic subsystems connected to each other using so-called power ports, or ports for short.

In previous works of ours, we developed port-Hamiltonian models for Euler equations (i.e., ideal fluid flow) using scalar-valued forms only in Rashad et al. (2021a, 2021b) which was later extended to Navier–Stokes equations in Califano et al. (2021); Rashad et al. (2021), to Fourier-Navier–Stokes in Califano et al. (2022), and fluid–structure interaction in Califano et al. (2022). In these extensions, the momentum balance was usually written in terms of scalar-valued forms with only the stress treated as a bundle-valued form. On the other hand, the port-Hamiltonian model we derive in this paper will represent both momentum and momentum flux as bundle-valued forms which will highlight more their geometric characteristics. This article builds on the philosophy of applying the port-Hamiltonian framework for developing decomposed dynamical models, as established in our aforementioned previous works on fluid mechanics. However, the current port-Hamiltonian formulation extends far beyond those earlier contributions, offering a more general approach that is broadly applicable to both solid and fluid mechanics. Furthermore, it is worth mentioning that in Rashad et al. (2023) we did not discuss at all the underlying port-Hamiltonian structure of nonlinear elasticity.

Other works in the port-Hamiltonian literature that treat fluid mechanics and simplified elastic models for beams and plates can be found in Cheng et al. (2023); Trivedi et al. (2016); Brugnoli et al. (2019, 2021); Macchelli and Melchiorri (2004) In contrast with these works that usually mathematically manipulate the governing equations of motion into a port-Hamiltonian form, we shall use in this work Hamiltonian reduction theory to derive the port-Hamiltonian models of continuum mechanics from first principles in a geometric coordinate-free manner.

Fig. 1

Convective port-Hamiltonian model of nonlinear elasticity showing the kinetic energy, stress power, and constitutive relations subsystems

3 Port-Hamiltonian Model of Nonlinear Elasticity

In this section, we present the port-Hamiltonian dynamic model of nonlinear elasticity in the material, spatial, and convective representations. The unique characteristic of the port-Hamiltonian paradigm that distinguishes it from the standard Hamiltonian approach in Simo et al. (1988); Lewis et al. (1986); Mazer and Ratiu (1989) is that the physical system is modeled as a network of energetic units interconnected by ports.

Each port $(f,e) \in V\times V^*$ consists of a pair of power conjugate variables $e\in V^*$ and $f\in V$, referred to, respectively, as effort and flow variables in port-based modeling. In our work, the space of flow variables V will be given by vector-valued forms while its dual space $V^*$ will be given by covector-valued pseudo-forms of complementary degree. The duality pairing of the port variables e and f is denoted by $\langle e | f \rangle _{M}:= \int _M f\ \dot{\wedge }\ e$ and characterizes the power flowing in the port (f, e). One has that $M=\mathcal {B}$ in case (f, e) are represented in the material or convective representation, while $M=\mathcal {S}$ in the case of the spatial representation. On some occasions, we will also represent the port variables using scalar valued forms. With a slight abuse of notation, we will also denote their corresponding duality pairing by $\langle e | f \rangle _{M}:= \int _M f\wedge e$, as it will be clear from the context.

In the standard Hamiltonian formalism, one follows a top-down approach by starting from a Hamiltonian functional that characterizes the elastic body’s kinetic energy, strain energy, and boundary conditions. The governing equations are represented using a Poisson structure such that the Hamiltonian is conserved (Simo et al. 1988). On the other hand, in the port-Hamiltonian formalism, we follow a bottom-up approach by representing each energetic subsystem separate from the others using a Dirac structure that characterizes the subsystem’s unique power balance. For nonlinear elasticity, its port-Hamiltonian model will consist of 1) a kinetic energy subsystem, 2) a stress power subsystem, and 3) a stress constitutive relation subsystem, as depicted in Fig. 1. By interconnecting these three subsystems to each other using power ports, the result is a decomposed model of the overall system that explicates the underlying energetic structure of the theory and emphasizes the natural duality of the port variables.

We shall present the three representations of the aforementioned energetic units using the exterior calculus formulation introduced earlier. We also show how one recovers the mass, momentum and energy balance laws from the presented port-Hamiltonian models. In fact, the port-Hamiltonian procedure we will present will allow us to derive the equations of motion from first principles using Hamiltonian reduction techniques. The reader is assumed to be familiar with the standard Hamiltonian framework of mechanics, e.g., in Marsden and Ratiu (2013); Holm et al. (2009).

3.1 Frechet and Variational Derivatives

Before proceeding, we introduce the Frechet derivative and variational derivatives of functionals of differential forms which will be used extensively in this paper. Consider any functional ${\mathscr {F}}:{\varOmega ^{k}(M;\mathbb {E})}\rightarrow {\mathbb {R}}$ with $\mathbb {E}$ denoting a vector bundle over the manifold M. The Frechet derivative $\textrm{D}_\alpha \mathscr {F}$ of $\mathscr {F}$ with respect to $\alpha \in \varOmega ^{k}(M;\mathbb {E})$ is defined for any $\delta \alpha \in \varOmega ^{k}(M;\mathbb {E})$ as the functional

$$\begin{aligned} \textrm{D}_\alpha \mathscr {F}[\alpha ,\delta \alpha ]:= \left. \frac{d}{ds}\right| _{s=0}\mathscr {F}[\alpha + s \delta \alpha ], \end{aligned}$$

(13)

while the variational derivative ${\delta _{\alpha }{\mathscr {F}}}:{\varOmega ^{k}(M;\mathbb {E})}\rightarrow {\varOmega ^{n-k}(M;\mathbb {E}^*)}$ is defined by

$$\begin{aligned} \int _M \delta \alpha \ \dot{\wedge }\ \delta _{\alpha }{\mathscr {F}}(\alpha ) = \textrm{D}_\alpha \mathscr {F}[\alpha ,\delta \alpha ], \end{aligned}$$

(14)

for any $\delta \alpha \in \varOmega ^{k}(M;\mathbb {E})$.

The above definitions can be trivially extended to define 1) partial Frechet derivatives, 2) partial variational derivatives, and 3) functionals of scalar-valued forms. Furthermore, we shall deal later with the case that $\alpha$ belongs to an infinite-dimensional manifold or a bundle, which will be essential for the port-Hamiltonian modeling procedure. For notational simplicity in these more complex cases, we shall opt for denoting $\textrm{D}_\alpha \mathscr {F}[\alpha ,\delta \alpha ]$ as $\textrm{D}_\alpha \mathscr {F}(\alpha ) \cdot \delta \alpha$ following (Marsden and Hughes 1994) and usually denote $\delta _{\alpha }{\mathscr {F}}(\alpha )$ simply by $\delta _{\alpha }{\mathscr {F}}$.

3.2 Kinetic Energy Subsystem

For an elastic body undergoing a motion described by ${\varphi _t}:{\mathcal {B}}\rightarrow {\mathcal {S}}$, its kinetic energy as a function of time is expressed using spatial, material, and convective variables, respectively, as

$$\begin{aligned} {E}_{\textrm{kin}}(t) = \int _\mathcal {S} \frac{1}{2}g(v_t,v_t) \mu _t = \int _\mathcal {B} \frac{1}{2}(g\circ \varphi _t)(\tilde{v}_t,\tilde{v}_t) \tilde{\mu }= \int _\mathcal {B} \frac{1}{2}\hat{g}_t({\hat{v}}_t,{\hat{v}}_t) \hat{\mu }, \end{aligned}$$

(15)

where $\tilde{\mu }= \hat{\mu }\in \varOmega ^{n}(\mathcal {B})$ denote the (time-independent) intrinsic mass form associated o the body manifold. The second equality in (15) follows from the change of variables formula, and the third one can be seen in local components:

$$(g\circ \varphi _t)(\tilde{v},\tilde{v}) = (g_{ij}\circ \varphi _t) \tilde{v}^i\tilde{v}^j = (g_{ij}\circ \varphi _t)F^i_I F^j_J {\hat{v}}^I{\hat{v}}^J = \hat{g}_{IJ}{\hat{v}}^I{\hat{v}}^J = \hat{g}({\hat{v}},{\hat{v}}).$$

In the spatial representation, the kinetic energy has state dependency on the spatial mass form and velocity field $(\mu _t,v_t)\in \varOmega ^{n}(\mathcal {S})\times \varGamma (T\mathcal {S})$ and parametric dependency on the Riemannian metric $g\in \mathcal {M}(\mathcal {S})$. In the material representation, it has state dependency on the configuration and material velocity field $(\varphi _t,\tilde{v}_t) \in T\mathscr {C}$ and parametric dependency on $(g,\tilde{\mu }) \in \mathcal {M}(\mathcal {S})\times \varOmega ^{n}(\mathcal {B})$. In the convective representation, it has state dependency on the convective metric and velocity field $(\hat{g}_t,{\hat{v}}_t)\in \mathcal {M}(\mathcal {B})\times \varGamma (T\mathcal {B})$ and parametric dependency on $\hat{\mu }\in \varOmega ^{n}(\mathcal {B})$.

The kinetic energy gives rise to three Lagrangian functionals on the state spaces detailed above. In what follows, we shall derive their Hamiltonian counterparts and the associated energetic structures that characterize conservation of energy. We start with the material representation exploiting the canonical Poisson structure on the cotangent bundle $T^*\mathscr {C}$. Then we present the Hamiltonian reduction procedure to the spatial representation, followed by the convective one. We perform this reduction by direct calculation similar to Simo et al. (1988); Mazer and Ratiu (1989) with two distinctions. First, instead of using tensor calculus, our formulation is based fully on exterior calculus which is (arguably) more mathematically elegant because it explicates the intrinsic relation between differential forms and integration. Second, instead of applying Hamiltonian reduction to the total energy of the system, we follow the bottom-up approach described earlier which simplifies the reduction process greatly as we deal with the kinetic energy separately.

3.2.1 Material Poisson Structure

Let ${\tilde{\mathcal {P}}}_{\textrm{k}}:=\mathcal {M}(\mathcal {S})\times \varOmega ^{n}(\mathcal {B})$ denote the space of parameters in the material representation and its state space given by the tangent bundle $T\mathscr {C}$. The material Lagrangian functional ${{\tilde{L}}_{\textrm{k}}}:{T\mathscr {C}\times {\tilde{\mathcal {P}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ is defined by

$$\begin{aligned} {\tilde{L}}_{\textrm{k}}[\varphi ,\tilde{v};g,\tilde{\mu }]:= \int _{\mathcal {B}}\frac{1}{2}(g\circ \varphi )(\tilde{v},\tilde{v}) \tilde{\mu }, \end{aligned}$$

(16)

which depends parametrically on $(g,\tilde{\mu }) \in {\tilde{\mathcal {P}}}_{\textrm{k}}$ while its time dependency is implicit through its state variables $(\varphi ,\tilde{v}) \in T\mathscr {C}$. In what follows, we will omit the time dependency and usually suppress the parametric dependency unless needed. It is important to note that $T\mathscr {C}$ is not a product space and thus caution should be made when defining variations (Simo et al. 1988).

In terms of the exterior calculus construction in Sect. 2.2, we can rewrite (16) as

$$\begin{aligned} {\tilde{L}}_{\textrm{k}}[\varphi ,\tilde{v}]:= \int _{\mathcal {B}}\frac{1}{2}\tilde{v}\ \dot{\wedge }\ \tilde{\star }_c\tilde{v}, \end{aligned}$$

(17)

where we consider $\tilde{v}\in \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})$ while the state dependency on $\varphi$ and parametric dependency on $g$ and $\tilde{\mu }$ have been absorbed in the material Hodge star operator ${\tilde{\star }_c}:{\varOmega _\varphi ^{k}(\mathcal {B};T\mathcal {S})}\rightarrow {\varOmega _\varphi ^{n-k}(\mathcal {B};T^*\mathcal {S})}$.

By applying a partial Legendre transformation on $\tilde{v}$, we define on the cotangent bundle $T^*\mathscr {C} =: {\tilde{\mathcal {X}}}_{\textrm{k}}$ the material Hamiltonian functional ${{\tilde{H}}_{\textrm{k}}}:{{\tilde{\mathcal {X}}}_{\textrm{k}}\times {\tilde{\mathcal {P}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ as

$$\begin{aligned} {\tilde{H}}_{\textrm{k}}[\varphi ,\tilde{\mathcal {M}}]:= \int _{\mathcal {B}}\frac{1}{2}\tilde{\star }_c^{-1}\tilde{\mathcal {M}}\ \dot{\wedge }\ \tilde{\mathcal {M}}, \end{aligned}$$

(18)

where $\tilde{\mathcal {M}}:= \tilde{\star }_c\tilde{v}\in \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})$ denotes the material momentum of the body and ${\tilde{\star }_c^{-1}}:{\varOmega _\varphi ^{k}(\mathcal {B};T^*\mathcal {S})}\rightarrow {\varOmega _\varphi ^{n-k}(\mathcal {B};T\mathcal {S})}$ is the inverse map of $\tilde{\star }_c$.

Remark 1

As shown in Rashad et al. (2023), $\tilde{\mathcal {M}}$ is locally expressed as $\tilde{\mathcal {M}}= \tilde{v}_i \tilde{\mu }\otimes e^i$ where $\tilde{v}_i:= (g_{ij}\circ \varphi ) \tilde{v}^j \in \varOmega ^{0}(\mathcal {B})$ denote the components of the covector velocity field $\tilde{v}^\flat$. Note that $\tilde{\mathcal {M}}$ is not a trivial covector-valued top form, i.e., $\tilde{\mathcal {M}}\ne \tilde{\mu }\otimes \tilde{v}^\flat$. Indeed, both expressions are equivalent when treated as a two-point tensor, cf. (Kanso et al. 2007, eq. 16). However, as a bundle-valued form its form part is given by $\tilde{v}_i \tilde{\mu }\in \varOmega ^{n}(\mathcal {B})$.

In the theory of Hamiltonian mechanics, the canonical equations of motion on ${\tilde{\mathcal {X}}}_{\textrm{k}}$ are given implicitly by $\dot{\tilde{\mathscr {F}}} = ^m\{\tilde{\mathscr {F}},{\tilde{H}}_{\textrm{k}}\}$, where ${\tilde{\mathscr {F}}}:{{\tilde{\mathcal {X}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ is an arbitrary functional and ${ ^m\{\cdot ,\cdot \}}:{C^\infty ({\tilde{\mathcal {X}}}_{\textrm{k}})\times C^\infty ({\tilde{\mathcal {X}}}_{\textrm{k}})}\rightarrow {\mathbb {R}}$ denotes the canonical Poisson bracket in the material representation, to be introduced shortly. To write the explicit expressions of the above canonical equations and bracket, it is essential to introduce the partial Frechet and variational derivatives of functionals on the cotangent bundle ${\tilde{\mathcal {X}}}_{\textrm{k}}$.

In our work, we identify tangent and cotangent spaces of $\mathscr {C}$ by $T_{\varphi }\mathscr {C}\cong \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})$ and $T_{\varphi }^*\mathscr {C}\cong \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})$ (Rashad et al. 2023). Consequently, we have that

$$\begin{aligned} T_{\tilde{\chi }}{\tilde{\mathcal {X}}}_{\textrm{k}}\cong&\ \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S}) \times \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S}) \times \varOmega _\varphi ^{0}(\partial \mathcal {B};T\mathcal {S}) \\ T^*_{\tilde{\chi }}{\tilde{\mathcal {X}}}_{\textrm{k}}\cong&\ \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S}) \times \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S}) \times \varOmega _\varphi ^{n-1}(\partial \mathcal {B};T^*\mathcal {S}). \end{aligned}$$

A functional ${\tilde{\mathscr {F}}}:{{\tilde{\mathcal {X}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ is said to have partial variational derivatives, if for any $(\varphi ,\tilde{\mathcal {M}}) \in {\tilde{\mathcal {X}}}_{\textrm{k}}$ there exist the bundle-valued forms

$$\delta _{\varphi }{\tilde{\mathscr {F}}} \in \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S}), \qquad \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}} \in \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S}), \qquad \delta ^\cup _{\varphi }{\tilde{\mathscr {F}}} \in \varOmega _\varphi ^{n-1}(\partial \mathcal {B};T^*\mathcal {S}), $$

that satisfy for any $(\delta \tilde{\varphi },\delta \tilde{\mathcal {M}}) \in \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})\times \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})$

$$\begin{aligned} \textrm{D}_\varphi \tilde{\mathscr {F}}(\varphi ,\tilde{\mathcal {M}}) \cdot \delta \tilde{\varphi } =&\int _{\mathcal {B}}\delta \tilde{\varphi } \ \dot{\wedge }\ \delta _{\varphi }{\tilde{\mathscr {F}}} + \int _{\partial \mathcal {B}} \delta \tilde{\varphi }|_{{\partial \mathcal {B}}} \ \dot{\wedge }\ \delta ^\cup _{\varphi }{\tilde{\mathscr {F}}} \end{aligned}$$

(19)

$$\begin{aligned} \textrm{D}_{\tilde{\mathcal {M}}}\tilde{\mathscr {F}}(\varphi ,\tilde{\mathcal {M}}) \cdot \delta \tilde{\mathcal {M}}=&\int _{\mathcal {B}}\delta \tilde{\mathcal {M}}\ \dot{\wedge }\ \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}}, \end{aligned}$$

(20)

where for any $\tilde{\alpha } \in \varOmega _\varphi ^{k}(\mathcal {B};T\mathcal {S})$, we denote by $\tilde{\alpha }|_{{\partial \mathcal {B}}}:= i_\text {f}^*(\tilde{\alpha }) \in \varOmega _\varphi ^{k}(\partial \mathcal {B};T\mathcal {S})$ the trace of the form part of $\tilde{\alpha }$ to the boundary where ${i}:{{\partial \mathcal {B}}}\rightarrow {\mathcal {B}}$ denotes the inclusion map and $i_\text {f}^*$ denotes the pullback of the form part of a bundle-valued form.

The partial Frechet derivative $\textrm{D}_{\tilde{\mathcal {M}}}\tilde{\mathscr {F}}$ is defined such that

$$\begin{aligned} \textrm{D}_{\tilde{\mathcal {M}}}\tilde{\mathscr {F}}(\varphi ,\tilde{\mathcal {M}}) \cdot \delta \tilde{\mathcal {M}}:= \left. \frac{d}{ds}\right| _{s=0}\tilde{\mathscr {F}}[\varphi ,\tilde{\mathcal {M}}+ s \delta \tilde{\mathcal {M}}]. \end{aligned}$$

(21)

Intuitively, $\textrm{D}_{\tilde{\mathcal {M}}}\tilde{\mathscr {F}}$ accounts to varying $\tilde{\mathcal {M}}$ as a covector in $T^*_\varphi \mathscr {C}$ while keeping $\varphi$ fixed, i.e., a fiber derivative. On the other hand, to define the partial Frechet derivative $\textrm{D}_\varphi \tilde{\mathscr {F}}$, one must intuitively fix the covector $\tilde{\mathcal {M}}$ while allowing the base point $\varphi$ to vary. This must be done with caution since ${\tilde{\mathcal {X}}}_{\textrm{k}}= T^*\mathscr {C}$ is not a product space similar to $T\mathscr {C}$. Let the variation $\delta \tilde{\varphi }\in T_\varphi \mathscr {C}\cong \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})$ be the tangent vector to the smooth curve $s \mapsto \varphi _s \in \mathscr {C}$ at $s=0$, i.e., $\delta \tilde{\varphi }:= \left. \frac{d}{ds}\right| _{s=0}\varphi _s$ with $\varphi _s|_{s=0} = \varphi$. The partial Frechet derivative $\textrm{D}_{\varphi }\tilde{\mathscr {F}}$ is defined such that

$$\begin{aligned} \textrm{D}_{\varphi }\tilde{\mathscr {F}}(\varphi ,\tilde{\mathcal {M}}) \cdot \delta \tilde{\varphi }:= \left. \frac{d}{ds}\right| _{s=0}\tilde{\mathscr {F}}[\varphi _s,\tilde{\mathcal {M}}_s], \end{aligned}$$

(22)

where $\tilde{\mathcal {M}}_s$ would be the induced variation on the covector in $T_{\varphi }^*\mathscr {C}$ due to varying $\varphi$. The details of the above construction can be found in Lewis et al. (1986); Mazer and Ratiu (1989) which we shall refer to later. Consequently, the rate of change of any ${\tilde{\mathscr {F}}}:{{\tilde{\mathcal {X}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ with respect to time is expressed as

$$\begin{aligned} \dot{\tilde{\mathscr {F}}} = \int _{\mathcal {B}}\partial _t \varphi \ \dot{\wedge }\ \delta _{\varphi }{\tilde{\mathscr {F}}} + \textrm{D}_t\tilde{\mathcal {M}}\ \dot{\wedge }\ \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}} + \int _{\partial \mathcal {B}} \partial _t\varphi |_{{\partial \mathcal {B}}} \ \dot{\wedge }\ \delta ^\cup _{\varphi }{\tilde{\mathscr {F}}}, \end{aligned}$$

(23)

with $\partial _t \varphi \in \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})$ and $\textrm{D}_t \tilde{\mathcal {M}}\in \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})$.

The canonical Poisson bracket on the cotangent bundle ${\tilde{\mathcal {X}}}_{\textrm{k}}$ then takes the form Lewis et al. (1986)

$$\begin{aligned} \begin{aligned} ^m\{\tilde{\mathscr {F}},\tilde{\mathscr {G}}\}:=&\int _{\mathcal {B}}\delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {G}}} \ \dot{\wedge }\ \delta _{\varphi }{\tilde{\mathscr {F}}} - \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}} \ \dot{\wedge }\ \delta _{\varphi }{\tilde{\mathscr {G}}} \\&+ \int _{\partial \mathcal {B}} \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {G}}}|_{\partial \mathcal {B}} \ \dot{\wedge }\ \delta ^\cup _{\varphi }{\tilde{\mathscr {F}}} - \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}}|_{\partial \mathcal {B}} \ \dot{\wedge }\ \delta ^\cup _{\varphi }{\tilde{\mathscr {G}}}. \end{aligned} \end{aligned}$$

(24)

Using the above construction, the canonical equations of motion are given explicitly by the following result.

Proposition 1

For the kinetic energy as a Hamiltonian functional with state variables $\tilde{\chi }:= (\varphi ,\tilde{\mathcal {M}}) \in {\tilde{\mathcal {X}}}_{\textrm{k}}$, the equations of motion in the material representation are given by:

$$\begin{aligned} \begin{pmatrix} \partial _t \varphi \\ \textrm{D}_t \tilde{\mathcal {M}} \end{pmatrix} = { \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}} \begin{pmatrix} \delta _{\varphi }{{\tilde{H}}_{\textrm{k}}}\\ \delta _{\tilde{\mathcal {M}}}{{\tilde{H}}_{\textrm{k}}} \end{pmatrix}, \qquad \partial _t \varphi |_{{\partial \mathcal {B}}} = \delta _{\tilde{\mathcal {M}}}{{\tilde{H}}_{\textrm{k}}}|_{{\partial \mathcal {B}}}, \qquad 0 = \delta ^\cup _{\varphi }{{\tilde{H}}_{\textrm{k}}}. \end{aligned}$$

(25)

The Hamiltonian functional (18) admits its rate of change such that along trajectories $(\varphi (t),\tilde{\mathcal {M}}(t))$ of (25), it holds that

$$\begin{aligned} {\dot{\tilde{H}}}_{\textrm{k}}= 0. \end{aligned}$$

(26)

Proof

First, by comparing (23) and (24), one can easily deduce the equations of motion (25). As for the energy balance (26), it follows immediately from the skew symmetry of the bracket (24) that ${\dot{\tilde{H}}}_{\textrm{k}}= ^m\{{\tilde{H}}_{\textrm{k}},{\tilde{H}}_{\textrm{k}}\} = 0$. $\square$

Corollary 1

The variational derivatives of the functional ${\tilde{H}}_{\textrm{k}}$ in (18) with respect to $\varphi \in \mathscr {C}$ and $\tilde{\mathcal {M}}\in \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})$ are given, respectively, by

$$\begin{aligned}\delta _{\varphi }{{\tilde{H}}_{\textrm{k}}} = 0 \in \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S}), \qquad \delta _{\tilde{\mathcal {M}}}{{\tilde{H}}_{\textrm{k}}} = \tilde{v}\in \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S}), \qquad \delta ^\cup _{\varphi }{{\tilde{H}}_{\textrm{k}}} = 0 \in \varOmega _\varphi ^{n-1}(\partial \mathcal {B};T^*\mathcal {S}).\end{aligned}$$

Consequently, one can rewrite the Hamilton equations (25) as

$$\partial _t \varphi = \tilde{v}, \quad D_t \tilde{\mathcal {M}}= 0, \quad \partial _t \varphi |_{{\partial \mathcal {B}}} = \tilde{v}|_{{\partial \mathcal {B}}},$$

where the first and third equations represent the kinematic definition of the material velocity while the second represents the conservation of momentum.

Proof

Second, by construction from the Legendre transformation, we have that $\delta _{\tilde{\mathcal {M}}}{{\tilde{H}}_{\textrm{k}}} = \tilde{v}\in \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S}).$

In Simo et al. (1988), the procedure for computing $\delta _{\varphi }{{\tilde{H}}_{\textrm{k}}}$ and $\delta ^\cup _{\varphi }{{\tilde{H}}_{\textrm{k}}}$ has been detailed. The construction is involved due to addressing the whole dynamic system at once. However, in our case where the Hamiltonian functional consists only of kinetic energy, one can see from (Simo et al. 1988, Prop. 3.1) that both $\delta _{\varphi }{{\tilde{H}}_{\textrm{k}}}$ and $\delta ^\cup _{\varphi }{{\tilde{H}}_{\textrm{k}}}$ are equal to 0. $\square$

In summary, the closed Hamiltonian system, defined by the tuple $({\tilde{\mathcal {X}}}_{\textrm{k}},{\tilde{H}}_{\textrm{k}}, ^m\{\cdot ,\cdot \})$, describes the conservation of energy and the corresponding evolution of the state $\tilde{\chi }\in {\tilde{\mathcal {X}}}_{\textrm{k}}$ as a conservative system isolated from any energy exchange with its surroundings.

Remark 2

Note that the Hamilton equations of motion (25) are valid for any continuum (whether solid or fluid) and for any Hamiltonian ${\tilde{H}}:{T^*\mathscr {C}}\rightarrow {\mathbb {R}}$. However, the condition $\delta ^\cup _{\varphi }{\tilde{H}} = 0$ immediately restricts the class of admissible functions such that the Hamiltonian is conserved. It will be shown later in the spatial representation that this condition amounts to having vanishing momentum flux at the image of the boundary in the ambient space. While this condition is naturally satisfied for solid mechanics in Corollary 1, it will have interesting implications for the case of fluid mechanics that we discuss in Sect. 4.

3.2.2 Material Dirac Structure

Dirac structures are geometric objects that generalize symplectic and Poisson structures and play a central role in the theory of port-Hamiltonian systems. Let $T{\tilde{\mathcal {X}}}_{\textrm{k}}\oplus T^*{\tilde{\mathcal {X}}}_{\textrm{k}}$ denote the Whitney sum bundle over ${\tilde{\mathcal {X}}}_{\textrm{k}}= T^*\mathscr {C}$, i.e., the bundle over the base ${\tilde{\mathcal {X}}}_{\textrm{k}}$ with fiber over $\tilde{\chi }:= (\varphi ,\tilde{\mathcal {M}})$ equal to $T_{\tilde{\chi }}{\tilde{\mathcal {X}}}_{\textrm{k}}\times T_{\tilde{\chi }}^*{\tilde{\mathcal {X}}}_{\textrm{k}}$. As per the construction in (19-23), the duality pairing between any $\tilde{f} :=(\tilde{f}_{\varphi },\tilde{f}_{\tilde{\mathcal {M}}},\tilde{f}_{\varphi }^\cup ) \in T_{\tilde{\chi }}{\tilde{\mathcal {X}}}_{\textrm{k}}$ and $\tilde{e} :=(\tilde{e}_{\varphi },\tilde{e}_{\tilde{\mathcal {M}}},\tilde{e}_{\varphi }^\cup ) \in T^*_{\tilde{\chi }}{\tilde{\mathcal {X}}}_{\textrm{k}}$ is expressed as

$$\langle \tilde{e} | \tilde{f} \rangle _{T{\tilde{\mathcal {X}}}_{\textrm{k}}}:= \langle \tilde{e}_{\varphi } | \tilde{f}_{\varphi } \rangle _{\mathcal {B}} + \langle \tilde{e}_{\tilde{\mathcal {M}}} | \tilde{f}_{\tilde{\mathcal {M}}} \rangle _{\mathcal {B}} + \langle \tilde{e}_{\varphi }^\cup | \tilde{f}_{\varphi }^\cup \rangle _{\partial \mathcal {B}}.$$

A Dirac structure on ${\tilde{\mathcal {X}}}_{\textrm{k}}$ is defined as the sub-bundle $\tilde{\mathcal {D}} \subset T{\tilde{\mathcal {X}}}_{\textrm{k}}\oplus T^*{\tilde{\mathcal {X}}}_{\textrm{k}}$ such that $\tilde{\mathcal {D}} = \tilde{\mathcal {D}}^\perp$, with $\tilde{\mathcal {D}}^\perp$ denoting the annihilator with respect to the symmetric bilinear form:

$$\langle \langle (\tilde{f} ^1,\tilde{e} ^1) | (\tilde{f} ^2,\tilde{e} ^2) \rangle \rangle := \langle {\tilde{e} ^1}|{\tilde{f} ^2}\rangle _{T{\tilde{\mathcal {X}}}_{\textrm{k}}} + \langle {\tilde{e} ^2}|{\tilde{f} ^1}\rangle _{T{\tilde{\mathcal {X}}}_{\textrm{k}}},$$

for all $(\tilde{f} ^1,\tilde{e} ^1), (\tilde{f} ^2,\tilde{e} ^2) \in T_{\tilde{\chi }}{\tilde{\mathcal {X}}}_{\textrm{k}}\times T_{\tilde{\chi }}^*{\tilde{\mathcal {X}}}_{\textrm{k}}.$

The canonical Dirac structure on the cotangent bundle ${\tilde{\mathcal {X}}}_{\textrm{k}}$ corresponding to the Poisson bracket (24) will be denoted by $\tilde{J}_\text {k} \subset T{\tilde{\mathcal {X}}}_{\textrm{k}}\oplus T^*{\tilde{\mathcal {X}}}_{\textrm{k}}$ whose fiber at any $\tilde{\chi }\in {\tilde{\mathcal {X}}}_{\textrm{k}}$ is given by:

$$\begin{aligned} \begin{aligned} \tilde{J}_\text {k}(\tilde{\chi }):=&\left\{ (\tilde{f}_{\text {k}},\tilde{e}_{\text {k}}) \in T_{\tilde{\chi }}{\tilde{\mathcal {X}}}_{\textrm{k}}\times T_{\tilde{\chi }}^*{\tilde{\mathcal {X}}}_{\textrm{k}}| \right. \\&\begin{pmatrix} \tilde{f}_{\varphi }\\ \tilde{f}_{\tilde{\mathcal {M}}} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \tilde{e}_{\varphi }\\ \tilde{e}_{\tilde{\mathcal {M}}} \end{pmatrix},\\&\left. \tilde{f}_{\varphi }^\cup = \tilde{e}_{\tilde{\mathcal {M}}}|_{{\partial \mathcal {B}}}, \qquad \qquad 0 = \tilde{e}_{\varphi }^\cup \right\} , \end{aligned} \end{aligned}$$

(27)

which characterizes the power balance $\langle \tilde{e}_{\text {k}} | \tilde{f}_{\text {k}} \rangle _{T{\tilde{\mathcal {X}}}_{\textrm{k}}} = 0$. The Hamilton equations in Prop. 1 immediately follow by setting $\tilde{e}_{\text {k}} = (\delta _{\varphi }{{\tilde{H}}_{\textrm{k}}},\delta _{\tilde{\mathcal {M}}}{{\tilde{H}}_{\textrm{k}}},\delta ^\cup _{\varphi }{{\tilde{H}}_{\textrm{k}}})$ as inputs and $\tilde{f}_{\text {k}} = (\partial _t\varphi ,\partial _t \tilde{\mathcal {M}}, \partial _t\varphi |_{{\partial \mathcal {B}}})$ as outputs yielding the conservation of the Hamiltonian (18).

Fig. 2

Material port-Hamiltonian model of the kinetic energy subsystem expressed in terms of block diagrams (left) and bond graphs (right)

Unlike Poisson structures that characterize conservative dynamical systems, Dirac structures can be used to characterize open dynamical systems that allow for nonzero energy exchange either through the boundary or within the spatial domain. Let $\tilde{\mathcal {F}}\in \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})$ be an external body force field and ${\tilde{\mathcal {I}}}_{\textrm{k}}:= \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S}) \times \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})$ denote the interaction space such that $(\tilde{v},\tilde{\mathcal {F}}) \in {\tilde{\mathcal {I}}}_{\textrm{k}}$ define an open external port used to model the interaction with other subsystems within the spatial domain.² This interaction is characterized by the power given by the duality pairing $\langle \tilde{\mathcal {F}} | \tilde{v} \rangle _{\mathcal {B}}$ and will be used later to model internal stress forces. In principle, one could add any body force, e.g., electrostatics or gravity in the same manner.

Now we can extend the canonical Dirac structure (27) with an interaction port to construct the Dirac structure ${\tilde{\mathcal {D}}}_{\textrm{k}}\subset {\tilde{\mathfrak {F}}}_{\textrm{k}}\oplus {\tilde{\mathfrak {F}}^*}_{\textrm{k}}$, with ${\tilde{\mathfrak {F}}}_{\textrm{k}}:= T{\tilde{\mathcal {X}}}_{\textrm{k}}\times \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})$ and its dual space ${\tilde{\mathfrak {F}}^*}_{\textrm{k}}:= T^*{\tilde{\mathcal {X}}}_{\textrm{k}}\times \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})$, whose fiber at any $\tilde{\chi }\in {\tilde{\mathcal {X}}}_{\textrm{k}}$ is defined as

$$\begin{aligned} \begin{aligned} {\tilde{\mathcal {D}}}_{\textrm{k}}(\tilde{\chi }) =&\left\{ ((\tilde{f}_{\text {k}}, \tilde{f}_{\text {d}}),(\tilde{e}_{\text {k}}, \tilde{e}_{\text {d}})) \in {\tilde{\mathfrak {F}}}_{\textrm{k}}(\tilde{\chi })\times {\tilde{\mathfrak {F}}^*}_{\textrm{k}}(\tilde{\chi }) | \right. \\&\begin{pmatrix} \tilde{f}_{\varphi }\\ \tilde{f}_{\tilde{\mathcal {M}}} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \tilde{e}_{\varphi }\\ \tilde{e}_{\tilde{\mathcal {M}}} \end{pmatrix} + \begin{pmatrix} 0\\ 1 \end{pmatrix} \tilde{e}_{\text {d}},\\&\qquad \tilde{f}_{\text {d}}= \begin{pmatrix} 0&1 \end{pmatrix} \begin{pmatrix} \tilde{e}_{\varphi }\\ \tilde{e}_{\tilde{\mathcal {M}}} \end{pmatrix}, \\&\quad \left. \tilde{f}_{\varphi }^\cup = \tilde{e}_{\tilde{\mathcal {M}}}|_{{\partial \mathcal {B}}}, \qquad \qquad 0 = \tilde{e}_{\varphi }^\cup \right\} , \end{aligned} \end{aligned}$$

(28)

with $\tilde{f}_{\text {k}} = (\tilde{f}_{\varphi },\tilde{f}_{\tilde{\mathcal {M}}},\tilde{f}_{\varphi }^\cup ) \in T_{\tilde{\chi }}{\tilde{\mathcal {X}}}_{\textrm{k}}$ and $\tilde{e}_{\text {k}} = (\tilde{e}_{\varphi },\tilde{e}_{\tilde{\mathcal {M}}},\tilde{e}_{\varphi }^\cup )\in T^*_{\tilde{\chi }}{\tilde{\mathcal {X}}}_{\textrm{k}}$. The Dirac structure (28) characterizes the power balance

$$\begin{aligned} \langle \tilde{e}_{\varphi } | \tilde{f}_{\varphi } \rangle _{\mathcal {B}} + \langle \tilde{e}_{\tilde{\mathcal {M}}} | \tilde{f}_{\tilde{\mathcal {M}}} \rangle _{\mathcal {B}} + \langle \tilde{e}_{\varphi }^\cup | \tilde{f}_{\varphi }^\cup \rangle _{{\partial \mathcal {B}}} = \langle \tilde{e}_{\text {d}} | \tilde{f}_{\text {d}} \rangle _{\mathcal {B}}. \end{aligned}$$

(29)

By setting the inputs to be $\tilde{e}_{\text {k}} = ({\delta _{\varphi }{{\tilde{H}}_{\textrm{k}}}},{\delta _{\tilde{\mathcal {M}}}{{\tilde{H}}_{\textrm{k}}}},{\delta ^\cup _{\varphi }{{\tilde{H}}_{\textrm{k}}}})$ and $\tilde{e}_{\text {d}} = \tilde{\mathcal {F}},$ and the outputs to be $\tilde{f}_{\text {k}}= ({\partial _t \varphi },{D_t \tilde{\mathcal {M}}},\partial _t \varphi |_{{\partial \mathcal {B}}})$ and $\tilde{f}_{\text {d}}= \tilde{v}$, one extends (25) to become a port-Hamiltonian system such that the energy balance (26) becomes

$$\begin{aligned} {\dot{\tilde{H}}}_{\textrm{k}}= \langle \tilde{\mathcal {F}} | \tilde{v} \rangle _{\mathcal {B}}, \end{aligned}$$

(30)

which states that the rate of change of kinetic energy of the elastic body is equal to the work done due to external stress forces. Consequently, the balance of momentum in Corollary 1 takes the form $D_t \tilde{\mathcal {M}}= \tilde{\mathcal {F}}$.

To summarize, with reference to Fig. 2, the open port-Hamiltonian system is given by the tuple $({\tilde{\mathcal {X}}}_{\textrm{k}},{\tilde{H}}_{\textrm{k}},{\tilde{\mathcal {D}}}_{\textrm{k}},{\tilde{\mathcal {I}}}_{\textrm{k}})$ where $({\tilde{\mathcal {X}}}_{\textrm{k}},{\tilde{H}}_{\textrm{k}})$ represent the storage of kinetic energy, ${\tilde{\mathcal {I}}}_{\textrm{k}}$ represents the power supplied due to stress, and ${\tilde{\mathcal {D}}}_{\textrm{k}}$ represents the balance laws of the system.

3.2.3 Spatial Dirac Structure

Now we turn attention to the spatial counterpart of the kinetic energy port-Hamiltonian system introduced above. We start by transforming the material variables to spatial ones by the nonlinear diffeomorphism

$$\begin{aligned} \begin{aligned} ^s\varPhi : {{\tilde{\mathcal {X}}}_{\textrm{k}}\times {\tilde{\mathcal {P}}}_{\textrm{k}}}&\rightarrow {{{\mathcal {X}}}_{\textrm{k}}\times {{\mathcal {P}}}_{\textrm{k}}}\\ {((\varphi ,\tilde{\mathcal {M}}),(g,\tilde{\mu }))}&\mapsto {((\mu ,\mathcal {M}),g)} \end{aligned} \end{aligned}$$

(31)

with $\mu := \varphi _*( \tilde{\mu })\in \varOmega ^{n}(\mathcal {S})$ denoting the extensive mass form and $\mathcal {M}:= \varphi _{\textrm{f},*}(\tilde{\mathcal {M}}) \in \varOmega ^{n}(\mathcal {S};T^*\mathcal {S})$ denoting the extensive momentum form in the spatial representation. We denote the spatial state space by ${{\mathcal {X}}}_{\textrm{k}}:= \varOmega ^{n}(\mathcal {S})\times \varOmega ^{n}(\mathcal {S};T^*\mathcal {S})$ and the spatial parameter space by ${{\mathcal {P}}}_{\textrm{k}}:= \mathcal {M}(\mathcal {S})$.

The spatial kinetic energy Hamiltonian ${{{H}}_{\textrm{k}}}:{{{\mathcal {X}}}_{\textrm{k}}\times {{\mathcal {P}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ is defined by ${{H}}_{\textrm{k}}[\mu ,\mathcal {M};g] = {\tilde{H}}_{\textrm{k}}[\varphi ,\tilde{\mathcal {M}};g,\tilde{\mu }]$ and is expressed as

$$\begin{aligned} {{H}}_{\textrm{k}}[\mu ,\mathcal {M}]:= \int _{\mathcal {S}}\frac{1}{2}\star _c^{-1}\mathcal {M}\ \dot{\wedge }\ \mathcal {M}, \end{aligned}$$

(32)

where the spatial Hodge star ${\star _c}:{\varOmega ^{k}(\mathcal {S};T\mathcal {S})}\rightarrow {\varOmega ^{n-k}(\mathcal {S};T^*\mathcal {S})}$ incorporates a state dependency on the mass form $\mu$ and a parametric dependency on the spatial metric $g$. Note that in contrast with the material and convective cases which were represented by bundle-valued forms, the spatial case is represented by a combination of scalar-valued and bundle-valued forms.

Remark 3

In local coordinates, we have that the spatial momentum is computed from its material counterpart (cf. . Remark 1) by

$$\mathcal {M}= \varphi _{\textrm{f},*}(\tilde{\mathcal {M}})= \varphi _*(\tilde{v}_i \tilde{\mu }) \otimes e^i = \varphi _*(\tilde{v}_i) \varphi _*(\tilde{\mu }) \otimes e^i = v_i\mu \otimes e^i,$$

and thus, $\mathcal {M}$ is identified with the trivial covector-valued form $\mu \otimes v^\flat$ which is in contrast with $\tilde{\mathcal {M}}$.

One can easily show the equivalence between ${{H}}_{\textrm{k}}$ and the standard kinetic energy expression in (15) by

$${{H}}_{\textrm{k}}= \int _{\mathcal {S}}\frac{1}{2}v\ \dot{\wedge }\ \mathcal {M}= \int _{\mathcal {S}}\frac{1}{2}v\ \dot{\wedge }\ (\mu \otimes v^\flat ) =\int _{\mathcal {S}}\frac{1}{2}v^\flat (v)\mu = \int _{\mathcal {S}}\frac{1}{2}g(v,v)\mu .$$

In general, the spatial counterpart ${\mathscr {F}}:{{{\mathcal {X}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ of any material functional ${\tilde{\mathscr {F}}}:{{\tilde{\mathcal {X}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ can be defined such that $\mathscr {F}:= \tilde{\mathscr {F}} \circ ^s\varPhi ^{-1}$, i.e.,

$$\begin{aligned} \mathscr {F}[\varphi _*( \tilde{\mu }),\varphi _{\textrm{f},*}(\tilde{\mathcal {M}})] = \tilde{\mathscr {F}}[\varphi ,\tilde{\mathcal {M}}]. \end{aligned}$$

(33)

The variational derivatives of the functional $\mathscr {F}$ with respect to $(\mu ,\mathcal {M})$ are the scalar-valued and bundle-valued forms $\delta _{\mu }{\mathscr {F}} \in \varOmega ^{0}(\mathcal {S})$ and $\delta _{\mathcal {M}}{\mathscr {F}} \in \varOmega ^{0}(\mathcal {S};T\mathcal {S}),$ that satisfy for any $\delta \mu \in \varOmega ^{n}(\mathcal {S}),\delta \mathcal {M}\in \varOmega ^{n}(\mathcal {S};T^*\mathcal {S})$

$$\begin{aligned} \textrm{D}_\mu \mathscr {F}(\mu ,\mathcal {M}) \cdot \delta \mu = \int _{\mathcal {S}}\delta \mu \wedge \delta _{\mu }{\mathscr {F}}, \qquad \textrm{D}_{\mathcal {M}}\mathscr {F}(\mu ,\mathcal {M}) \cdot \delta \mathcal {M}= \int _{\mathcal {S}}\delta \mathcal {M}\ \dot{\wedge }\ \delta _{\mathcal {M}}{\mathscr {F}}. \nonumber \\ \end{aligned}$$

(34)

The rate of change of $\mathscr {F}$ with respect to time is then expressed as

$$\begin{aligned} \dot{\mathscr {F}} = \langle \delta _{\mu }{\mathscr {F}} | \partial _t \mu \rangle _{\mathcal {S}} + \langle \delta _{\mathcal {M}}{\mathscr {F}} | \partial _t \mathcal {M} \rangle _{\mathcal {S}}, \end{aligned}$$

(35)

with $\partial _t \mu \in \varOmega ^{n}(\mathcal {S})$ and $\partial _t \mathcal {M}\in \varOmega ^{n}(\mathcal {S};T^*\mathcal {S})$.

The pushforward and pullback maps of (31) are given by the following important result which is essential for the derivation of the spatial Dirac structure. We recall first a number of operators for bundle-valued forms from Rashad et al. (2023). For any vector bundle $\mathbb {E}_\mathcal {S}$ over $\mathcal {S}$, we denote by ${\varphi _{\textrm{f}}^*}:{\varOmega ^{k}(\mathcal {S};\mathbb {E}_\mathcal {S})}\rightarrow {\varOmega _\varphi ^{k}(\mathcal {B};\mathbb {E}_\mathcal {S})}$ the pullback of the form part of a bundle-valued form, transforming it from the spatial representation to the material one, and denote by $\varphi _{\textrm{f},*}$ its inverse (i.e., the pushforward map). For any $u \in \varGamma (T\mathcal {S})$, we denote by ${\iota _u}:{\varOmega ^{k}(\mathcal {S};\mathbb {E}_\mathcal {S})}\rightarrow {\varOmega ^{k-1}(\mathcal {S};\mathbb {E}_\mathcal {S})}$ the interior product of a covector-valued form defined by inserting u into its form part and we denote by $\mathcal {L}_{u}{(\cdot )}$ the Lie derivative of any tensor field along u.

Proposition 2

Let $(\varphi ,\tilde{\mathcal {M}}) \in {\tilde{\mathcal {X}}}_{\textrm{k}}$ and $(\mu ,\mathcal {M})\in {{\mathcal {X}}}_{\textrm{k}}$ be related by the diffeomorphism $^s\varPhi$ in (31).

i) The pushforward map of $^s\varPhi$ is given by

$$\begin{aligned} \begin{aligned} ^s\varPhi _* : {T{\tilde{\mathcal {X}}}_{\textrm{k}}}&\rightarrow {T{{\mathcal {X}}}_{\textrm{k}}}\\ {(\varphi ,\tilde{\mathcal {M}},\delta \tilde{\varphi },\delta \tilde{\mathcal {M}},\delta \tilde{\varphi }|_{{\partial \mathcal {B}}})}&\mapsto {(\mu ,\mathcal {M},\delta \mu ,\delta \mathcal {M})} \end{aligned} \end{aligned}$$

with

$$\delta \mu = -\mathcal {L}_{\delta \varphi }{\mu } \in \varOmega ^{n}(\mathcal {S}),\qquad \qquad \delta \mathcal {M}= -\textrm{d}_{\nabla }(\iota _{\delta \varphi }\mathcal {M}) + \varphi _{\textrm{f},*}(\delta \tilde{\mathcal {M}}) \in \varOmega ^{n}(\mathcal {S};T^*\mathcal {S}),$$

where $\delta \varphi = \varphi _{\textrm{f},*}(\delta \tilde{\varphi })\in \varOmega ^{0}(\mathcal {S};T\mathcal {S}) \cong \varGamma (T\mathcal {S})$.

ii) The pullback map of $^s\varPhi$ is given by

$$\begin{aligned} \begin{aligned} ^s\varPhi ^* : {T^*{{\mathcal {X}}}_{\textrm{k}}}&\rightarrow {T^*{\tilde{\mathcal {X}}}_{\textrm{k}}}\\ {(\mu ,\mathcal {M},\delta _{\mu }{\mathscr {F}},\delta _{\mathcal {M}}{\mathscr {F}})}&\mapsto {(\varphi ,\tilde{\mathcal {M}},\delta _{\varphi }{\tilde{\mathscr {F}}}, \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}}, \delta ^\cup _{\varphi }{\tilde{\mathscr {F}}})} \end{aligned} \end{aligned}$$

with

$$\begin{aligned} \delta _{\varphi }{\tilde{\mathscr {F}}} =&\varphi _{\textrm{f}}^*(\mu \otimes (\textrm{d}\delta _{\mu }{\mathscr {F}} + \nabla \delta _{\mathcal {M}}{\mathscr {F}} \ \dot{\wedge }\ v^\flat )) & \in \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S}) \\ \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}} =&\varphi _{\textrm{f}}^*(\delta _{\mathcal {M}}{\mathscr {F}}) & \in \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})\\ \delta ^\cup _{\varphi }{\tilde{\mathscr {F}}} =&-\varphi _{\textrm{f}}^*(\delta _{\mu }{\mathscr {F}} \mu + \iota _{\delta _{\mathcal {M}}{\mathscr {F}}}\mathcal {M})|_{{\partial \mathcal {B}}} & \in \varOmega _\varphi ^{n-1}(\partial \mathcal {B};T^*\mathcal {S}). \end{aligned}$$

Proof

See Appendix A.1. $\square$

From the Hamiltonian reduction theory, the diffeomorphism (31) induces a Poisson bracket on ${{\mathcal {X}}}_{\textrm{k}}$ using the canonical Poisson bracket (24) on ${\tilde{\mathcal {X}}}_{\textrm{k}}$ such that

$$\begin{aligned} ^s\{\mathscr {F},\mathscr {G}\} \circ ^s \varPhi := ^m\{\mathscr {F}\circ ^s \varPhi ,\mathscr {G}\circ ^s \varPhi \}. \end{aligned}$$

(36)

Similarly, $^s\varPhi$ induces from (28) a spatial Dirac structure that extends the Poisson structure (36) with an interaction port. This spatial Dirac structure will be denoted by ${{\mathcal {D}}}_{\textrm{k}}$ and is derived in the following result.

Theorem 1

The spatial Dirac structure corresponding to (28) under the diffeomorphism (31) is the sub-bundle ${{\mathcal {D}}}_{\textrm{k}}\subset {{\mathfrak {F}}}_{\textrm{k}}\oplus {{\mathfrak {F}}^*}_{\textrm{k}}$ with

$$\begin{aligned} {{\mathfrak {F}}}_{\textrm{k}}:=& ^s \varPhi _*(T{\tilde{\mathcal {X}}}_{\textrm{k}}) \times \varphi _{\textrm{f},*}(\varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})) = T{{\mathcal {X}}}_{\textrm{k}}\times \varOmega ^{0}(\mathcal {S};T\mathcal {S})\\ {{\mathfrak {F}}^*}_{\textrm{k}}:=&( ^s \varPhi ^*)^{-1}(T^*{\tilde{\mathcal {X}}}_{\textrm{k}}) \times \varphi _{\textrm{f},*}(\varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})) = T^*{{\mathcal {X}}}_{\textrm{k}}\times \varOmega ^{n}(\mathcal {S};T^*\mathcal {S}), \end{aligned}$$

whose fiber at any ${\chi }:= (\mu ,\mathcal {M})\in {{\mathcal {X}}}_{\textrm{k}}$ is defined by

$$\begin{aligned} \begin{aligned} {{\mathcal {D}}}_{\textrm{k}}({\chi }):=&\left\{ (({f}_{\textrm{k}},{f}_{\textrm{d}}),({e}_{\textrm{k}},{e}_{\textrm{d}})) \in {{\mathfrak {F}}}_{\textrm{k}}({\chi })\times {{\mathfrak {F}}^*}_{\textrm{k}}({\chi }) | \right. \\&\begin{pmatrix} {f}_{\mu }\\ {f}_{\mathcal {M}} \end{pmatrix} = \begin{pmatrix} 0 & - \textrm{d}\iota _{(\cdot )}\mu \\ -\mu \otimes \textrm{d}(\cdot ) & -\mathcal {L}_{(\cdot )}{\mathcal {M}} \end{pmatrix} \begin{pmatrix} {e}_{\mu }\\ {e}_{\mathcal {M}} \end{pmatrix} + \begin{pmatrix} 0\\ 1 \end{pmatrix} {e}_{\textrm{d}},\\&\qquad {f}_{\textrm{d}}= \begin{pmatrix} 0&1 \end{pmatrix} \begin{pmatrix} {e}_{\mu }\\ {e}_{\mathcal {M}} \end{pmatrix}, \\&\quad \left. 0 = ({e}_{\mu } \mu + \iota _{{e}_{\mathcal {M}}} \mathcal {M})|_{{\partial \mathcal {S}}} \right\} , \end{aligned} \end{aligned}$$

(37)

with ${f}_{\textrm{k}} = ({f}_{\mu },{f}_{\mathcal {M}}) \in T_{\chi }{{\mathcal {X}}}_{\textrm{k}}$ and ${e}_{\textrm{k}} = ({e}_{\mu },{e}_{\mathcal {M}})\in T^*_{\chi }{{\mathcal {X}}}_{\textrm{k}}$. The Dirac structure (37) characterizes the power balance

$$\begin{aligned} \langle {e}_{\mu } | {f}_{\mu } \rangle _{\mathcal {S}} + \langle {e}_{\mathcal {M}} | {f}_{\mathcal {M}} \rangle _{\mathcal {S}} = \langle {e}_{\textrm{d}} | {f}_{\textrm{d}} \rangle _{\mathcal {S}}. \end{aligned}$$

(38)

Proof

The proof follows by deriving the expression of (36) from (24) using the results of Prop.2. Using the change of variables formula and the expressions of $\delta _{\varphi }{\tilde{\mathscr {G}}}$ and $\delta ^\cup _{\varphi }{\tilde{\mathscr {G}}}$ in Prop.2 (ii), one has that

$$\begin{aligned} \int _{\mathcal {B}}\delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}} \ \dot{\wedge }\ \delta _{\varphi }{\tilde{\mathscr {G}}} =&\int _{\mathcal {S}}\delta _{\mathcal {M}}{\mathscr {F}} \ \dot{\wedge }\ \mu \otimes (\textrm{d}\delta _{\mu }{\mathscr {G}} + \nabla \delta _{\mathcal {M}}{\mathscr {G}} \ \dot{\wedge }\ v^\flat ),\end{aligned}$$

(39)

$$\begin{aligned} \int _{\partial \mathcal {B}} \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}}|_{\partial \mathcal {B}} \ \dot{\wedge }\ \delta ^\cup _{\varphi }{\tilde{\mathscr {G}}} =&-\int _{{\partial \mathcal {S}}} \delta _{\mathcal {M}}{\mathscr {F}}|_{{\partial \mathcal {S}}} \ \dot{\wedge }\ (\delta _{\mu }{\mathscr {G}} \mu + \iota _{\delta _{\mathcal {M}}{\mathscr {G}}}\mathcal {M})|_{{\partial \mathcal {S}}}. \end{aligned}$$

(40)

Furthermore, using the change of variables formula and the steps in the proof of Prop.2 (ii), one can show that

$$\begin{aligned} & \int _{\mathcal {B}}\delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {G}}} \ \dot{\wedge }\ \delta _{\varphi }{\tilde{\mathscr {F}}} + \int _{\partial \mathcal {B}} \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {G}}}|_{\partial \mathcal {B}} \ \dot{\wedge }\ \delta ^\cup _{\varphi }{\tilde{\mathscr {F}}} \nonumber \\ & \quad = - \int _{\mathcal {S}}\delta _{\mu }{\mathscr {F}} \wedge \mathcal {L}_{\delta _{\mathcal {M}}{\mathscr {G}}}{\mu } + \delta _{\mathcal {M}}{\mathscr {F}} \ \dot{\wedge }\ \textrm{d}_{\nabla }\iota _{\delta _{\mathcal {M}}{\mathscr {G}}}\mathcal {M}, \end{aligned}$$

(41)

where (96) and (97) were used in reverse with $\delta _{\mathcal {M}}{\mathscr {G}}$ instead of $\delta \varphi$. Substituting (39-41) in (24), and using the Lie derivative identity (Gilbert and Vanneste 2023)

$$\begin{aligned} \mathcal {L}_{\eta }{\mathcal {M}} = \textrm{d}_{\nabla }(\iota _\eta \mathcal {M}) + \mu \otimes (\nabla \eta \ \dot{\wedge }\ v^\flat ), \qquad \forall \eta \in \varGamma (T\mathcal {S}), \end{aligned}$$

(42)

allow us to express the spatial Poisson bracket (36) as

$$\begin{aligned} ^s\{\mathscr {F},\mathscr {G}\} =&-\int _{\mathcal {S}}\delta _{\mu }{\mathscr {F}} \wedge \mathcal {L}_{\delta _{\mathcal {M}}{\mathscr {G}}}{\mu } + \delta _{\mathcal {M}}{\mathscr {F}} \ \dot{\wedge }\ \mu \otimes \textrm{d}\delta _{\mu }{\mathscr {G}} \\&\quad + \delta _{\mathcal {M}}{\mathscr {F}} \ \dot{\wedge }\ (\textrm{d}_{\nabla }\iota _{\delta _{\mathcal {M}}{\mathscr {G}}}\mathcal {M}+ \mu \otimes ( \nabla \delta _{\mathcal {M}}{\mathscr {G}} \ \dot{\wedge }\ v^\flat )) \\&\quad + \int _{{\partial \mathcal {S}}} \delta _{\mathcal {M}}{\mathscr {F}}|_{{\partial \mathcal {S}}} \ \dot{\wedge }\ (\delta _{\mu }{\mathscr {G}} \mu + \iota _{\delta _{\mathcal {M}}{\mathscr {G}}}\mathcal {M})|_{{\partial \mathcal {S}}}\\ =&-\int _{\mathcal {S}}\delta _{\mu }{\mathscr {F}} \wedge \textrm{d}\iota _{\delta _{\mathcal {M}}{\mathscr {G}}}{\mu } + \delta _{\mathcal {M}}{\mathscr {F}} \ \dot{\wedge }\ \mu \otimes \textrm{d}\delta _{\mu }{\mathscr {G}} + \delta _{\mathcal {M}}{\mathscr {F}} \ \dot{\wedge }\ \mathcal {L}_{\delta _{\mathcal {M}}{\mathscr {G}}}{\mathcal {M}}\\&\quad + \int _{{\partial \mathcal {S}}} \delta _{\mathcal {M}}{\mathscr {F}}|_{{\partial \mathcal {S}}} \ \dot{\wedge }\ (\delta _{\mu }{\mathscr {G}} \mu + \iota _{\delta _{\mathcal {M}}{\mathscr {G}}}\mathcal {M})|_{{\partial \mathcal {S}}}. \end{aligned}$$

Now by letting $\mathscr {F}=\mathscr {G}$ one has that $^s\{\mathscr {G},\mathscr {G}\} =0$ and by introducing $, {e}_{\mu }:= \delta _{\mu }{\mathscr {G}}, {e}_{\mathcal {M}}:= \delta _{\mathcal {M}}{\mathscr {G}}$, and

$${f}_{\mu }:=-\textrm{d}\iota _{{e}_{\mathcal {M}}}\mu ,\qquad {f}_{\mathcal {M}}:= - \mu \otimes \textrm{d}{e}_{\mu } - \mathcal {L}_{{e}_{\mathcal {M}}}{\mathcal {M}},\qquad ({e}_{\mu } \mu + \iota _{{e}_{\mathcal {M}}} \mathcal {M})|_{{\partial \mathcal {S}}} = 0,$$

then the above Poisson bracket expression becomes $\langle {e}_{\mu } | {f}_{\mu } \rangle _{\mathcal {S}} + \langle {e}_{\mathcal {M}} | {f}_{\mathcal {M}} \rangle _{\mathcal {S}} = 0$. Finally, by introducing ${f}_{\textrm{d}}:= \varphi _{\textrm{f},*}(\tilde{f}_{\text {d}}),{e}_{\textrm{d}}:= \varphi _{\textrm{f},*}({\tilde{e}}_{\textrm{d}})$, the interaction port $(\tilde{f}_{\text {d}},\tilde{e}_{\text {d}})$ in (28) can be written using the change of variables formula as

$$\int _{\mathcal {B}}\tilde{f}_{\text {d}}\ \dot{\wedge }\ \tilde{e}_{\text {d}} = \int _{\mathcal {S}}\varphi _{\textrm{f},*}(\tilde{f}_{\text {d}})\ \dot{\wedge }\ \varphi _{\textrm{f},*}({\tilde{e}}_{\textrm{d}}) = \int _{\mathcal {S}}{f}_{\textrm{d}}\ \dot{\wedge }\ {e}_{\textrm{d}},$$

which concludes the transformation of (28) into (37). $\square$

Proposition 3

For the kinetic energy as a Hamiltonian functional with state variables ${\chi }:= (\mu ,\mathcal {M}) \in {{\mathcal {X}}}_{\textrm{k}}$, the equations of motion in the spatial representation are given by:

$$\begin{aligned} \begin{pmatrix} \partial _t \mu \\ \partial _t \mathcal {M} \end{pmatrix}&= \begin{pmatrix} 0 & - \textrm{d}\iota _{(\cdot )} \mu \\ - \mu \otimes \textrm{d}(\cdot ) & - \mathcal {L}_{(\cdot )}{\mathcal {M}} \end{pmatrix} \begin{pmatrix} \delta _{\mu }{{{H}}_{\textrm{k}}}\\ \delta _{\mathcal {M}}{{{H}}_{\textrm{k}}} \end{pmatrix} + \begin{pmatrix} 0\\ 1 \end{pmatrix} {{\mathcal {F}}}, \end{aligned}$$

(43)

$$\begin{aligned} v&= \left( 0 \quad 1\right) \begin{pmatrix} \delta _{\mu }{{{H}}_{\textrm{k}}}\\ \delta _{\mathcal {M}}{{{H}}_{\textrm{k}}} \end{pmatrix}, \end{aligned}$$

(44)

with ${\mathcal {F}}\in \varOmega ^{n}(\mathcal {S};T^*\mathcal {S})$ denoting the external body force field in the spatial representation. The variational derivatives of the functional ${{H}}_{\textrm{k}}$ in (32) with respect to $\mu \in \varOmega ^{n}(\mathcal {S})$ and $\mathcal {M}\in \varOmega ^{n}(\mathcal {S};T^*\mathcal {S})$ are given, respectively, by

$$\delta _{\mu }{{{H}}_{\textrm{k}}} = - \frac{1}{2}\iota _{v}v^\flat \in \varOmega ^{0}(\mathcal {S}), \qquad \qquad \delta _{\mathcal {M}}{{{H}}_{\textrm{k}}} = v\in \varOmega ^{0}(\mathcal {S};T\mathcal {S}),$$

where $v= \star _c^{-1}\mathcal {M}$ is the spatial velocity field. Furthermore, the state variables are constrained on the boundary ${\partial \mathcal {S}}$ such that the extensive momentum flux vanishes, i.e., $\iota _v\mathcal {M}|_{{\partial \mathcal {S}}} = 0 \in \varOmega ^{n-1}(\partial \mathcal {S};T^*\mathcal {S})$.

The Hamiltonian functional (32) admits its rate of change such that along trajectories $(\mu (t),\mathcal {M}(t))$ of (43-44), it holds that

$$\begin{aligned} {\dot{{H}}}_{\textrm{k}}= \langle {\mathcal {F}} | v \rangle _{\mathcal {S}}. \end{aligned}$$

(45)

Proof

(i) The proof of (43-44) follows from the expression of ${{\mathcal {D}}}_{\textrm{k}}({\chi })$ in (37) by setting the inputs $({e}_{\mu },{e}_{\mathcal {M}}, {e}_{\textrm{d}}) = ({\delta _{\mu }{{{H}}_{\textrm{k}}}},{\delta _{\mathcal {M}}{{{H}}_{\textrm{k}}}},{\mathcal {F}})$ which leads to the outputs $({f}_{\mu },{f}_{\mathcal {M}}, {f}_{\textrm{d}}) = ({\partial _t \mu },{\partial _t \mathcal {M}},\delta _{\mathcal {M}}{{{H}}_{\textrm{k}}})$. Consequently, using (35) for ${{H}}_{\textrm{k}}$, the energy balance (38) becomes ${\dot{{H}}}_{\textrm{k}}= \langle {\mathcal {F}} | \delta _{\mathcal {M}}{{{H}}_{\textrm{k}}} \rangle _{\mathcal {S}}$.

(ii) The proof of the variational derivatives goes as follows.

First, the Hamiltonian (32) is rewritten to show explicitly its dependence on $\mu$ and $\mathcal {M}$ as

$${{H}}_{\textrm{k}}[\mu ,\mathcal {M}]:= \int _{\mathcal {S}}\frac{1}{2 \mu } g^{ij} \mathcal {M}_i\mathcal {M}_j,$$

where $g^{ij}\in \varOmega ^{0}(\mathcal {S})$ are the components of the inverse metric $g^{-1}$ and $\mathcal {M}_i\in \varOmega ^{n}(\mathcal {S})$ denotes the top-form components of $\mathcal {M}$. Now consider the tangent vector to the curve $s \mapsto \mu _s \in \varOmega ^{n}(\mathcal {S})$ denoted by $\delta \mu := \left. \frac{d}{ds}\right| _{s=0}\mu _s \in \varOmega ^{n}(\mathcal {S}),$ with $\mu _s|_{s=0} = \mu$. The variational derivative $\delta _{\mu }{{{H}}_{\textrm{k}}}\in \varOmega ^{0}(\mathcal {S})$ is defined such that

$$\left. \frac{d}{ds}\right| _{s=0}{{H}}_{\textrm{k}}[\mu _s,\mathcal {M}] = \langle \delta _{\mu }{{{H}}_{\textrm{k}}} | \delta \mu \rangle _{\mathcal {S}},$$

where the LHS can be expressed as

$$\left. \frac{d}{ds}\right| _{s=0}{{H}}_{\textrm{k}}[\mu _s,\mathcal {M}] = \int _{\mathcal {S}}\frac{1}{2} \left. \frac{d}{ds}\right| _{s=0}(\frac{1}{\mu }) g^{ij} \mathcal {M}_i\mathcal {M}_j.$$

Using the identity $\left. \frac{d}{ds}\right| _{s=0}(\frac{1}{\mu }) = - \frac{1}{\mu ^2}\left. \frac{d}{ds}\right| _{s=0}(\mu _s)$, the definition of $\delta \mu$ and the fact that $\mathcal {M}_i = g_{ij} v^j \mu$ leads to

$$\left. \frac{d}{ds}\right| _{s=0}{{H}}_{\textrm{k}}[\mu _s,\mathcal {M}] = - \int _{\mathcal {S}}\frac{1}{2} g_{km} v^k v^m \delta \mu = - \int _{\mathcal {S}}\frac{1}{2} g(v,v) \delta \mu .$$

Using exterior calculus, one has that $g(v,v) = \iota _{v}v^\flat \in \varOmega ^{0}(\mathcal {S})$, which concludes the derivation of $\delta _{\mu }{{{H}}_{\textrm{k}}}$. As for $\delta _{\mathcal {M}}{{{H}}_{\textrm{k}}} = v$, one can also follow similar steps in deriving it as above or alternatively use the relation $\delta _{\tilde{\mathcal {M}}}{{\tilde{H}}_{\textrm{k}}} = \varphi _{\textrm{f}}^*(\delta _{\mathcal {M}}{{{H}}_{\textrm{k}}})$ derived in Prop. 2. Consequently, using the expression of $\delta _{\tilde{\mathcal {M}}}{{\tilde{H}}_{\textrm{k}}}$ it follows that $\delta _{\mathcal {M}}{{{H}}_{\textrm{k}}} = v$.

(iii) Finally, the boundary constraint in (37) can be rewritten using the expressions of the variational derivatives as

$$\begin{aligned} 0&= (\delta _{\mu }{{{H}}_{\textrm{k}}} \mu + \iota _{\delta _{\mathcal {M}}{{{H}}_{\textrm{k}}}} \mathcal {M})|_{{\partial \mathcal {S}}} = (-\frac{1}{2}\iota _vv^\flat \mu + \iota _v\mu \otimes v^\flat )|_{{\partial \mathcal {S}}}\\&= (-\frac{1}{2}\iota _v\mu \otimes v^\flat + \iota _v\mu \otimes v^\flat )|_{{\partial \mathcal {S}}} = (\frac{1}{2}\iota _v\mu \otimes v^\flat )|_{{\partial \mathcal {S}}} = \frac{1}{2}\iota _v\mathcal {M}|_{{\partial \mathcal {S}}}, \end{aligned}$$

which concludes the proof. $\square$

Corollary 2

The port-Hamiltonian dynamic equations (43) incorporate the conservation of mass and momentum laws which can be expressed, respectively, in:

(i) Advection formulation:

$$\partial _t \mu = -\mathcal {L}_{v}{\mu }, \qquad \partial _t \mathcal {M}= -\mathcal {L}_{v}{\mathcal {M}} + \frac{1}{2}\mu \otimes \textrm{d}|v|_{g}^2 + {\mathcal {F}},$$

with $|v|_{g}:= \sqrt{g(v,v)}$ denoting the norm with respect to $g$.

(ii) Conservation formulation:

$$\partial _t \mu = -\textrm{d}\iota _{v}{\mu }, \qquad \partial _t \mathcal {M}= -\textrm{d}_{\nabla }\iota _{v}\mathcal {M}+ {\mathcal {F}}.$$

Proof

The proof of the mass balance in (ii) and the momentum balance in (i) follows immediately from substituting the expressions of the variational derivatives, whereas the mass balance in (i) is derived using Cartan’s homotopy formula $\mathcal {L}_{v} = \textrm{d}\iota _v + \iota _v\textrm{d}$, while the momentum balance in (ii) follows from using the identity $\nabla v\ \dot{\wedge }\ v^\flat = \frac{1}{2}\textrm{d}|v|_{g}^2$ and the Lie derivative identity (42). $\square$

3.2.4 Convective Dirac Structure

Finally, we conclude by presenting the convective representation of the kinetic energy subsystem which follows exactly the same line of thought presented above for the reduction to the spatial representation. For compactness purposes, the proofs of this section will be shorter versions of their counterparts in the spatial representation.

The transformation of the material variables to the convective ones is carried out using the nonlinear diffeomorphism

$$\begin{aligned} \begin{aligned} ^c\varPhi : {{\tilde{\mathcal {X}}}_{\textrm{k}}\times {\tilde{\mathcal {P}}}_{\textrm{k}}}&\rightarrow {{\hat{\mathcal {X}}}_{\textrm{k}}\times {\hat{\mathcal {P}}}_{\textrm{k}}}\\ {((\varphi ,\tilde{\mathcal {M}}),(g,\tilde{\mu }))}&\mapsto {((\hat{g},\hat{\mathcal {M}}),\hat{\mu })} \end{aligned} \end{aligned}$$

(46)

with $\hat{g}:= \varphi ^* (g) \in \mathcal {M}(\mathcal {B})$ denoting the convective metric and $\hat{\mathcal {M}}:= \varphi _{\textrm{v}}^*(\tilde{\mathcal {M}}) \in \varOmega ^{n}(\mathcal {B};T^*\mathcal {B})$ denoting the extensive momentum in the convective representation, while $\hat{\mu }= \tilde{\mu }\in \varOmega ^{n}(\mathcal {B})$. We denote the convective state space by ${\hat{\mathcal {X}}}_{\textrm{k}}:= \mathcal {M}(\mathcal {B})\times \varOmega ^{n}(\mathcal {B};T^*\mathcal {B})$ and the convective parameter space by ${\hat{\mathcal {P}}}_{\textrm{k}}:=\varOmega ^{n}(\mathcal {B})$.

The convective kinetic energy Hamiltonian ${{\hat{H}}_{\textrm{k}}}:{{\hat{\mathcal {X}}}_{\textrm{k}}\times {\hat{\mathcal {P}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ is defined by ${\hat{H}}_{\textrm{k}}[\hat{g},\hat{\mathcal {M}};\hat{\mu }] = {\tilde{H}}_{\textrm{k}}[\varphi ,\tilde{\mathcal {M}};g,\tilde{\mu }]$ and is expressed as

$$\begin{aligned} {\hat{H}}_{\textrm{k}}[\hat{g},\hat{\mathcal {M}}]:= \int _{\mathcal {B}}\frac{1}{2}\hat{\star }_c^{-1}\hat{\mathcal {M}}\ \dot{\wedge }\ \hat{\mathcal {M}}, \end{aligned}$$

(47)

where the convective Hodge star ${\hat{\star }_c}:{\varOmega ^{k}(\mathcal {B};T\mathcal {B})}\rightarrow {\varOmega ^{n-k}(\mathcal {B};T^*\mathcal {B})}$ incorporates a state dependency on the convective metric $\hat{g}$ and a parametric dependency on the mass form $\hat{\mu }$.

Remark 4

In local coordinates, we have that the convective metric expressed as $\hat{g}= F^i_I F^j_J g_{ij}\circ \varphi E^I \otimes E^J$, with $F^i_I\in \varOmega ^{0}(\mathcal {B})$ denoting the local components of the deformation gradient $F:= T\varphi$. The convective momentum is computed from its material counterpart (cf. . Remark 1) by

$$\hat{\mathcal {M}}= \varphi _{\textrm{v}}^*(\tilde{\mathcal {M}})= \tilde{v}_i \tilde{\mu }\otimes \varphi ^*{(e^i)} = \tilde{v}_i \tilde{\mu }\otimes F^i_I E^I = {\hat{v}}_I \hat{\mu }\otimes E^I,$$

and thus, $\hat{\mathcal {M}}$ is identified with the trivial covector-valued form $\hat{\mu }\otimes \hat{v}^\flat$ similar to the spatial momentum $\mathcal {M}$.

In general, we denote by ${\hat{\mathscr {F}}:= \tilde{\mathscr {F}} \circ ^c\varPhi ^{-1}}:{{\hat{\mathcal {X}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ the convective counterpart of any material functional ${\tilde{\mathscr {F}}}:{{\tilde{\mathcal {X}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ defined such that:

$$\hat{\mathscr {F}}[\varphi ^* (g),\varphi _{\textrm{v}}^*(\tilde{\mathcal {M}})] = \tilde{\mathscr {F}}[\varphi ,\tilde{\mathcal {M}}].$$

In our work, we identify tangent and cotangent spaces of $\mathcal {M}(\mathcal {B})$ by Rashad et al. (2023) $T_{\hat{g}}\mathcal {M}(\mathcal {B})\cong \varOmega _{\textrm{sym}}^{1}(\mathcal {B};T^*\mathcal {B})\subset \varOmega ^{1}(\mathcal {B};T^*\mathcal {B})$ and $T_{\hat{g}}^*\mathcal {M}(\mathcal {B})\cong \varOmega _{\textrm{sym}}^{n-1}(\mathcal {B};T\mathcal {B})\subset \varOmega ^{n-1}(\mathcal {B};T\mathcal {B})$. The restriction to symmetric subspaces comes from the symmetric nature of the Riemannian metric $\hat{g}$. More details on the technical constructions of these spaces are discussed in Sect. 3.3. The variational derivatives of the functional $\hat{\mathscr {F}}$ with respect to $(\hat{g},\hat{\mathcal {M}})$ are the bundle-valued forms $\delta _{\hat{g}}{\hat{\mathscr {F}}} \in \varOmega _{\textrm{sym}}^{n-1}(\mathcal {B};T\mathcal {B})$ and $\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}} \in \varOmega ^{0}(\mathcal {B};T\mathcal {B}),$ that satisfy for any $\delta \hat{g}\in \varOmega _{\textrm{sym}}^{1}(\mathcal {B};T^*\mathcal {B}),\delta \hat{\mathcal {M}}\in \varOmega ^{n}(\mathcal {B};T^*\mathcal {B})$

$$ \textrm{D}_{\hat{g}}\hat{\mathscr {F}}(\hat{g},\hat{\mathcal {M}}) \cdot \delta \hat{g}= \int _{\mathcal {B}}\delta \hat{g}\ \dot{\wedge }\ \delta _{\hat{g}}{\hat{\mathscr {F}}}, \qquad \qquad \textrm{D}_{\hat{\mathcal {M}}}\hat{\mathscr {F}}(\hat{g},\hat{\mathcal {M}}) \cdot \delta \hat{\mathcal {M}}= \int _{\mathcal {B}}\delta \hat{\mathcal {M}}\ \dot{\wedge }\ \delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}}.$$

The rate of change of $\hat{\mathscr {F}}$ with respect to time is then expressed as

$$\begin{aligned} \dot{\hat{\mathscr {F}}} = \langle \delta _{\hat{g}}{\hat{\mathscr {F}}} | \partial _t \hat{g} \rangle _{\mathcal {B}} + \langle \delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}} | \partial _t \hat{\mathcal {M}} \rangle _{\mathcal {B}}, \end{aligned}$$

(48)

with $\partial _t \hat{g}\in \varOmega _{\textrm{sym}}^{1}(\mathcal {B};T^*\mathcal {B})$ and $\partial _t \hat{\mathcal {M}}\in \varOmega ^{n}(\mathcal {B};T^*\mathcal {B})$.

Let ${\varphi _{\textrm{v}}^*}:{\varOmega _\varphi ^{k}(\mathcal {B};T\mathcal {S})}\rightarrow {\varOmega ^{k}(\mathcal {B};T\mathcal {B})}$ denote the pullback of the value part of a vector-valued form, transforming it from the material representation to the convective one, and $\varphi _{\textrm{v},*}$ denote its inverse. We shall also use the same notation for covector-valued forms. Then the convective counterpart of Prop. 2 can be stated as follows.

Proposition 4

Let $(\varphi ,\tilde{\mathcal {M}}) \in {\tilde{\mathcal {X}}}_{\textrm{k}}$ and $(\hat{g},\hat{\mathcal {M}})\in {\hat{\mathcal {X}}}_{\textrm{k}}$ be related by the diffeomorphism $^c\varPhi$ in (46).

i) The pushforward map of $^c\varPhi$ is given by

$$\begin{aligned} \begin{aligned} ^c\varPhi _* : {T{\tilde{\mathcal {X}}}_{\textrm{k}}}&\rightarrow {T{\hat{\mathcal {X}}}_{\textrm{k}}}\\ {(\varphi ,\tilde{\mathcal {M}},\delta \tilde{\varphi },\delta \tilde{\mathcal {M}},\delta \tilde{\varphi }|_{{\partial \mathcal {B}}})}&\mapsto {(\hat{g},\hat{\mathcal {M}},\delta \hat{g},\delta \hat{\mathcal {M}})} \end{aligned} \end{aligned}$$

$$\delta \hat{g}= \mathcal {L}_{\delta \hat{\varphi }}{\hat{g}} \in \varOmega _{\textrm{sym}}^{1}(\mathcal {B};T\mathcal {B}),\qquad \delta \hat{\mathcal {M}}= \hat{\mu }\otimes (\hat{\nabla }\delta \hat{\varphi }\ \dot{\wedge }\ \hat{v}^\flat ) + \varphi _{\textrm{v}}^*(\delta \tilde{\mathcal {M}}) \in \varOmega ^{n}(\mathcal {B};T^*\mathcal {B}),$$

with $\delta \hat{\varphi }:= \varphi _{\textrm{v}}^*(\delta \tilde{\varphi }) = T \varphi ^{-1} \circ \delta \tilde{\varphi }\in \varOmega ^{0}(\mathcal {B};T\mathcal {B}) \cong \varGamma (T\mathcal {B})$.

ii) The pullback map of $^c\varPhi$ is given by

$$\begin{aligned} \begin{aligned} ^c\varPhi ^* : {T^*{\hat{\mathcal {X}}}_{\textrm{k}}}&\rightarrow {T^*{\tilde{\mathcal {X}}}_{\textrm{k}}}\\ {(\hat{g},\hat{\mathcal {M}},\delta _{\hat{g}}{\hat{\mathscr {F}}},\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}})}&\mapsto {(\varphi ,\tilde{\mathcal {M}},\delta _{\varphi }{\tilde{\mathscr {F}}}, \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}}, \delta ^\cup _{\varphi }{\tilde{\mathscr {F}}})} \end{aligned} \end{aligned}$$

with

$$\begin{aligned} \delta _{\varphi }{\tilde{\mathscr {F}}} =&- \varphi _{\textrm{v},*} (\hat{\textrm{d}}_{\hat{\nabla }}(2 (\delta _{\hat{g}}{\hat{\mathscr {F}}}) ^\flat + \iota _{\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}}}\hat{\mathcal {M}})) & \in \varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S}) \\ \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}} =&\varphi _{\textrm{v},*}(\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}}) & \in \varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})\\ \delta ^\cup _{\varphi }{\tilde{\mathscr {F}}} =&\varphi _{\textrm{v},*}(2 (\delta _{\hat{g}}{\hat{\mathscr {F}}}) ^\flat + \iota _{\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}}}\hat{\mathcal {M}})|_{{\partial \mathcal {B}}} & \in \varOmega _\varphi ^{n-1}(\partial \mathcal {B};T^*\mathcal {S}). \end{aligned}$$

Proof

See Appendix A.2. $\square$

Theorem 2

The convective Dirac structure corresponding to (28) under the diffeomorphism (46) is the sub-bundle ${\hat{\mathcal {D}}}_{\textrm{k}}\subset {\hat{\mathfrak {F}}}_{\textrm{k}}\oplus {\hat{\mathfrak {F}}^*}_{\textrm{k}}$ with

$$\begin{aligned} {\hat{\mathfrak {F}}}_{\textrm{k}}:=& ^c \varPhi _*(T{\tilde{\mathcal {X}}}_{\textrm{k}}) \times \varphi _{\textrm{v}}^*(\varOmega _\varphi ^{0}(\mathcal {B};T\mathcal {S})) = T{\hat{\mathcal {X}}}_{\textrm{k}}\times \varOmega ^{0}(\mathcal {B};T\mathcal {B})\\ {\hat{\mathfrak {F}}^*}_{\textrm{k}}:=&( ^c \varPhi ^*)^{-1}(T^*{\tilde{\mathcal {X}}}_{\textrm{k}}) \times \varphi _{\textrm{v}}^*(\varOmega _\varphi ^{n}(\mathcal {B};T^*\mathcal {S})) = T^*{\hat{\mathcal {X}}}_{\textrm{k}}\times \varOmega ^{n}(\mathcal {B};T^*\mathcal {B}), \end{aligned}$$

whose fiber at any $\hat{\chi }:= (\hat{g},\hat{\mathcal {M}})\in {\hat{\mathcal {X}}}_{\textrm{k}}$ is defined by

$$\begin{aligned} \begin{aligned} {\hat{\mathcal {D}}}_{\textrm{k}}(\hat{\chi }):= \{ (({\hat{f}}_{\textrm{k}},&{\hat{f}}_{\textrm{d}}),({\hat{e}}_{\textrm{k}},{\hat{e}}_{\textrm{d}} )) \in {\hat{\mathfrak {F}}}_{\textrm{k}}(\hat{\chi })\times {\hat{\mathfrak {F}}^*}_{\textrm{k}}(\hat{\chi }) | \\ \begin{pmatrix} \hat{f}_{\hat{g}}\\ \hat{f}_{\hat{\mathcal {M}}} \end{pmatrix}&= \begin{pmatrix} 0 & 2\ \textrm{sym}\circ \hat{g}\circ \hat{\textrm{d}}_{\hat{\nabla }} \\ 2\ \hat{\textrm{d}}_{\hat{\nabla }}\circ \hat{g} & \mathcal {L}_{(\cdot )}{\hat{\mathcal {M}}} \end{pmatrix} \begin{pmatrix} \hat{e}_{\hat{g}}\\ \hat{e}_{\hat{\mathcal {M}}} \end{pmatrix} + \begin{pmatrix} 0\\ 1 \end{pmatrix} {\hat{e}}_{\textrm{d}},\\ {\hat{f}}_{\textrm{d}}&= \begin{pmatrix} 0&1 \end{pmatrix} \begin{pmatrix} \hat{e}_{\hat{g}}\\ \hat{e}_{\hat{\mathcal {M}}} \end{pmatrix}, \\ 0&= (2 \hat{g}\cdot \hat{e}_{\hat{g}} + \iota _{\hat{e}_{\hat{\mathcal {M}}}}\hat{\mathcal {M}})|_{{\partial \mathcal {B}}} \}, \end{aligned} \end{aligned}$$

(49)

with ${\hat{f}}_{\textrm{k}} = (\hat{f}_{\hat{g}},\hat{f}_{\hat{\mathcal {M}}}) \in T_{\hat{\chi }}{\hat{\mathcal {X}}}_{\textrm{k}}$ and ${\hat{e}}_{\textrm{k}} = (\hat{e}_{\hat{g}},\hat{e}_{\hat{\mathcal {M}}})\in T^*_{\hat{\chi }}{\hat{\mathcal {X}}}_{\textrm{k}}$. The Dirac structure (49) characterizes the power balance

$$\begin{aligned} \langle \hat{e}_{\hat{g}} | \hat{f}_{\hat{g}} \rangle _{\mathcal {B}} + \langle \hat{e}_{\hat{\mathcal {M}}} | \hat{f}_{\hat{\mathcal {M}}} \rangle _{\mathcal {B}} = \langle {\hat{e}}_{\textrm{d}} | {\hat{f}}_{\textrm{d}} \rangle _{\mathcal {B}}. \end{aligned}$$

(50)

Proof

Using the expressions of $\delta _{\varphi }{\tilde{\mathscr {G}}}$ and $\delta ^\cup _{\varphi }{\tilde{\mathscr {G}}}$ in Prop.4 (ii) in addition to the duality of $\varphi _{\textrm{v},*}$ and $\varphi _{\textrm{v}}^*$, one can rewrite the different integrals in (24) as

$$\begin{aligned} & \int _{\mathcal {B}}\delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}} \ \dot{\wedge }\ \delta _{\varphi }{\tilde{\mathscr {G}}} + \int _{\partial \mathcal {B}} \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}}|_{\partial \mathcal {B}} \ \dot{\wedge }\ \delta ^\cup _{\varphi }{\tilde{\mathscr {G}}}\nonumber \\ & \quad = -\int _{\mathcal {B}}\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}} \ \dot{\wedge }\ \hat{\textrm{d}}_{\hat{\nabla }}\mathcal {E}_{\hat{\mathscr {G}}} + \int _{\partial \mathcal {B}} \delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}}|_{{\partial \mathcal {B}}} \ \dot{\wedge }\ \mathcal {E}_{\hat{\mathscr {G}}}|_{{\partial \mathcal {B}}}, \end{aligned}$$

(51)

where $\mathcal {E}_{\hat{\mathscr {G}}}:= 2 (\delta _{\hat{g}}{\hat{\mathscr {G}}}) ^\flat + \iota _{\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {G}}}}\hat{\mathcal {M}}\in \varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B})$. Furthermore, using the steps in the proof of Prop.4 (ii) in reverse, one can show that

$$\begin{aligned} \int _{\mathcal {B}}\delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {G}}} \ \dot{\wedge }\ \delta _{\varphi }{\tilde{\mathscr {F}}} +&\int _{\partial \mathcal {B}} \delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {G}}}|_{\partial \mathcal {B}} \ \dot{\wedge }\ \delta ^\cup _{\varphi }{\tilde{\mathscr {F}}} = \int _{\mathcal {B}}\hat{\nabla }\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {G}}} \ \dot{\wedge }\ (2 \hat{g}\delta _{\hat{g}}{\hat{\mathscr {F}}} + \iota _{\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}}}\hat{\mathcal {M}})\nonumber \\ =&\int _{\mathcal {B}}\delta _{\hat{g}}{\hat{\mathscr {F}}} \ \dot{\wedge }\ 2 \textrm{sym}(\hat{g}\cdot \hat{\nabla }\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {G}}}) + \delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}} \ \dot{\wedge }\ \hat{\mu }\otimes (\hat{\nabla }\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {G}}} \ \dot{\wedge }\ \hat{v}^\flat ). \end{aligned}$$

(52)

Let $^c\{\hat{\mathscr {F}},\hat{\mathscr {G}}\}$ denote the convective representation of the canonical Poisson bracket defined similar to (36). By substituting (51,52) in (24) and using the convective version of (42), one can express this convective Poisson bracket as

$$\begin{aligned} ^c\{\hat{\mathscr {F}},\hat{\mathscr {G}}\} =&\int _{\mathcal {B}}\delta _{\hat{g}}{\hat{\mathscr {F}}} \ \dot{\wedge }\ 2 \textrm{sym}(\hat{g}\cdot \hat{\nabla }\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {G}}}) + \delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}} \ \dot{\wedge }\ \hat{\mu }\otimes (\hat{\nabla }\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {G}}} \ \dot{\wedge }\ \hat{v}^\flat ) \\&+ \delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}} \ \dot{\wedge }\ \hat{\textrm{d}}_{\hat{\nabla }}((\delta _{\hat{g}}{\hat{\mathscr {G}}}) ^\flat + \iota _{\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {G}}}}\hat{\mathcal {M}}) - \int _{\partial \mathcal {B}} \delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}}|_{{\partial \mathcal {B}}} \ \dot{\wedge }\ \mathcal {E}_{\hat{\mathscr {G}}}|_{{\partial \mathcal {B}}}\\ =&\int _{\mathcal {B}}\delta _{\hat{g}}{\hat{\mathscr {F}}} \ \dot{\wedge }\ 2 \textrm{sym}(\hat{g}\cdot \hat{\nabla }\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {G}}}) +\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}} \ \dot{\wedge }\ (\hat{\textrm{d}}_{\hat{\nabla }}(\delta _{\hat{g}}{\hat{\mathscr {G}}}) ^\flat + \mathcal {L}_{\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {G}}}}{\hat{\mathcal {M}}}) \\&- \int _{\partial \mathcal {B}} \delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}}|_{{\partial \mathcal {B}}} \ \dot{\wedge }\ (2 (\delta _{\hat{g}}{\hat{\mathscr {G}}}) ^\flat + \iota _{\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {G}}}}\hat{\mathcal {M}})|_{{\partial \mathcal {B}}}. \end{aligned}$$

Finally, following the same line of thought of the proof of Th. 1, one can transform the Poisson bracket $^c\{\hat{\mathscr {F}},\hat{\mathscr {G}}\}$ into (49). $\square$

Proposition 5

For the kinetic energy as a Hamiltonian functional with state variables $\hat{\chi }:= (\hat{g},\hat{\mathcal {M}}) \in {\hat{\mathcal {X}}}_{\textrm{k}}$, the equations of motion in the convective representation are given by:

$$\begin{aligned} \begin{pmatrix} \partial _t \hat{g}\\ \partial _t \hat{\mathcal {M}} \end{pmatrix}&= \begin{pmatrix} 0 & 2\ \textrm{sym}\circ \hat{g}\circ \hat{\textrm{d}}_{\hat{\nabla }} \\ 2\ \hat{\textrm{d}}_{\hat{\nabla }}\circ \hat{g} & \mathcal {L}_{(\cdot )}{\hat{\mathcal {M}}} \end{pmatrix} \begin{pmatrix} \delta _{\hat{g}}{{\hat{H}}_{\textrm{k}}}\\ \delta _{\hat{\mathcal {M}}}{{\hat{H}}_{\textrm{k}}} \end{pmatrix}+ \begin{pmatrix} 0\\ 1 \end{pmatrix} {\hat{\mathcal {F}}}, \end{aligned}$$

(53)

$$\begin{aligned} {\hat{v}}&= \begin{pmatrix} 0&1 \end{pmatrix} \begin{pmatrix} \delta _{\hat{g}}{{\hat{H}}_{\textrm{k}}}\\ \delta _{\hat{\mathcal {M}}}{{\hat{H}}_{\textrm{k}}} \end{pmatrix}. \end{aligned}$$

(54)

with $\hat{\mathcal {F}}\in \varOmega ^{n}(\mathcal {B};T^*\mathcal {B})$ denoting the external body force field in the convective representation. The variational derivatives of the functional ${\hat{H}}_{\textrm{k}}$ in (47) with respect to $\hat{g}\in \mathcal {M}(\mathcal {B})$ and $\hat{\mathcal {M}}\in \varOmega ^{n}(\mathcal {B};T^*\mathcal {B})$ are given, respectively, by

$$\delta _{\hat{g}}{{\hat{H}}_{\textrm{k}}} = - \frac{1}{2}\iota _{{\hat{v}}}\hat{\mu }\otimes {\hat{v}}\in \varOmega ^{n-1}(\mathcal {B};T\mathcal {B}), \qquad \qquad \delta _{\hat{\mathcal {M}}}{{\hat{H}}_{\textrm{k}}} = {\hat{v}}\in \varOmega ^{0}(\mathcal {B};T\mathcal {B}),$$

where ${\hat{v}}= \hat{\star }_c^{-1}\hat{\mathcal {M}}$ is the convective velocity field.

The Hamiltonian functional (47) admits its rate of change such that along trajectories $(\hat{g}(t),\hat{\mathcal {M}}(t))$ of (53-54), it holds that

$$\begin{aligned} {\dot{\hat{H}}}_{\textrm{k}}= \langle \hat{\mathcal {F}} | {\hat{v}} \rangle _{\mathcal {B}}. \end{aligned}$$

(55)

Proof

i) The proof of (53-54) follows from (49) by setting the inputs $(\hat{e}_{\hat{g}},\hat{e}_{\hat{\mathcal {M}}}, {\hat{e}}_{\textrm{d}}) = ({\delta _{\hat{g}}{{\hat{H}}_{\textrm{k}}}},{\delta _{\hat{\mathcal {M}}}{{\hat{H}}_{\textrm{k}}}},\hat{\mathcal {F}})$ which leads to the outputs $(\hat{f}_{\hat{g}},\hat{f}_{\hat{\mathcal {M}}}, {\hat{f}}_{\textrm{d}}) = ({\partial _t \hat{g}},{\partial _t \hat{\mathcal {M}}},\delta _{\hat{\mathcal {M}}}{{\hat{H}}_{\textrm{k}}})$. Consequently, using (48), ${\dot{\hat{H}}}_{\textrm{k}}= \langle \hat{\mathcal {F}} | \delta _{\hat{\mathcal {M}}}{{\hat{H}}_{\textrm{k}}} \rangle _{\mathcal {B}}$ follows from (50).

ii) The Hamiltonian (47) can be rewritten to show explicitly its dependence on $\hat{g}$ and $\hat{\mathcal {M}}$ as

$${\hat{H}}_{\textrm{k}}[\hat{g},\hat{\mathcal {M}}]:= \int _{\mathcal {B}}\frac{1}{2 \hat{\mu }} \hat{g}^{IJ} \hat{\mathcal {M}}_I\hat{\mathcal {M}}_J,$$

where $\hat{g}^{IJ}\in \varOmega ^{0}(\mathcal {B})$ are the components of the inverse metric $\hat{g}^{-1}$ such that $\hat{g}^{IJ}\hat{g}_{JK} = \delta ^I_K$, with $\hat{g}= \hat{g}_{JK} E^J\otimes E^K$. Furthermore, $\hat{\mathcal {M}}_J\in \varOmega ^{n}(\mathcal {B})$ denotes the top-form components of $\hat{\mathcal {M}}$ such that $\hat{\mathcal {M}}= \hat{\mathcal {M}}_J\otimes E^J$.

Consider the tangent vector to the curve $s \mapsto \hat{g}_s \in \mathcal {M}(\mathcal {B})$ denoted by

$${\delta \hat{g}}:= \left. \frac{d}{ds}\right| _{s=0}\hat{g}_s \in T_{\hat{g}}\mathcal {M}(\mathcal {B})\cong \varOmega ^{1}(\mathcal {B};T^*\mathcal {B}),$$

with $\hat{g}_s|_{s=0} = \hat{g}$. The variational derivative $\delta _{\hat{g}}{{\hat{H}}_{\textrm{k}}}\in \varOmega ^{n-1}(\mathcal {B};T\mathcal {B})$ is then defined implicitly such that

$$\left. \frac{d}{ds}\right| _{s=0}{\hat{H}}_{\textrm{k}}[\hat{g}_s,\hat{\mathcal {M}}] = \langle \delta _{\hat{g}}{{\hat{H}}_{\textrm{k}}} | {\delta \hat{g}} \rangle _{\mathcal {B}}.$$

Using the expression of ${\hat{H}}_{\textrm{k}}$ in (i), one can express the LHS as

$$\left. \frac{d}{ds}\right| _{s=0}{\hat{H}}_{\textrm{k}}[\hat{g}_s,\hat{\mathcal {M}}] = \int _{\mathcal {B}}\frac{1}{2 \hat{\mu }} \left. \frac{d}{ds}\right| _{s=0}(\hat{g}_s^{IJ}) \hat{\mathcal {M}}_I\hat{\mathcal {M}}_J.$$

Using the identity $\left. \frac{d}{ds}\right| _{s=0}(\hat{g}_s^{IJ}) = - \hat{g}^{IK}\hat{g}^{JL}\left. \frac{d}{ds}\right| _{s=0}((\hat{g}_s)_{KL})$, the definition of ${\delta \hat{g}}$ and the fact that $\hat{\mathcal {M}}_I = \hat{g}_{IJ} {\hat{v}}^J \hat{\mu }$ leads to

$$\left. \frac{d}{ds}\right| _{s=0}{\hat{H}}_{\textrm{k}}[\hat{g}_s,\hat{\mathcal {M}}] = - \int _{\mathcal {B}}\frac{1}{2 \hat{\mu }} \hat{g}^{IK}\hat{g}^{JL} \hat{\mathcal {M}}_I \hat{\mathcal {M}}_J {\delta \hat{g}}_{KL} = - \int _{\mathcal {B}}\frac{1}{2}{\hat{v}}^K {\hat{v}}^L {\delta \hat{g}}_{KL} \hat{\mu }.$$

The above expression can be written in tensor and exterior, respectively, as

$$\left. \frac{d}{ds}\right| _{s=0}{\hat{H}}_{\textrm{k}}[\hat{g}_s,\hat{\mathcal {M}}] = - \int _{\mathcal {B}}(\frac{1}{2}{\hat{v}}\otimes {\hat{v}}):{\delta \hat{g}}\ \hat{\mu }= - \int _{\mathcal {B}}{\delta \hat{g}}\ \dot{\wedge }\ \hat{\star }_c(\frac{1}{2}{\hat{v}}\otimes {\hat{v}}).$$

Using the identity $\hat{\star }_c({\hat{v}}\otimes {\hat{v}}) = \iota _{{\hat{v}}}\hat{\mu }\otimes {\hat{v}}$ concludes the derivation of $\delta _{\hat{g}}{{\hat{H}}_{\textrm{k}}}$. As for $\delta _{\hat{\mathcal {M}}}{{\hat{H}}_{\textrm{k}}} = {\hat{v}}$, it is straightforward following the steps in proof of Prop. 3.

(iii) Finally, the boundary constraint in (49) can be rewritten using the expressions of the variational derivatives as

$$\begin{aligned} 0 = (2 (\delta _{\hat{g}}{{\hat{H}}_{\textrm{k}}})^\flat + \iota _{{\hat{v}}} \hat{\mathcal {M}})|_{{\partial \mathcal {B}}} = (\iota _{{\hat{v}}}\hat{\mu }\otimes \hat{v}^\flat - \iota _{\hat{v}}\hat{\mathcal {M}})|_{{\partial \mathcal {B}}}, \end{aligned}$$

which is naturally satisfied by the state variables since $\iota _{{\hat{v}}}\hat{\mu }\otimes \hat{v}^\flat = \iota _{\hat{v}}\hat{\mathcal {M}}$. $\square$

Corollary 3

The port-Hamiltonian dynamic equations (53) incorporate the convective metric advection and conservation of momentum laws which can be expressed, respectively, in:

(i) Advection formulation:

$$\partial _t \hat{g}= \mathcal {L}_{{\hat{v}}}{\hat{g}}, \qquad \partial _t \hat{\mathcal {M}}= \mathcal {L}_{{\hat{v}}}{\hat{\mathcal {M}}} - \hat{\textrm{d}}_{\hat{\nabla }}\iota _{{\hat{v}}}\hat{\mathcal {M}}+ \hat{\mathcal {F}}.$$

(i) Conservation formulation:

$$\partial _t \hat{g}= 2\ \textrm{sym}(\hat{\nabla }\hat{v}^\flat ), \qquad \partial _t \hat{\mathcal {M}}= \frac{1}{2}\hat{\mu }\otimes \textrm{d}|{\hat{v}}|_{\hat{g}}^2 + \hat{\mathcal {F}},$$

with $|{\hat{v}}|_{\hat{g}}:= \sqrt{\hat{g}({\hat{v}},{\hat{v}})}$ denoting the norm with respect to $\hat{g}$.

Proof

The proof of the advection equation for $\hat{g}$ in (ii) and the momentum balance in (i) follows immediately from substituting the expressions of the variational derivatives and using the fact that $\iota _{\hat{v}}\hat{\mathcal {M}}= \iota _{{\hat{v}}}\hat{\mu }\otimes \hat{v}^\flat$, whereas the advection equation in (i) follows from the identity $\frac{1}{2}\mathcal {L}_{{\hat{v}}}{\hat{g}} = \textrm{sym}(\hat{\nabla }\hat{v}^\flat )$, while the momentum balance in (ii) follows from the identity $\hat{\nabla }{\hat{v}}\ \dot{\wedge }\ \hat{v}^\flat = \frac{1}{2}\textrm{d}|{\hat{v}}|_{\hat{g}}^2$ and the convective counterpart of (42). $\square$

Remark 5

We conclude this section by a number of remarks:

(i) It is important to note that the expression of ${\tilde{\mathcal {D}}}_{\textrm{k}}(\tilde{\chi })$ in (28) is constant and independent from the base point $\tilde{\chi }$ which is not the case for ${{\mathcal {D}}}_{\textrm{k}}({\chi })$ in (37) and ${\hat{\mathcal {D}}}_{\textrm{k}}(\hat{\chi })$ in (49). Thus, in contrast to its material counterpart ${\tilde{\mathcal {D}}}_{\textrm{k}}$, the spatial and convective Dirac structures are said to be modulated by the states ${\chi }= (\mu ,\mathcal {M})$ and $\hat{\chi }= (\hat{g},\hat{\mathcal {M}})$, respectively, which is a consequence of the Hamiltonian reduction procedure.

(ii) Another interesting observation is that in the spatial dynamics of nonlinear elasticity in Prop. 3, the momentum flux $\iota _{v}\mathcal {M}\in \varOmega ^{n-1}(\mathcal {S};T^*\mathcal {S})$ is constrained to vanish on the boundary ${\partial \mathcal {S}}:= \varphi ({\partial \mathcal {B}})$ of the embedded body in the ambient space, which is unlike the convective dynamics in Prop.5. This is a consequence of the fact that the boundary constraint in (49) is naturally satisfied by the kinetic energy Hamiltonian (47) because it naturally encloses the matter particles for all time. It would be interesting to investigate further the implications of such difference between the convective and spatial representation on the analysis of these dynamic equations.

3.3 Stress Power Subsystem

Fig. 3

Convective port-Hamiltonian model of the stress power Stokes–Dirac structure expressed in terms of bond graphs (left) and block diagrams (right)

The second subsystem of the port-Hamiltonian model for continuum mechanics characterizes the internal stress acting on the elastic body. We shall only present the convective representation of this port-Hamiltonian model depicted in Fig. 3. Its spatial and material counterparts are easily deducible.

In general, one has on part of the boundary, which we denote by ${\partial \mathcal {B}}_1$, the velocity is an input (i.e., boundary data), and traction is an output. On the other part of the boundary, which we denote by ${\partial \mathcal {B}}_2$, the traction is an input and the velocity in an output. One has that ${\partial \mathcal {B}}= {\partial \mathcal {B}}_1 \cup {\partial \mathcal {B}}_2$, ${\partial \mathcal {B}}_1 \cap {\partial \mathcal {B}}_2 = \emptyset$.

The pairing of the convective velocity ${\hat{v}}\in \varOmega ^{0}(\mathcal {B};T\mathcal {B})$ and stress $\hat{\mathcal {T}}\in \varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B})$ on the boundary of the body is defined as

$$\begin{aligned} \mathscr {P}_{st}:= \int _{\partial \mathcal {B}} i^*({\hat{v}}\ \dot{\wedge }\ \hat{\mathcal {T}}) = \int _{\partial \mathcal {B}} {\hat{v}}|_{\partial \mathcal {B}}\ \dot{\wedge }\ \hat{\mathcal {T}}|_{\partial \mathcal {B}} = \int _{\partial \mathcal {B}_1} {\hat{v}}|_{\partial \mathcal {B}_1}\ \dot{\wedge }\ \hat{\mathcal {T}}|_{\partial \mathcal {B}_1} + \int _{\partial \mathcal {B}_2} {\hat{v}}|_{\partial \mathcal {B}_2}\ \dot{\wedge }\ \hat{\mathcal {T}}|_{\partial \mathcal {B}_2}, \end{aligned}$$

where ${\hat{v}}|_{\partial \mathcal {B}}:= i_{\textrm{f}}^*({\hat{v}}) \in \varOmega ^{0}(\partial \mathcal {B};T\mathcal {B})$ and $\hat{\mathcal {T}}|_{\partial \mathcal {B}}:= i_{\textrm{f}}^*(\hat{\mathcal {T}}) \in \varOmega ^{n-1}(\partial \mathcal {B};T^*\mathcal {B})$ denote, respectively, the (partial) pullback of the convective velocity and stress on the boundary under the body inclusion map ${i}:{{\partial \mathcal {B}}}\rightarrow {\mathcal {B}}$. One has that ${\hat{v}}|_{\partial \mathcal {B}}$ represent the boundary’s velocity and $\hat{\mathcal {T}}|_{\partial \mathcal {B}}$ the traction on the boundary.

Furthermore, from the integration by parts formula for bundle-valued forms Rashad et al. (2023), one has that

$$\int _{\mathcal {B}}\hat{\nabla }{\hat{v}}\ \dot{\wedge }\ \hat{\mathcal {T}}+ {\hat{v}}\ \dot{\wedge }\ \hat{\textrm{d}}_{\hat{\nabla }}\hat{\mathcal {T}}= \int _{\mathcal {B}}\textrm{d}({\hat{v}}\ \dot{\wedge }\ \hat{\mathcal {T}}) = \int _{\partial \mathcal {B}} i^*({\hat{v}}\ \dot{\wedge }\ \hat{\mathcal {T}}).$$

Combining the above expressions together yields the stress power balance

$$\begin{aligned} \int _{\mathcal {B}}\hat{\nabla }{\hat{v}}\ \dot{\wedge }\ \hat{\mathcal {T}}+ {\hat{v}}\ \dot{\wedge }\ \hat{\textrm{d}}_{\hat{\nabla }}\hat{\mathcal {T}}= \int _{\partial \mathcal {B}_1} {\hat{v}}|_{\partial \mathcal {B}_1}\ \dot{\wedge }\ \hat{\mathcal {T}}|_{\partial \mathcal {B}_1} + \int _{\partial \mathcal {B}_2} {\hat{v}}|_{\partial \mathcal {B}_2}\ \dot{\wedge }\ \hat{\mathcal {T}}|_{\partial \mathcal {B}_2}. \end{aligned}$$

(56)

The above power balance (56) can be encompassed into a Stokes–Dirac³ structure denoted by ${\hat{\mathcal {D}}}_{\textrm{s}}\subset {\hat{\mathfrak {F}}}_{\textrm{s}}\times {\hat{\mathfrak {F}}^*}_{\textrm{s}}$, with ${\hat{\mathfrak {F}}}_{\textrm{s}}:= \varOmega ^{0}(\mathcal {B};T\mathcal {B})\times \varOmega ^{1}(\mathcal {B};T\mathcal {B})\times \varOmega ^{0}(\partial \mathcal {B};T\mathcal {B})$ and ${\hat{\mathfrak {F}}^*}_{\textrm{s}}$ denoting its dual. We define this (constant) Stokes–Dirac structure ${\hat{\mathcal {D}}}_{\textrm{s}}$ as

$$\begin{aligned} \begin{aligned} {\hat{\mathcal {D}}}_{\textrm{s}}:=&\left\{ (({\hat{f}}_{\textrm{d}}, {\hat{f}}_{\textrm{s}}, \hat{f}_\partial ),({\hat{e}}_{\textrm{d}}, {\hat{e}}_{\textrm{s}}, \hat{e}_\partial )) \in {\hat{\mathfrak {F}}}_{\textrm{s}}\times {\hat{\mathfrak {F}}^*}_{\textrm{s}}| \right. \\&\quad {\hat{f}}_{\textrm{s}}= \hat{\nabla }{\hat{f}}_{\textrm{d}},\qquad \qquad \ \ {\hat{e}}_{\textrm{d}}= \hat{\textrm{d}}_{\hat{\nabla }}{\hat{e}}_{\textrm{s}}\\&\quad {\hat{f}}_{\textrm{s}}|_{{\partial \mathcal {B}}_1}=\hat{f}_\partial ^1,\qquad \qquad {\hat{e}}_{\textrm{s}}|_{{\partial \mathcal {B}}_2}=\hat{e}_\partial ^2\\&\quad \left. \hat{f}_\partial ^2={\hat{f}}_{\textrm{s}}|_{{\partial \mathcal {B}}_2},\qquad \qquad \hat{e}_\partial ^1= {\hat{e}}_{\textrm{s}}|_{{\partial \mathcal {B}}_1} \right\} , \end{aligned} \end{aligned}$$

(57)

with $\hat{f}_\partial = (\hat{f}_\partial ^1,\hat{f}_\partial ^2) \in \varOmega ^{0}(\partial \mathcal {B}_{1};T\mathcal {B})\times \varOmega ^{0}(\partial \mathcal {B}_{2};T\mathcal {B})$ and $\hat{e}_\partial = (\hat{e}_\partial ^1,\hat{e}_\partial ^2) \in \varOmega ^{n-1}(\partial \mathcal {B}_{1};T^*\mathcal {B})\times \varOmega ^{n-1}(\partial \mathcal {B}_{2};T^*\mathcal {B})$. The Dirac structure (57) characterizes the power balance

$$\begin{aligned} \langle {\hat{e}}_{\textrm{d}} | {\hat{f}}_{\textrm{d}} \rangle _{\mathcal {B}}+ \langle {\hat{e}}_{\textrm{s}} | {\hat{f}}_{\textrm{s}} \rangle _{\mathcal {B}}= \langle \hat{e}_\partial ^1 | \hat{f}_\partial ^1 \rangle _{{\partial \mathcal {B}}_1} + \langle \hat{e}_\partial ^2 | \hat{f}_\partial ^2 \rangle _{{\partial \mathcal {B}}_2}, \end{aligned}$$

(58)

which is equivalent to (56) by setting the inputs $({\hat{f}}_{\textrm{d}},{\hat{e}}_{\textrm{s}},\hat{f}_\partial ^1,\hat{e}_\partial ^2) = ({\hat{v}},\hat{\mathcal {T}},{\hat{v}}|_{\partial \mathcal {B}_1},\hat{\mathcal {T}}|_{\partial \mathcal {B}_2})$ which yields the outputs $({\hat{e}}_{\textrm{d}}, {\hat{f}}_{\textrm{s}}, \hat{e}_\partial ^1, \hat{f}_\partial ^2)=(\hat{\textrm{d}}_{\hat{\nabla }}\hat{\mathcal {T}}, \hat{\nabla }{\hat{v}}, \hat{\mathcal {T}}|_{\partial \mathcal {B}_1}, {\hat{v}}|_{\partial \mathcal {B}_2})$.

Finally, it is straightforward to transform ${\hat{\mathcal {D}}}_{\textrm{s}}$ into a material one ${\tilde{\mathcal {D}}}_{\textrm{s}}$ (or spatial one ${{\mathcal {D}}}_{\textrm{s}}$) which will have identical structure to (57) with $\hat{\textrm{d}}_{\hat{\nabla }}$ replaced by its material counterpart $\tilde{\textrm{d}}_{\tilde{\nabla }}$ (or spatial counterpart $\textrm{d}_{\nabla }$).

3.4 Decomposition of Power Ports

So far we constructed the port-Hamiltonian subsystems that describe storage of kinetic energy and incorporate traction and velocity boundary conditions into the dynamics. The last step is to add the elastic body’s constitutive model to the power port $(\hat{\nabla }{\hat{v}},\hat{\mathcal {T}})$.

Fig. 4

Symmetric–asymmetric decomposition of $(\hat{\nabla }{\hat{v}},\hat{\mathcal {T}})$ and the volumetric–deviatoric decomposition of $(\hat{\varepsilon },{\hat{\mathcal {T}}}_{\textrm{sym}})$

Before doing so, we will discuss an essential decomposition of the power port $(\hat{\nabla }{\hat{v}},\hat{\mathcal {T}})$ that will split it into a symmetric part $(\hat{\varepsilon },{\hat{\mathcal {T}}}_{\textrm{sym}})$ and an asymmetric part $(\hat{w},{\hat{\mathcal {T}}}_{\textrm{asy}})$. By doing so, we will factor out the rigid body motion information included in the flow variable $\hat{\nabla }{\hat{v}}$ and the constitutive relation will be added to the symmetric port $(\hat{\varepsilon },{\hat{\mathcal {T}}}_{\textrm{sym}})$ only. Furthermore, we will also discuss another decomposition of the port $(\hat{\varepsilon },{\hat{\mathcal {T}}}_{\textrm{sym}})$ which will split it into a volumetric part $({\hat{\varepsilon }}_{\textrm{vol}},{\hat{\mathcal {T}}}_{\textrm{vol}})$ and a deviatoric part $({\hat{\varepsilon }}_{\textrm{dev}},{\hat{\mathcal {T}}}_{\textrm{dev}})$. Such decomposition is optional but very useful to distinguish the material’s response due to volumetric deformation from its response to isochoric processes (e.g., shear deformation). The two aforementioned decompositions will be implemented in a power consistent manner using Dirac structures, as illustrated in Fig. 4. The technical details of the aforementioned decompositions are discussed in Appendix A.3.

3.4.1 Symmetric–Asymmetric Decomposition

Based on the fact that $\varOmega ^{1}(\mathcal {B};T\mathcal {B}) \cong \varGamma (T^1_1\mathcal {B})$ and utilizing the Riemannian metric structure of $(\mathcal {B},\hat{g})$, one can define the projection map

$$\begin{aligned} {{\hat{\pi }}_{\textrm{sym}}}:{\varOmega ^{1}(\mathcal {B};T\mathcal {B})}\rightarrow {\varOmega ^{1}(\mathcal {B};T\mathcal {B})}, \qquad \qquad {\hat{\pi }}_{\textrm{sym}}:=\hat{g}^{-1}\circ \textrm{sym}\circ \hat{g}, \end{aligned}$$

(59)

where $\textrm{sym}$ denotes the symmetrization operation of 2-covariant tensor fields. Consequently, one can decompose the space of vector-valued one forms such that

$$\varOmega ^{1}(\mathcal {B};T\mathcal {B}) = \varOmega _{\textrm{sym}}^{1}(\mathcal {B};T\mathcal {B}) \oplus \varOmega _{\textrm{asy}}^{1}(\mathcal {B};T\mathcal {B}),$$

where $\varOmega _{\textrm{sym}}^{1}(\mathcal {B};T\mathcal {B}):= \textrm{im}({\hat{\pi }}_{\textrm{sym}})$ and $\varOmega _{\textrm{asy}}^{1}(\mathcal {B};T\mathcal {B}):= \textrm{ker}({\hat{\pi }}_{\textrm{sym}})$. Applying the above decomposition on the velocity gradient $\hat{\nabla }{\hat{v}}$ allows us to express it as

$$\hat{\nabla }{\hat{v}}= \hat{\varepsilon }+ \hat{w}, \qquad \qquad \hat{\varepsilon }:= {\hat{\pi }}_{\textrm{sym}}(\hat{\nabla }{\hat{v}}), \qquad \hat{w}:= {\hat{\pi }}_{\textrm{asy}}(\hat{\nabla }{\hat{v}}):= \hat{\nabla }{\hat{v}}- {\hat{\pi }}_{\textrm{sym}}(\hat{\nabla }{\hat{v}}),$$

where $\hat{\varepsilon }\in \varOmega _{\textrm{sym}}^{1}(\mathcal {B};T\mathcal {B})$ represents the convective rate of strain field while $\hat{w}\in \varOmega _{\textrm{asy}}^{1}(\mathcal {B};T\mathcal {B})$ represents the convective vorticity field⁴.

While the velocity field ${\hat{v}}$ includes information on how the elastic body deforms as well as the translation and rotation (i.e. rigid body motion), it performs in the ambient space; the deformation is described only by the symmetric component $\hat{\varepsilon }$ and not $\hat{w}$. This can be seen by inspecting the covariant velocity gradient $\hat{\nabla }\hat{v}^\flat := \hat{g}\cdot \hat{\nabla }{\hat{v}}\in \varOmega ^{1}(\mathcal {B};T^*\mathcal {B})$ which can be written as (Rashad et al. 2023, Prop. 1)

$$\hat{\nabla }\hat{v}^\flat = \textrm{sym}(\hat{\nabla }\hat{v}^\flat ) + \textrm{skew}(\hat{\nabla }\hat{v}^\flat ) = \frac{1}{2}\mathcal {L}_{{\hat{v}}}{\hat{g}} + \frac{1}{2}\textrm{d}\hat{v}^\flat ,.$$

For rigid body motion, one has that $\mathcal {L}_{{\hat{v}}}{\hat{g}} = \varphi ^*(\mathcal {L}_{v}{g}) = 0$, i.e., ${\hat{v}}$ and $v$ are Killing vector fields. Consequently, only the symmetric part of $\hat{\nabla }\hat{v}^\flat$ and $\nabla v^\flat$ describes the deformation of the body. Thus, the operation $\textrm{sym}\circ \hat{\nabla }$ factors out the rigid body motion information included in $\hat{v}^\flat$. The rate of strain tensor $\hat{\varepsilon }$ plays an important role in finite-strain elasticity theory as will be shown in the coming section.

By duality, one can decompose the stress tensor field $\hat{\mathcal {T}}\in \varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B})$ as the sum of ${\hat{\mathcal {T}}}_{\textrm{sym}} \in \varOmega _{\textrm{sym}}^{n-1}(\mathcal {B};T^*\mathcal {B})$ and ${\hat{\mathcal {T}}}_{\textrm{asy}} \in \varOmega _{\textrm{asy}}^{n-1}(\mathcal {B};T^*\mathcal {B})$. With reference to Fig. 4, the above symmetric–asymmetric decomposition of the power port $(\hat{\nabla }{\hat{v}},\hat{\mathcal {T}})$ can be encompassed in the Dirac structure $\hat{\mathcal {D}}_\text {s/a}(\hat{g})$, which is modulated by $\hat{g}$ and defined as the relation corresponding to the following map:

$$\begin{aligned} \begin{pmatrix} \hat{\varepsilon }\\ \hat{w}\\ \hat{\mathcal {T}} \end{pmatrix} = \begin{pmatrix} 0 & 0 & {\hat{\pi }}_{\textrm{sym}}\\ 0 & 0 & {\hat{\pi }}_{\textrm{asy}}\\ \mathbbm {1} & \mathbbm {1} & 0 \end{pmatrix} \begin{pmatrix} {\hat{\mathcal {T}}}_{\textrm{sym}}\\ {\hat{\mathcal {T}}}_{\textrm{asy}}\\ \hat{\nabla }{\hat{v}} \end{pmatrix}, \end{aligned}$$

(60)

which characterizes the power balance

$$\begin{aligned} \langle \hat{\mathcal {T}} | \hat{\nabla }{\hat{v}} \rangle _{\mathcal {B}} = \langle {\hat{\mathcal {T}}}_{\textrm{sym}} | \hat{\varepsilon } \rangle _{\mathcal {B}} + \langle {\hat{\mathcal {T}}}_{\textrm{asy}} | \hat{w} \rangle _{\mathcal {B}}. \end{aligned}$$

(61)

Remark 6

Note that if one assumes that $\hat{\mathcal {T}}\in \varOmega _{\textrm{sym}}^{n-1}(\mathcal {B};T^*\mathcal {B})$, then the power balance (61) simplifies to

$$\begin{aligned} \langle \hat{\mathcal {T}} | \hat{\nabla }{\hat{v}} \rangle _{\mathcal {B}} = \langle \hat{\mathcal {T}} | \hat{\varepsilon } \rangle _{\mathcal {B}}, \end{aligned}$$

(62)

which leads to the condition stated in Kanso et al. (2007); Gilbert and Vanneste (2023):

$$\begin{aligned} (\hat{\alpha } \otimes \hat{\beta }^\sharp )\ \dot{\wedge }\ \hat{\mathcal {T}}= (\hat{\beta } \otimes \hat{\alpha }^\sharp )\ \dot{\wedge }\ \hat{\mathcal {T}}, \qquad \forall \hat{\alpha },\hat{\beta } \in \varOmega ^{1}(\mathcal {B}). \end{aligned}$$

In the formulation above, we have shown that it follows naturally from duality as a consequence of factoring out rigid body motions. In numerical methods, assuming $\hat{\mathcal {T}}\in \varOmega _{\textrm{sym}}^{n-1}(\mathcal {B};T^*\mathcal {B})$ corresponds to the strong imposition of symmetry while using the formulation in (60) with the vorticity $\hat{w}$ as a Lagrange multiplier corresponds to the weak imposition of symmetry. The interested reader is referred to Arnold (2018).

Remark 7

(Spatial and Material counterparts of $\hat{\mathcal {D}}_{s/a}$) The above construction can be identically repeated for the spatial representation by defining the projection map ${\pi }_{\textrm{sym}}:= g^{-1}\circ \textrm{sym}\circ g$ to factor out rigid body motions. This is used in Sect. 4 when we discuss fluids. On the other hand, the above construction does not hold for the material representation. An immediate reason is the fact that the symmetrization operation of two-point tensors is not defined simply because its two indices belong to different spaces. This has a strong relation to the axiom of material frame independence of constitutive laws. The interested reader can refer to (Rashad et al. 2023, Sect. 6.5) and (Marsden and Hughes 1994, Ch.3).

3.4.2 Volumetric–Deviatoric Decomposition

Now we turn attention to the volumetric–deviatoric decomposition of the power port $(\hat{\varepsilon },\hat{\mathcal {T}})$ which can be constructed along the same line of thought as above and is detailed in Appendix A.3. For notational simplicity, we assume the symmetry $\hat{\mathcal {T}}$ to be enforced strongly (cf. Remark 6).

Let ${{\hat{\pi }}_{\textrm{vol}}}:{\varOmega ^{1}(\mathcal {B};T\mathcal {B})}\rightarrow {\varOmega _{\textrm{vol}}^{1}(\mathcal {B};T\mathcal {B})}$ denote the volumetric projection map defined by ${\hat{\pi }}_{\textrm{vol}}(\alpha ): = \frac{1}{n}\textrm{tr}(\alpha )I_n,$ for any $\alpha \in \varOmega ^{1}(\mathcal {B};T\mathcal {B})$, with ${\textrm{tr}}:{\varOmega ^{1}(\mathcal {B};T\mathcal {B})}\rightarrow {\varOmega ^{0}(\mathcal {B})}$ denoting the trace map. Furthermore, let ${{\hat{\pi }}_{\textrm{dev}}}:{\varOmega ^{1}(\mathcal {B};T\mathcal {B})}\rightarrow {\varOmega _{\textrm{dev}}^{1}(\mathcal {B};T\mathcal {B})}$ denote the deviatoric projection map defined by ${\hat{\pi }}_{\textrm{dev}}(\alpha ):= \alpha - {\hat{\pi }}_{\textrm{vol}}(\alpha )$. The volumetric–deviatoric decomposition is implemented by the Dirac structure $\hat{\mathcal {D}}_{v/d}$ defined as the relation corresponding to the following map:

$$\begin{aligned} \begin{pmatrix} {\hat{\varepsilon }}_{\textrm{vol}}\\ {\hat{\varepsilon }}_{\textrm{dev}}\\ \hat{\mathcal {T}} \end{pmatrix} = \begin{pmatrix} 0 & 0 & \textrm{tr}\circ {\hat{\pi }}_{\textrm{vol}}\\ 0 & 0 & {\hat{\pi }}_{\textrm{dev}}\\ \textrm{tr}^* & \mathbbm {1} & 0 \end{pmatrix} \begin{pmatrix} {\hat{\mathcal {T}}}_{\textrm{vol}}\\ {\hat{\mathcal {T}}}_{\textrm{dev}}\\ \hat{\varepsilon } \end{pmatrix}, \end{aligned}$$

(63)

where ${\hat{\varepsilon }}_{\textrm{vol}}\in \varOmega ^{0}(\mathcal {B})$ and ${\hat{\mathcal {T}}}_{\textrm{vol}} \in \varOmega ^{n}(\mathcal {B})$ are the volumetric components and ${\hat{\varepsilon }}_{\textrm{dev}}\in \varOmega _{\textrm{dev}}^{1}(\mathcal {B};T\mathcal {B})$ and ${\hat{\mathcal {T}}}_{\textrm{dev}} \in \varOmega _{\textrm{dev}}^{n-1}(\mathcal {B};T^*\mathcal {B})$ are the deviatoric components of the rate of strain and stress tensors, respectively, while ${\textrm{tr}^*}:{\varOmega ^{n}(\mathcal {B})}\rightarrow {\varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B})}$ denotes the dual trace map.

The volumetric component of the rate of strain ${\hat{\varepsilon }}_{\textrm{vol}}$ is an important quantity that describes the volumetric deformation of the elastic body. In fact, one can show that is equal to the divergence of the convective velocity:

$$\begin{aligned} \widehat{\textrm{div}}({\hat{v}}) = \textrm{tr}(\hat{\nabla }\hat{v}^\flat ) = \textrm{tr}(\hat{\varepsilon }) + \textrm{tr}(\hat{w}) = {\hat{\varepsilon }}_{\textrm{vol}}, \end{aligned}$$

(64)

which follows since the trace of $\hat{w}$ is equal to zero due to its skew-symmetric nature.

Finally, the power balance characterized by $\hat{\mathcal {D}}_{v/d}$ is given by

$$\begin{aligned} \langle \hat{\mathcal {T}} | \hat{\varepsilon } \rangle _{\mathcal {B}} = \langle {\hat{\mathcal {T}}}_{\textrm{vol}} | {\hat{\varepsilon }}_{\textrm{vol}} \rangle _{\mathcal {B}} + \langle {\hat{\mathcal {T}}}_{\textrm{dev}} | {\hat{\varepsilon }}_{\textrm{dev}} \rangle _{\mathcal {B}}, \end{aligned}$$

(65)

where the first pairing in the RHS is in terms of scalar-valued forms.

3.5 Intensive Rougee Stress Tensor Field

The last point to consider is the relation between the (extensive) stress variables we introduced so far as covector-valued forms in $\varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B})$ and their usual representation as (intensive) second-rank tensor fields. The identification between the two is achieved using the metric structure represented by the Hodge star $\hat{\star }_c$. Let $\hat{\tau }:= \hat{\star }_c^{-1}\hat{\mathcal {T}}\in \varOmega ^{1}(\mathcal {B};T\mathcal {B}) \cong \varGamma (T^1_1\mathcal {B})$ denote the intensive version of $\hat{\mathcal {T}}\in \varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B})$. In local coordinates for $n=3$, one can see that Rashad et al. (2023)

$$\begin{aligned} \hat{\mathcal {T}}_{AIJ} = \hat{\tau }_A^B \hat{\mu }_{BIJ} \in \varOmega ^{n-1}(\mathcal {B}), \end{aligned}$$

(66)

which indicates that $\hat{\mathcal {T}}$ contains two types of information: one related to the intrinsic mass form $\hat{\mu }$ and another related to the intensive (i.e., mass and volume independent) stress which should be provided by the elastic properties of the body. Indeed, this combination in (66) is what makes the form part and value part explicit in the exterior calculus formulation. However, for the purposes of defining the stress–strain constitutive relation, it suffices to provide an expression relating $\hat{\tau }$ and $\hat{\varepsilon }$.

The condition that $\hat{\mathcal {T}}\in \varOmega _{\textrm{sym}}^{n-1}(\mathcal {B};T^*\mathcal {B})$ is equivalent to the symmetrization of the 2-contravariant tensor $\hat{\tau }^\sharp \in \varGamma (T^2_0\mathcal {B})$, i.e., $\hat{\tau }^{IJ} = \hat{\tau }^{JI}$ which is the standard stress-symmetry requirement in tensor calculus. Furthermore, the expression of $\hat{\mathcal {T}}$ in (63):

$$\hat{\mathcal {T}}= \textrm{tr}^*({\hat{\mathcal {T}}}_{\textrm{vol}}) + {\hat{\mathcal {T}}}_{\textrm{dev}}$$

is equivalent to

$$\begin{aligned} \hat{\tau }= {\hat{\tau }}_{\textrm{vol}} I_n + {\hat{\tau }}_{\textrm{dev}}, \end{aligned}$$

(67)

with ${\hat{\tau }}_{\textrm{vol}} = \star ^{-1} {\hat{\mathcal {T}}}_{\textrm{vol}} \in \varOmega ^{0}(\mathcal {B})$ and ${\hat{\tau }}_{\textrm{dev}} = \hat{\star }_c^{-1}{\hat{\mathcal {T}}}_{\textrm{dev}} \in \varOmega ^{1}(\mathcal {B};T\mathcal {B})$ denoting the volumetric and deviatoric components of $\hat{\tau }$, respectively. The power balance (65) can then be rewritten in exterior and tensor calculus, respectively, as

$$\begin{aligned} \int _{\mathcal {B}}\hat{\varepsilon }\ \dot{\wedge }\ \hat{\star }_c\hat{\tau }= \int _{\mathcal {B}}{\hat{\varepsilon }}_{\textrm{vol}}\wedge \star {\hat{\tau }}_{\textrm{vol}}+ {\hat{\varepsilon }}_{\textrm{dev}}\ \dot{\wedge }\ \hat{\star }_c{\hat{\tau }}_{\textrm{dev}}= \int _{\mathcal {B}}({\hat{\varepsilon }}_{\textrm{vol}}{\hat{\tau }}_{\textrm{vol}}+ {\hat{\varepsilon }}_{\textrm{dev}}: {\hat{\tau }}_{\textrm{dev}}) \hat{\mu }. \end{aligned}$$

(68)

Note that in the above power balance the inner product of ${\hat{\varepsilon }}_{\textrm{vol}}$ and ${\hat{\tau }}_{\textrm{vol}}$ is in terms of scalar-valued forms with the Hodge star $\hat{\star }$ including mass information such that $\hat{\star } \hat{\mu }= \mathbbm {1}$ where $\mathbbm {1}$ denotes the identity function on $\mathcal {B}$. In Sect. 5, we will discuss a number of constitutive models to close the relation between $\hat{\varepsilon }$ and $\hat{\tau }$.

3.6 Summary of the Convective port-Hamiltonian Dynamics

Fig. 5

Convective port-Hamiltonian structure of nonlinear elasticity showing its constituting subsystems in bond graph notation. The boundary ports are depicted without causality information for simplicity

In conclusion, the port-Hamiltonian dynamics we constructed so far is characterized by the storage of kinetic energy, one port for traction boundary input, one port for velocity boundary input, and the open distributed port $(\hat{\varepsilon },\hat{\star }_c\hat{\tau })$ which will be closed by the stress constitutive law. At the heart of the model are the three Dirac structures ${\hat{\mathcal {D}}}_{\textrm{k}}, {\hat{\mathcal {D}}}_{\textrm{s}}, \hat{\mathcal {D}}_{s/a}$ which we combine into ${\hat{\mathcal {D}}}_{\textrm{solid}}$ and will characterize the combined power balances (50,58,62). With reference to Fig. 5, the resulting port-Hamiltonian model in the convective representation is summarized for convenience by the following theorem.

Theorem 3

The equations of motion governing the convective metric and extensive momentum variables $\hat{\chi }:= (\hat{g},\hat{\mathcal {M}}) \in \mathcal {M}(\mathcal {B})\times \varOmega ^{n}(\mathcal {B};T^*\mathcal {B}) =: {\hat{\mathcal {X}}}_{\textrm{k}}$ are represented in the port-Hamiltonian framework as

$$\begin{aligned} \begin{pmatrix} \partial _t \hat{g}\\ \partial _t \hat{\mathcal {M}} \end{pmatrix}&= \begin{pmatrix} 0 & 2\ \textrm{sym}\circ \hat{g}\circ \hat{\textrm{d}}_{\hat{\nabla }} \\ 2\ \hat{\textrm{d}}_{\hat{\nabla }}\circ \hat{g} & \mathcal {L}_{(\cdot )}{\hat{\mathcal {M}}} \end{pmatrix} \begin{pmatrix} \delta _{\hat{g}}{{\hat{H}}_{\textrm{k}}}\\ \delta _{\hat{\mathcal {M}}}{{\hat{H}}_{\textrm{k}}} \end{pmatrix} + \begin{pmatrix} 0\\ \hat{\textrm{d}}_{\hat{\nabla }} \end{pmatrix} \hat{\star }_c\hat{\tau }, \end{aligned}$$

(69)

$$\begin{aligned} \hat{\varepsilon }&= \left( 0 \quad {\hat{\pi }}_{\textrm{sym}} \circ \hat{\nabla }\right) \begin{pmatrix} \delta _{\hat{g}}{{\hat{H}}_{\textrm{k}}}\\ \delta _{\hat{\mathcal {M}}}{{\hat{H}}_{\textrm{k}}} \end{pmatrix} \end{aligned}$$

(70)

where the Hamiltonian functional ${{\hat{H}}_{\textrm{k}}}:{{\hat{\mathcal {X}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ is given by (47), the variational derivatives of ${\hat{H}}_{\textrm{k}}$ are given by Prop. 5, and the projection map ${\hat{\pi }}_{\textrm{sym}}$ is given by (59).. The symmetry of the stress $\hat{\tau }\in \varOmega _{\textrm{sym}}^{1}(\mathcal {B};T\mathcal {B})$ is assumed, while the boundary data are given by ${\hat{v}}|_{\partial \mathcal {B}}$ and $\hat{\mathcal {T}}|_{\partial \mathcal {B}}$, with $\hat{\mathcal {T}}= \hat{\star }_c\hat{\tau }\in \varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B})$.

The rate of change of ${\hat{H}}_{\textrm{k}}$ along trajectories of (69-70) is given by

$$\begin{aligned} {\dot{\hat{H}}}_{\textrm{k}}= \int _{\partial \mathcal {B}} {\hat{v}}|_{\partial \mathcal {B}}\ \dot{\wedge }\ \hat{\mathcal {T}}|_{\partial \mathcal {B}} - \int _{\mathcal {B}}\hat{\varepsilon }\ \dot{\wedge }\ \hat{\star }_c\hat{\tau }. \end{aligned}$$

(71)

We finally conclude by a number of insights that the port-Hamiltonian model above reveals, in contrast with the tensor calculus PDEs in Sect. 2.1.

• 1. The formulation of nonlinear elasticity above explicates i) the geometric structure using coordinate-free expressions, ii) the topological structures using exterior calculus, and iii) the energetic structure using the Dirac structures and power ports.

• 2. There is a clear role of duality between effort and flow variables that make up the power ports.

• 3. There is a clear identification of the skew-symmetric operator ${\hat{J}(\hat{\chi })}:{T^* {\hat{\mathcal {X}}}_{\textrm{k}}}\rightarrow {T {\hat{\mathcal {X}}}_{\textrm{k}}}$ which represents the Poisson structure, associated to the convective Poisson bracket $^c\{\cdot ,\cdot \}$ which can be defined similar to (36). The skew symmetry of $\hat{J}(\hat{\chi })$ is directly related to the conservation of energy characterized by the Poisson bracket.

• 4. There is a clear separation between the constitutive relations and the Dirac structures. One key consequence is that one distinguishes between the nonlinearities of each, for example as in the state dependency of ${\hat{\mathcal {D}}}_{\textrm{k}}(\hat{\chi })$.

• 5. It shows that the evolution of the dynamic equation (in the convective and spatial case) is independent of the configuration ${\varphi }:{\mathcal {B}}\rightarrow {\mathcal {S}}$. In fact, one case reconstructs it by solving the ordinary differential equation

  $$\begin{aligned} \partial _t \varphi ^{-1}(x,t) = - {\hat{v}}(\varphi ^{-1}(x,t),t). \end{aligned}$$

  (72)

• 6. The boundary conditions and how they affect the power balance can be easily identified. In fact, it shows that one does not need to impose displacement boundary conditions but rather velocity boundary conditions. This is also a crucial point when interconnecting flexible and rigid bodies that leads either to displacement constraints or velocity constraints. Indeed, both types of boundary conditions are one to one as shown in (72). However, it has several consequences on numerical methods. We shall come back to this point in Sect. 5.4 when discussing incompressibility.

• 7. The equivalence between the balance laws corresponding to the port-Hamiltonian model (69-70), cf. Corollary 3, and their tensor calculus counterparts (9-12) can be found in Rashad et al. (2023).

4 Port-Hamiltonian Model of Fluid Flow

Nonlinear elasticity has a lot of similarity with fluid mechanics. A compressible fluid can be considered a special case of visco-elastic materials (Marsden and Hughes 1994). Even though the governing equations of motion and the Hamiltonian reduction process of both dynamical systems are identical, there a few fundamental differences between fluids and solids from the geometric point of view. In this small section, we shall highlight these similarities and differences by presenting a port-Hamiltonian formulation of the dynamics of fluid flow.

4.1 Geometric Formulation

In principle, the geometric setting of fluids can also be described by a matter space $\mathcal {B}$ that is embedded in an ambient space $\mathscr {A}$ along the same line of Sect. 2.1. In this setting, one should consider an additional control volume subdomain of $\mathscr {A}$ where the dynamics will be described as in Califano et al. (2022). However, it is more common in the literature to premise the geometric formulation of fluids by assuming the fluid particles to form a continuum filling the whole spatial domain, denoted by $\mathcal {S}$, which we assume to have a Riemannian structure (Arnold 1965; Marsden and Weinstein 1982). The configuration of the fluid is given by the diffeomorphism ${\varphi }:{\mathcal {S}}\rightarrow {\mathcal {S}}$ with the configuration space given by $\mathscr {C}:= \text {Diff}(\mathcal {S})$, the space of diffeomorphisms on $\mathcal {S}$. The configuration space has the structure of a Frechet–Lie group (Modin et al. 2011).

Similar to solids, the material velocity of the fluid is represented as a vector on the infinite-dimensional manifold $\mathscr {C}$ and one can define the spatial and convective velocity fields $v, {\hat{v}}\in \varGamma (T\mathcal {S})$ similar to (1). The Lie group structure of $\mathscr {C}$ allows us to interpret (1) as a pushforward of the material velocity $\tilde{v}\in T_\varphi \mathscr {C}$ to the left (convective) and right (spatial) Lie algebras, respectively.

The group structure of $\mathscr {C}$ and the diffeomorphic property of $\varphi$ imply physically that the fluid particles remain within $\mathcal {S}$ for all time. This is only the case if $\mathcal {S}$ is a compact manifold without a boundary or that it has a boundary ${\partial \mathcal {S}}$ that is impermeable, i.e., the fluid velocity is always parallel to the boundary (Rashad et al. 2021a). This boundary impermeability condition is expressed in exterior calculus as $\iota _{v}\omega _{g}|_{{\partial \mathcal {S}}} = 0 \in \varOmega ^{n-1}({\partial \mathcal {S}})$, where $\omega _{g}\in \varOmega ^{n}(\mathcal {S})$ denotes the volume form of the spatial domain associated to the metric g. In vector calculus, this condition corresponds to the vanishing of the normal component of the velocity at the boundary. All related works on the Hamiltonian formulation of fluids (e.g., Marsden and Weinstein (1983); Marsden et al. (1984); Holm et al. (1986)) are limited to this case of impermeable or no boundary. However, we shall show next how this assumption can be relaxed using Dirac structures during the Hamiltonian reduction procedure.

4.2 Port-Hamiltonian Model

The port-Hamiltonian modeling process proceeds in the same manner of Sect. 3:

• 1. The material equations of motion are given by the canonical Poisson structure (24) of the cotangent bundle $T^*\mathscr {C}$ such that the kinetic energy functional (18) is conserved.

• 2. Following the Hamiltonian reduction procedure in Theorem 1, one can transform the canonical Poisson structure to a Lie–Poisson structure in the spatial representation.

• 3. By extending the Poisson structure to a Dirac structure, one reaches the port-Hamiltonian dynamical system (43-44) representing the kinetic energy subsystem.

• 4. Then one adds the stress and velocity boundary conditions via the stress subsystem, as in Sect. 3.2, and factorizes the rigid body motions, as in Sect. 3.3.

• 5. Finally, one adds the constitutive relation between the rate of strain and stress variables $\varepsilon ,\tau \in \varOmega ^{1}(\mathcal {S};T\mathcal {S})$ based on the class of fluids to be modeled.

One prominent difference between fluids and solids is that the former are allowed to have momentum flux through the spatial image ${\partial \mathcal {S}}$ of the boundary, i.e., $\iota _v\mathcal {M}|_{{\partial \mathcal {S}}} \ne 0$. In contrast with the Hamiltonian formalism, one can account for such energy exchange through the boundary in Step 3 by extending the Dirac structure (37) with the boundary interaction port $(f_\partial ,e_\partial )\in \varOmega ^{0}(\partial \mathcal {S};T\mathcal {S})\times \varOmega ^{n-1}(\partial \mathcal {S};T^*\mathcal {S})$ with $f_\partial := {e}_{\mathcal {M}}|_{{\partial \mathcal {S}}}$ and $e_\partial := -({e}_{\mu }\mu + \iota _{{e}_{\mathcal {M}}}\mathcal {M})|_{{\partial \mathcal {S}}}$. Consequently, the power balance (38) then becomes

$$\begin{aligned} \langle {e}_{\mu } | {f}_{\mu } \rangle _{\mathcal {S}} + \langle {e}_{\mathcal {M}} | {f}_{\mathcal {M}} \rangle _{\mathcal {S}} = \langle {e}_{\textrm{d}} | {f}_{\textrm{d}} \rangle _{\mathcal {S}} + \langle e_\partial | f_\partial \rangle _{{\partial \mathcal {S}}}, \end{aligned}$$

(73)

and as a result the energy exchange due to momentum flux via ${\partial \mathcal {S}}$ appears in (45).

Interesting, but not surprisingly, the vanishing condition of $\iota _v\mathcal {M}|_{{\partial \mathcal {S}}}$ in Prop. 3 is equivalent to the boundary impermeability condition since $\iota _v\mathcal {M}|_{{\partial \mathcal {S}}} = \iota _{v}\mu \otimes v^\flat$. This in fact shows that the group structure of $\mathscr {C}$ and the diffeomorphic property of $\varphi$ are pertinent to the Hamiltonian formulation which necessitates such zero exchange of energy through the boundary (Rashad et al. 2021a).

Fig. 6

Spatial port-Hamiltonian structure of fluid flow showing its constituting subsystems in bond graph notation. The boundary ports are depicted without causality information for simplicity

In conclusion, the port-Hamiltonian representation of fluid dynamics, summarized in Fig. 6, is characterized by the storage of kinetic energy, boundary ports for traction, momentum flux, and velocity on ${\partial \mathcal {S}}$, and the open distributed port $(\varepsilon ,\star _c\tau )$ which will be closed by the stress constitutive law. At the heart of the model are the three Dirac structures ${{\mathcal {D}}}_{\textrm{k}}, {{\mathcal {D}}}_{\textrm{s}}, {\mathcal {D}}_{s/a}$ which we combine into ${{\mathcal {D}}}_{\textrm{fluid}}$. The resulting port-Hamiltonian model in the spatial representation is summarized by the following theorem.

Theorem 4

The equations of motion governing the extensive mass and momentum variables $\chi := (\mu ,\mathcal {M}) \in \varOmega ^{n}(\mathcal {S})\times \varOmega ^{n}(\mathcal {S};T^*\mathcal {S}) =: {{\mathcal {X}}}_{\textrm{k}}$ is represented in the port-Hamiltonian framework as

$$\begin{aligned} \begin{pmatrix} \partial _t \mu \\ \partial _t \mathcal {M} \end{pmatrix}&= \begin{pmatrix} 0 & - \textrm{d}\iota _{(\cdot )} \mu \\ - \mu \otimes \textrm{d}(\cdot ) & - \mathcal {L}_{(\cdot )}{\mathcal {M}} \end{pmatrix} \begin{pmatrix} \delta _{\mu }{{{H}}_{\textrm{k}}}\\ \delta _{\mathcal {M}}{{{H}}_{\textrm{k}}} \end{pmatrix} + \begin{pmatrix} 0\\ \textrm{d}_{\nabla } \end{pmatrix} {\star _c\tau }, \end{aligned}$$

(74)

$$\begin{aligned} \varepsilon&= \left( 0 \quad {\pi }_{\textrm{sym}}\circ \nabla \right) \begin{pmatrix} \delta _{\mu }{{{H}}_{\textrm{k}}}\\ \delta _{\mathcal {M}}{{{H}}_{\textrm{k}}} \end{pmatrix}, \end{aligned}$$

(75)

where the Hamiltonian functional ${{{H}}_{\textrm{k}}}:{{{\mathcal {X}}}_{\textrm{k}}}\rightarrow {\mathbb {R}}$ is given by (32), the variational derivatives of ${{H}}_{\textrm{k}}$ are given by Prop. 3, and the projection map ${{\pi }_{\textrm{sym}}}:{\varOmega ^{1}(\mathcal {S};T\mathcal {S})}\rightarrow {\varOmega ^{1}(\mathcal {S};T\mathcal {S})}$ is given by ${\pi }_{\textrm{sym}}:= g^{-1}\circ \textrm{sym}\circ g$. The symmetry of the stress $\tau \in \varOmega _{\textrm{sym}}^{1}(\mathcal {S};T\mathcal {S})$ is assumed, while the boundary data are given by $v|_{{\partial \mathcal {S}}}$ and $(\mathcal {T}-\frac{1}{2}\iota _{v}\mathcal {M})|_{{\partial \mathcal {S}}}$, with $\hat{\mathcal {T}}= \hat{\star }_c\hat{\tau }\in \varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B})$ and $v= \star _c^{-1}\mathcal {M}\in \varOmega ^{0}(\mathcal {S};T\mathcal {S})$.

The rate of change of ${{H}}_{\textrm{k}}$ along trajectories of (74-75) is given by

$${\dot{{H}}}_{\textrm{k}}= \int _{{\partial \mathcal {S}}} v|_{{\partial \mathcal {S}}} \ \dot{\wedge }\ (\mathcal {T}- \frac{1}{2}\iota _v\mathcal {M})|_{{\partial \mathcal {S}}} - \int _{\mathcal {S}}\varepsilon \ \dot{\wedge }\ \star _c\tau .$$

Finally, the equivalence between the balance laws corresponding to the port-Hamiltonian model (74-75), cf. Corollary 2, and their tensor calculus counterparts (6-8) can be found in Rashad et al. (2023).

5 Constitutive Relations

Now we turn attention to the constitutive relation between the rate of strain and stress variables, which are both symmetric tensor fields. We aim in this section to highlight how different classes of constitutive relations are incorporated in the port-Hamiltonian framework in addition to the importance of the convective representation for constitutive modeling. In particular, we aim to emphasize the role of $\mathcal {M}(\mathcal {B})$ in defining constitutive equations that are geometrically consistent, attributed to the work of Rougée (2006).

In general, neglecting thermodynamic effects, the constitutive law is abstractly a relation between the stress and rate of strain variables $\hat{\tau }, \hat{\varepsilon }\in \varOmega _{\textrm{sym}}^{1}(\mathcal {B};T\mathcal {B})$ and some time integral variable ${\hat{\zeta }}\in \mathfrak {D}$ of $\hat{\varepsilon }$ that characterizes deformation. Different choices for the construction of ${\hat{\zeta }}$ and its corresponding space $\mathfrak {D}$ will be discussed in detail in this section. Using $\hat{\varepsilon }$, one can rewrite the evolution equation for the convective metric (11) as

$$\begin{aligned} \partial _t \hat{g}= 2 \hat{g}\cdot \hat{\varepsilon }\in T_{\hat{g}}\mathcal {M}(\mathcal {B})\cong \varOmega ^{1}(\mathcal {B};T^*\mathcal {B}), \end{aligned}$$

(76)

which identifies vectors on $\mathcal {M}(\mathcal {B})$ with the covariant version of $\hat{\varepsilon }$. Therefore, it is expected that ${\hat{\zeta }}$ to be constructed from the metric state $\hat{g}$. The constitutive law will be denoted by

$$\mathscr {R}(\hat{\tau }, \hat{\varepsilon }, {\hat{\zeta }}) = 0.$$

In general, $\mathscr {R}$ depends also on $X\in \mathcal {B}$ to incorporate non-homogeneous effects in the elastic material. Without loss of generality, we shall only focus on the homogeneous case in this work.

From an energetic perspective, the power port $(\hat{\varepsilon },\hat{\star }_c\hat{\tau })$ has to be closed by either 1) an energy storage unit characterized by a relation between $\hat{\tau }$ and ${\hat{\zeta }}$, 2) an energy dissipation unit characterized by a relation between $\hat{\tau }$ and $\hat{\varepsilon }$, or 3) a combination of both. The three cases above classify elastic, viscous, and visco-elastic materials, respectively. The relations could be either linear or nonlinear, as well as static or dynamic, as in hyper-elasticity or hypo-elasticity, respectively.

In what follows, we shall only focus on the case of (homogeneous) hyper-elastic materials in Secs. 5.1-5.4 and viscous fluid flow in Sect. 5.5.

5.1 General Hyper-elasticity

The power port $(\hat{\varepsilon },\hat{\star }_c\hat{\tau })$ is closed for hyper-elastic materials by an energy storage unit characterizing the functional

$$\hat{\varPsi }[{\hat{\zeta }}] = \int _{\mathcal {B}}\hat{\psi }({\hat{\zeta }}) \hat{\mu },$$

where ${\hat{\psi }}:{\mathfrak {D}}\rightarrow {C^\infty (\mathcal {B})}$ is its specific storage energy function (which is independent of mass and volume).

Proposition 6

Let the variable ${\hat{\zeta }}\in \mathfrak {D}$ be constructed from the metric state $\hat{g}$ by some function ${f}:{\mathcal {M}(\mathcal {B})}\rightarrow {\mathfrak {D}}$ and let ${Tf}:{T_{\hat{g}} \mathcal {M}(\mathcal {B})}\rightarrow {T_{{\hat{\zeta }}}\mathfrak {D}}$ denote its tangent map. Then, the stress constitutive law takes the following generalized form of the Doyle–Ericskon formula

$$\begin{aligned} \hat{\tau }= 2 \hat{\star }_c^{-1}\hat{f}^*\left( \hat{\star }_c\frac{\partial \hat{\psi }}{\partial {\hat{\zeta }}}({\hat{\zeta }})\right) , \end{aligned}$$

(77)

where ${\hat{f}^*}:{T^*_{{\hat{\zeta }}}\mathfrak {D}}\rightarrow {\varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B})}$ denotes the dual of ${\hat{f}:= Tf \circ \hat{g}}:{\varOmega ^{1}(\mathcal {B};T\mathcal {B})}\rightarrow {T_{{\hat{\zeta }}}\mathfrak {D}}$ with respect to the duality pairing $\langle \cdot | \cdot \rangle _{\mathcal {B}}$.

Proof

Let $\delta _{\hat{\zeta }}\hat{\varPsi } \in T^*_{{\hat{\zeta }}}\mathfrak {D}$ denote the variational derivative of $\hat{\varPsi }$ with respect to ${\hat{\zeta }}$. The rate of change of $\hat{\varPsi }$ is then given by

$$\dot{\hat{\varPsi }} = \int _{\mathcal {B}}\partial _t {\hat{\zeta }}\ \dot{\wedge }\ \delta _{\hat{\zeta }}\hat{\varPsi } = \int _{\mathcal {B}}Tf(\partial _t \hat{g})\ \dot{\wedge }\ \delta _{\hat{\zeta }}\hat{\varPsi } = \int _{\mathcal {B}}2 Tf(\hat{g}\hat{\varepsilon })\ \dot{\wedge }\ \delta _{\hat{\zeta }}\hat{\varPsi } = \int _{\mathcal {B}}2 \hat{\varepsilon }\ \dot{\wedge }\ \hat{f}^*(\delta _{\hat{\zeta }}\hat{\varPsi } ).$$

A necessary condition for the energy balance $\dot{\hat{\varPsi }} = \int _{\mathcal {B}}\hat{\varepsilon }\ \dot{\wedge }\ \hat{\star }_c\hat{\tau },$ to hold and be covariant⁵ is that $\hat{\varPsi }$ should be independent of derivatives of ${\hat{\zeta }}$ (Marsden and Hughes 1994, Ch.3). Consequently, the variational derivative is identified with the partial derivative of the top form $\hat{\psi }({\hat{\zeta }}) \hat{\mu }$ leading to the expression $\delta _{\hat{\zeta }}\hat{\varPsi } = \hat{\star }_c\frac{\partial \hat{\psi }}{\partial {\hat{\zeta }}}$. $\square$

Using the definition of $\dot{\hat{\varPsi }} = \langle \delta _{\hat{\zeta }}\hat{\varPsi } | \partial _t {\hat{\zeta }} \rangle _{\mathcal {B}}$, the total energy balance of the port-Hamiltonian system of nonlinear elasticity (69-70) after adding the constitutive relation (77) can be expressed as

$$\frac{\textrm{d}}{\textrm{d}t} ({\hat{H}}_{\textrm{k}}+ \hat{\varPsi }) = \int _{{\partial \mathcal {B}}} {\hat{v}}|_{{\partial \mathcal {B}}} \ \dot{\wedge }\ \hat{\mathcal {T}}|_{{\partial \mathcal {B}}}.$$

The above power balance states that the rate of change of total (kinetic plus strain) energy within the elastic body’s domain $\mathcal {B}$ is equal to the power flow through the boundary ${\partial \mathcal {B}}$, due to stress.

The most natural choice for the quantity ${\hat{\zeta }}$ is $\hat{g}\in \mathcal {M}(\mathcal {B})$ itself since the space of Riemannian metrics serves as a natural space of deformation (Rashad et al. 2023; Stramigioli 2024). In this case, ${\hat{\zeta }}$ represents a state of deformation and one has $f = \text {id}$ and $\hat{f} = \hat{f}^* = \hat{g}$ which simplifies the stress law (77) to be

$$\begin{aligned} \hat{\tau }= 2 \hat{g}\frac{\partial \hat{\psi }}{\partial \hat{g}}(\hat{g}), \end{aligned}$$

(78)

which is the standard Doyle–Erickson formula.

The alternative choice which we shall adopt is to define ${\hat{\zeta }}$ as a strain variable represented by a vector-valued 1-form, i.e., $\mathfrak {D} = \varOmega ^{1}(\mathcal {B};T\mathcal {B})$. In the vast literature of nonlinear hyper-elasticity, there have been several propositions for defining the strain ${\hat{\zeta }}$ and several propositions for defining the storage energy $\hat{\psi }$ leading to numerous expressions for the stress $\hat{\tau }$ which very quickly ramifies when we consider their spatial and material counterparts along with tensorial variants. Thanks to the work of Rougée (2006) and Fiala (2011), taking the (non-Euclidean) geometric nature of $\mathcal {M}(\mathcal {B})$ into account leads to the fact that there is only one natural way to define the strain as the relative deformation between two states in $\mathcal {M}(\mathcal {B})$ as discussed next.

5.2 Finite-Strain Hyper-elasticity

Fig. 7

Deformation process in finite-strain theory

At this point, it is essential to introduce the concept of the reference metric $G\in \mathcal {M}(\mathcal {B})$ which will characterize the equilibrium reference state of deformation of the body. This equilibrium state also corresponds to the minimum of the storage elastic energy $\hat{\psi }$ considered as a function on $\mathcal {M}(\mathcal {B})$. Given G, one can associate with it the reference configuration ${{\varphi }_{\textrm{ref}}}:{\mathcal {B}}\rightarrow {\mathscr {A}}$ such that $G = {\varphi }_{\textrm{ref}}^*(g)$. With reference to Fig. 7, the deformation process will be characterized by the curve $c_{\hat{g}}(t)$ in $\mathcal {M}(\mathcal {B})$ starting at the initial state $\hat{g}_0:= c_{\hat{g}}(0) \in \mathcal {M}(\mathcal {B})$. If the deformation process starts from the reference configuration, then ${\varphi }_{\textrm{ref}}= \varphi _0$ and $\hat{g}_0 = G$, but in general the elastic body could be in a stressed state at $t=0$.

The strain variable ${\hat{\zeta }}_t$ at any time t is defined as the geodesic connecting $\hat{g}_t$ to the reference state G. It has been shown in Fiala (2011) that this geodesic is characterized by $\partial _t \hat{\varepsilon }= 0$, i.e., $\hat{\varepsilon }$ is constant along the geodesic, which leads to the following exponential and logarithm maps of $\mathcal {M}(\mathcal {B})$ (Fiala 2011)

$$\begin{aligned} \begin{aligned} \text {Exp}_G : {T_{G}\mathcal {M}(\mathcal {B})}&\rightarrow {\mathcal {M}(\mathcal {B})}\\ {\delta G}&\mapsto {G \textrm{exp}(G^{-1} \delta G)} \end{aligned} \qquad \qquad \begin{aligned} \text {Log}_G : {\mathcal {M}(\mathcal {B})}&\rightarrow {T_G\mathcal {M}(\mathcal {B})}\\ {\hat{g}}&\mapsto { G\textrm{ln}(G^{-1}\hat{g})} \end{aligned} \end{aligned}$$

(79)

where $\textrm{exp}$ and $\textrm{ln}$ denote the matrix exponential and natural logarithm maps, respectively, applied to mixed second-rank tensors. By introducing $\hat{C}_t:= G^{-1} \hat{g}_t$ as shorthand notation, we define the strain variable to the vector-valued 1-form given by

$$\begin{aligned} {\hat{\zeta }}_t:= \frac{1}{2}\ln (\hat{C}_t) \in \varOmega ^{1}(\mathcal {B};T\mathcal {B}), \end{aligned}$$

(80)

which leads to the nonlinear map $f(\hat{g}): = \frac{1}{2}\ln (G^{-1}\hat{g})$.

Remark 8

It is important to note that the expressions of (79) and (80) are mere notation for expressing the geodesic flow on $\mathcal {M}(\mathcal {B})$. Thus, one should not interpret $G^{-1} \hat{g}$ as an index raising of $\hat{g}$ by G since both metrics are two different points on the manifold $\mathcal {M}(\mathcal {B})$. In Fiala (2011), this distinction between thee geometric objects and their algebraic representation was made by denoting the latter with bold.

From the definition of the strain (80), one has that

$${\hat{\zeta }}_t = 0 \implies \hat{C}_t = I \implies \hat{g}_t = G,$$

stating that G is the zero-stress undeformed state of the body. However, it should warned that $\mathcal {M}(\mathcal {B})$ does not have this privileged state intrinsically. In fact, it has been shown by Fiala (2011) that this is a consequence of the geodesic completeness of $\mathcal {M}(\mathcal {B})$, i.e., any two deformation states can be connected through a geodesic. From the construction, we presented here, it can be seen that G is actually a property of the elastic body’s material. Some composite materials might even have multiple zero-stress states, such as bielastic structures used in soft-robotics.

As reviewed by Neff et al. (2015), the logarithmic strain measure (80) possesses a number of remarkable properties. Perhaps the most important one is that it additively separates dilation from pure distortion, as shown in the following important result.

Theorem 5

Consider the volumetric–deviatoric decomposition of the strain tensor

$$\begin{aligned} {\hat{\zeta }}= \frac{1}{n} {{\hat{\zeta }}}_{\textrm{vol}}I_n + {{\hat{\zeta }}}_{\textrm{dev}}, \end{aligned}$$

(81)

where ${{\hat{\zeta }}}_{\textrm{vol}}:= \textrm{tr}({\hat{\zeta }}) \in \varOmega ^{0}(\mathcal {B})$ and ${{\hat{\zeta }}}_{\textrm{dev}}:= {\hat{\zeta }}- \frac{1}{n} {{\hat{\zeta }}}_{\textrm{vol}}I_n \in \varOmega _{\textrm{dev}}^{1}(\mathcal {B};T\mathcal {B})$. Then,

$$\partial _t {{\hat{\zeta }}}_{\textrm{vol}}= {\hat{\varepsilon }}_{\textrm{vol}},$$

where ${\hat{\varepsilon }}_{\textrm{vol}}:=\textrm{tr}(\hat{\varepsilon }) \in \varOmega ^{0}(\mathcal {B})$ is the volumetric component of the rate of strain tensor.

Proof

First, we rewrite (80) as $\alpha = \textrm{ln}(\hat{C})$ with $\alpha := 2{\hat{\zeta }}$. Using the power series definition of the matrix exponential, we have that

$$\hat{C} = \textrm{exp}(\alpha ) = I + \alpha + \frac{1}{2!} \alpha ^2 + \frac{1}{3!} \alpha ^3 + \cdots , $$

and similarly using the properties of the matrix exponential

$$\begin{aligned} \hat{C}^{-1} = \textrm{exp}(-\alpha ) = I - \alpha + \frac{1}{2!} \alpha ^2 - \frac{1}{3!} \alpha ^3 + \cdots . \end{aligned}$$

(82)

From the definition of $\hat{C} = G^{-1}\hat{g}$, it follows that that $\dot{\hat{C}} = 2 G^{-1}\hat{g}\hat{\varepsilon }= 2 \hat{C} \hat{\varepsilon },$ and consequently $2 \textrm{tr}(\hat{\varepsilon }) = \textrm{tr}(\hat{C}^{-1}\dot{\hat{C}})$.

The rate of change of $\hat{C}$ can be expressed using the product rule as

$$\begin{aligned} \dot{\hat{C}} = \partial _t \textrm{exp}(\alpha ) = \dot{\alpha } + \frac{1}{2!} (\alpha \dot{\alpha } + \dot{\alpha }\alpha ) + \frac{1}{3!}(\dot{\alpha }\alpha ^2 + \alpha \dot{\alpha }\alpha + \alpha ^2 \dot{\alpha }) + \cdots . \end{aligned}$$

(83)

By multiplying the expressions of $\hat{C}^{-1}$ and $\dot{\hat{C}}$ in (82) and (83), then expanding the product of sums, keeping only the fourth-order terms and lower, one can show that

$$\hat{C}^{-1} \dot{\hat{C}} = \dot{\alpha } + \frac{1}{2!}( \dot{\alpha }\alpha - \alpha \dot{\alpha })+ \frac{1}{3!}( \dot{\alpha }\alpha ^2 - 2\alpha \dot{\alpha }\alpha + \alpha ^2 \dot{\alpha }) +\frac{1}{4!}( \dot{\alpha }\alpha ^3 - 3\alpha \dot{\alpha }\alpha ^2 + 3\alpha ^2\dot{\alpha } \alpha - \alpha ^3 \dot{\alpha }).$$

Taking the trace of the above expression and using the cyclic property of the trace will lead to the vanishing of all terms in the RHS except $\dot{\alpha }$ leading to

$$2 \textrm{tr}(\hat{\varepsilon }) = \textrm{tr}(\hat{C}^{-1}\dot{\hat{C}}) = \textrm{tr}(\dot{\alpha }) = 2 \textrm{tr}(\dot{{\hat{\zeta }}}).$$

Due to the commutativity of the trace with time differentiation, then $\textrm{tr}(\dot{{\hat{\zeta }}}) = \partial _t\textrm{tr}({\hat{\zeta }}) = \partial _t {{\hat{\zeta }}}_{\textrm{vol}}$ which concludes the proof. $\square$

Corollary 4

The rate of change of the deviatoric part of ${\hat{\zeta }}$ is not equal to ${\hat{\varepsilon }}_{\textrm{dev}}$ in general but can be computed by

$$\partial _t {{\hat{\zeta }}}_{\textrm{dev}}= \partial _t {\hat{\zeta }}- \frac{1}{n} {\hat{\varepsilon }}_{\textrm{vol}}.$$

Proof

For $\partial _t {{\hat{\zeta }}}_{\textrm{dev}}= {\hat{\varepsilon }}_{\textrm{dev}}$ to hold, then one must have that $\partial _t {\hat{\zeta }}= \hat{\varepsilon }$ to hold. As shown earlier in the proof of Prop. 6, one has that

$$\partial _t {\hat{\zeta }}= 2 Tf \circ \hat{g}\cdot \hat{\varepsilon }, $$

which for the logarithmic strain definition $f(\hat{g}) = \frac{1}{2}\ln (G^{-1}\hat{g})$ one does not have that $Tf = \frac{1}{2}\hat{g}^{-1}$. $\square$

The above decomposition of ${\hat{\zeta }}$ implies that the constitutive closure problem can be done separately for the volumetric deformation ${{\hat{\zeta }}}_{\textrm{vol}}$ and the isochoric deformation ${{\hat{\zeta }}}_{\textrm{dev}}$. Thus, one can express $\hat{\varPsi }$ as

$$\hat{\varPsi }[{{\hat{\zeta }}}_{\textrm{vol}},{{\hat{\zeta }}}_{\textrm{dev}}] = \int _{\mathcal {B}}({\hat{\psi }}_{\textrm{vol}}({{\hat{\zeta }}}_{\textrm{vol}}) + {\hat{\psi }}_{\textrm{dev}}({{\hat{\zeta }}}_{\textrm{dev}}))\hat{\mu },$$

such that

$$\dot{\hat{\varPsi }} = \int _{\mathcal {B}}\partial _t{{\hat{\zeta }}}_{\textrm{vol}}\wedge \hat{\star } \frac{\partial {\hat{\psi }}_{\textrm{vol}}}{\partial {{\hat{\zeta }}}_{\textrm{vol}}} ( {{\hat{\zeta }}}_{\textrm{vol}}) + \partial _t{{\hat{\zeta }}}_{\textrm{dev}}\ \dot{\wedge }\ \hat{\star }_c\frac{\partial {\hat{\psi }}_{\textrm{dev}}}{\partial {{\hat{\zeta }}}_{\textrm{dev}}} ( {{\hat{\zeta }}}_{\textrm{dev}}),$$

provided the nonlinearity of $\hat{\psi }$ does not couple ${{\hat{\zeta }}}_{\textrm{vol}}$ and ${{\hat{\zeta }}}_{\textrm{dev}}$.

The decomposition above is quite favorable in practice since the volumetric part ${{\hat{\zeta }}}_{\textrm{vol}}$ distinguishes compressible and incompressible materials from each other, as will be discussed later. The remarkable point that Theorem 5 highlights is that ${{\hat{\zeta }}}_{\textrm{vol}}$ is an intrinsic quantity that is independent of the equilibrium configuration G since it is the time integral of another intrinsic quantity. This result is consistent from a geometric point of view since it follows from (64) that

$$\begin{aligned} \partial _t {{\hat{\zeta }}}_{\textrm{vol}}= {\hat{\varepsilon }}_{\textrm{vol}}= \textrm{tr}(\hat{\varepsilon }) = \widehat{\textrm{div}}({\hat{v}}) = \varphi ^*(\textrm{div}(v)), \end{aligned}$$

(84)

which shows that the volumetric deformation ${{\hat{\zeta }}}_{\textrm{vol}}\in \varOmega ^{0}(\mathcal {B})$ is simply the integral of the vector field’s divergence, in the convective representation, which is indeed an intrinsic quantity independent of G.

If we consider the volumetric–deviatoric decomposition of the stress in (67), then $\hat{\tau }$ is computed as the combination of the scalar function ${\hat{\tau }}_{\textrm{vol}}\in \varOmega ^{0}(\mathcal {B})$ and the vector-valued 1 form ${\hat{\tau }}_{\textrm{dev}}\in \varOmega _{\textrm{dev}}^{1}(\mathcal {B};T\mathcal {B})$. However, due to Corollary 4, one cannot in general compute ${\hat{\tau }}_{\textrm{vol}}$ from the gradient of ${\hat{\psi }}_{\textrm{vol}}$ and ${\hat{\tau }}_{\textrm{dev}}$ from the gradient of ${\hat{\psi }}_{\textrm{dev}}$ separately. Although this might seem counter-intuitive, it is in face a natural consequence of the (intrinsic) nonlinearity of $\hat{f}^*$ in the Doyle–Erickson formula (77). In other words, even though one can decompose the power ports $(\hat{\varepsilon },\hat{\star }_c\hat{\tau })$ and $(\partial _t {\hat{\zeta }},\delta _{{\hat{\zeta }}}{\hat{\varPsi }})$ each into two, the nonlinearity of the transformation by $\hat{f}$ and $\hat{f}^*$ mixes these two components.

In the special case of isotropic hyper-elasticity, luckily one can have the aforementioned separation (Sansour 2001). Consider for example the classic Hencky quadratic strain energy (Neff et al. 2016) which is expressed in exterior calculus as

$$\begin{aligned} \hat{\varPsi }[{{\hat{\zeta }}}_{\textrm{vol}},{{\hat{\zeta }}}_{\textrm{dev}}] = \int _{\mathcal {B}}\frac{1}{2}\left( \kappa {{\hat{\zeta }}}_{\textrm{vol}}\wedge \hat{\star } {{\hat{\zeta }}}_{\textrm{vol}}+ 2 \theta {{\hat{\zeta }}}_{\textrm{dev}}\ \dot{\wedge }\ \hat{\star }_c{{\hat{\zeta }}}_{\textrm{dev}}\right) , \end{aligned}$$

(85)

where $\kappa ,\theta \in C^\infty (\mathcal {B})$ denote the bulk and shear moduli respectively. The stress constitutive law, summarized in Fig. 8, takes then the form

$$\begin{aligned} {\hat{\tau }}_{\textrm{vol}}= \frac{\partial {\hat{\psi }}_{\textrm{vol}}}{\partial {{\hat{\zeta }}}_{\textrm{vol}}} = \kappa {{\hat{\zeta }}}_{\textrm{vol}}\qquad \qquad {\hat{\tau }}_{\textrm{dev}}= \frac{\partial {\hat{\psi }}_{\textrm{dev}}}{\partial {{\hat{\zeta }}}_{\textrm{dev}}} = 2 \theta {{\hat{\zeta }}}_{\textrm{dev}}. \end{aligned}$$

(86)

It has been shown in Anand (1979, 1986) that the isotropic Hencky constitutive model above agrees with experiments for moderate principal stretch values. One can extend the application range to high principal stretches by using a nonlinear strain energy $\hat{\psi }$ combined with the nonlinearity of ${\hat{\zeta }}$, e.g., the exponentiated Hencky strain energy (Neff et al. 2015, 2016). The extension of the isotropic Hencky model to the general anisotropic case and to isotropic plasticity can be found in Sansour (2001).

Fig. 8

Convective constitutive relations for finite-strain isotropic hyper-elasticity

5.3 Infinitesimal-Strain Hyper-elasticity

There are numerous hyper-elasticity models in the literature that rely on linear definitions of the strain measure ${\hat{\zeta }}$ to represent finite-strains. However, as shown by Fiala (2011), linear definitions of ${\hat{\zeta }}$ are not consistent with the (non-Euclidean) geometry of $\mathcal {M}(\mathcal {B})$ unless one assumes the deformations to be infinitesimally small.

The theory of small (infinitesimal) strains represents the deformation of the elastic body $\mathcal {B}$ linearized about the reference metric $G\in \mathcal {M}(\mathcal {B})$. Then, the deformation process characterized by the curve $c_{\hat{g}}(t)$ in $\mathcal {M}(\mathcal {B})$ is approximated by a trajectory in $T_G\mathcal {M}(\mathcal {B})$. Consequently, the dynamics takes place in the vector space $T_G\mathcal {M}(\mathcal {B})$ which we identify with $\varOmega ^{1}(\mathcal {B};T\mathcal {B})$.

At a time instant t, the strain tensor is then defined as

$$\begin{aligned} {\hat{\zeta }}_t:= \frac{1}{2}(\hat{C}_t - I), \qquad \qquad \hat{C}_t:= G^{-1} \hat{g}_t\in \varOmega ^{1}(\mathcal {B};T\mathcal {B}), \end{aligned}$$

(87)

where $I \in \varOmega ^{1}(\mathcal {B};T\mathcal {B})\cong \varGamma (T^1_1\mathcal {B})$ denotes the identity tensor. Consequently, using (76), one has that

$$\partial _t {\hat{\zeta }}_t:= \frac{1}{2}\partial _t\hat{C}_t = \frac{1}{2}G^{-1} \partial _t\hat{g}_t = G^{-1} \hat{g}_t \hat{\varepsilon }= \hat{C}_t \hat{\varepsilon }.$$

Thus, for the case of infinitesimal strain, $\hat{f}$ in Prop. 6 is given by the linear map $\hat{f} = \hat{C}$. One can in fact show that linearizing the logarithmic strain measure (80) at $\hat{C} = I$ using Taylor series leads to the above linear strain definition.

Using the linear strain measure above, one can then define different classes of hyper-elastic models by choosing either linear or nonlinear strain energy functions. For instance, the St. Venant–Kirchhoff model has a linear elastic strain energy of the form

$$\hat{\varPsi }[{\hat{\zeta }}] = \int _{\mathcal {B}}\frac{1}{2}{\hat{\zeta }}\ \dot{\wedge }\ \hat{\star }_c\mathbb {E} {\hat{\zeta }},$$

where ${\mathbb {E}}:{\varOmega ^{1}(\mathcal {B};T\mathcal {B})}\rightarrow {\varOmega ^{1}(\mathcal {B};T\mathcal {B})}$ denotes the fourth rank elasticity tensor, which is independent of ${\hat{\zeta }}$ but could vary at different points in $\mathcal {B}$ in the non-homogeneous case. The resulting linear stress–strain constitutive law takes the form $\hat{\tau }= \hat{g}G^{-1} \mathbb {E} {\hat{\zeta }},$ which can be reformulated as the 2-contravariant tensor $\hat{\tau }^\sharp := \hat{g}^{-1}\hat{\tau }$ with local components $\hat{\tau }^{AB}$:

$$\hat{\tau }^\sharp = \hat{\mathbb {E}}(\hat{g}- G) \in \varGamma (T^2_0\mathcal {B}), \qquad \qquad \hat{\tau }^{AB} = \hat{\mathbb {E}}^{ABIJ}(\hat{g}_{IJ} - G_{IJ}) \in C^\infty (\mathcal {B}),$$

with $\hat{\mathbb {E}}:= \frac{1}{2}G^{-1} \mathbb {E} G^{-1}$.

The constitutive law above represents a generalized Hooke’s law with the elasticity tensor $\hat{\mathbb {E}}$ having 21 independent components due to its (major and minor) symmetric properties. The presence of material symmetries reduces the independent components of $\hat{\mathbb {E}}$ with the simplest case being homogeneous isotropic hyper-elasticity with only 2 independent components $\kappa ,\theta$ denoting the bulk and shear modulus. In this simple case, the strain energy is also given by (85) with the stress volumetric and deviatoric components given by (86). While this linear stress–strain constitutive model is physically valid only in the (infinitesimally) small-strain regime, it accounts for arbitrary rigid body motions represented in the nonlinearity of the port-Hamiltonian model in (69-70).

Nonlinear models for isotropic hyper-elasticity, such as the Mooney–Rivlin and neo-Hookean models, are usually expressed in terms of rotational invariants of (87):

$$I_1({\hat{\zeta }}): = \textrm{tr}({\hat{\zeta }}), \qquad I_2({\hat{\zeta }}): = \frac{1}{2}(\textrm{tr}({\hat{\zeta }})^2 - \textrm{tr}({\hat{\zeta }}^2)), \qquad I_3({\hat{\zeta }}): = \det ({\hat{\zeta }}).$$

The strain functional for the (compressible) Mooney–Rivlin model takes the form

$$\begin{aligned} \hat{\varPsi }[{\hat{\zeta }}] = \int _{\mathcal {B}}(c_1 I_1({\hat{\zeta }}) + c_2 I_2({\hat{\zeta }}) + c_3 I_3({\hat{\zeta }}))\hat{\mu }, \end{aligned}$$

(88)

where the material parameters $c_1,c_2\in C^\infty (\mathcal {B})$ characterize distortion response while $c_3\in C^\infty (\mathcal {B})$ characterizes dilatation response. For the special case $c_2 = 0$, one gets the neo-Hookean model.

It is interesting to note that both the Mooney–Rivlin model (88) and the quadratic Hencky model (85) constitute a nonlinear closure relation for the port $(\hat{\varepsilon },\hat{\star }_c\hat{\tau })$. While the Mooney–Rivlin relies on a linear strain definition and a nonlinear strain energy, the Hencky model has a linear strain energy and a (geometrically) nonlinear strain definition. Even though both models have a comparable number of material parameters to fit, the Hencky model has been shown to outperform the Mooney–Rivlin and neo-Hookean models for moderate strains (Neff et al. 2015). This again highlights the importance of taking the geometry of $\mathcal {M}(\mathcal {B})$ into account for defining finite strain.

Remark 9

We emphasize again that in principle, one does need to define strain at all (cf. (Truesdell et al. 1966, Sect. 43)) by defining the strain energy in terms of the state of deformation $\hat{g}$ as in (78). Consequently, the stress law can be perceived as a covector field on $\mathcal {M}(\mathcal {B})$ Rougée (2006). However, the concept of strain and its geometric nature as a mixed tensor comes with numerous useful properties (such as the trace, volumetric–deviatoric decomposition, rotational invariants, and eigenvalues) that are not available to the deformation state $\hat{g}$ intrinsically due to its fully covariant nature.

Remark 10

The pushforward of the mixed tensor field $\hat{C}_t$ onto the reference configuration ${\varphi }_{\textrm{ref}}(\mathcal {B})$ corresponds to the Cauchy–Green deformation tensor, whereas the pushforward of the strain measures (87) and (80) onto ${\varphi }_{\textrm{ref}}(\mathcal {B})$ correspond, respectively, to standard Cauchy–Green strain tensor and the logarithmic strain tensor introduced by Hencky extending the work of Becker, cf. Martin et al. (2018) for a historic review. Furthermore, let the reference mass density function be denoted by ${\rho }_{\textrm{ref}} \in C^\infty (\mathcal {B})$ such that $\hat{\mu }= {\rho }_{\textrm{ref}} \omega _G$ with $\omega _G \in \varOmega ^{n}(\mathcal {B})$ denoting the volume form corresponding to the metric G. Then, the pushforward of the stress ${\rho }_{\textrm{ref}} \hat{\tau }$ onto ${\varphi }_{\textrm{ref}}(\mathcal {B})$ corresponds to the second Piola-Kirchhoff stress tensor represented as a mixed tensor field.

5.4 Incompressible Hyper-elasticity

Fig. 9

Convective constitutive relations of incompressible isotropic hyper-elasticity

Certain elastic materials can be characterized, in certain strain regimes, as being incompressible. A prototypical example is rubber. Such incompressibility is characterized intrinsically, in the spatial and convective representations, respectively, by the geometric conditions

$$\mathcal {L}_{v}{\omega _{g}} = 0 \implies \textrm{div}(v) = 0, \qquad \qquad \mathcal {L}_{{\hat{v}}}{\hat{\omega }_{\hat{g}}} = 0 \implies \widehat{\textrm{div}}({\hat{v}}) = 0,$$

where $\omega _{g}\in \varOmega ^{n}(\mathcal {S}), \hat{\omega }_{\hat{g}}\in \varOmega ^{n}(\mathcal {B})$ denote the volume forms associated with the Riemannian metrics $g$ and $\hat{g}$, respectively.

From the volumetric–deviatoric decomposition (81) of the logarithmic strain measure discussed earlier, (84) implies that the incompressibility constraint is equivalently expressed as

$${\hat{\varepsilon }}_{\textrm{vol}}= 0 \implies {{\hat{\zeta }}}_{\textrm{vol}}= \text {constant},$$

for any finite-strain deformation. Consequently, the volumetric stress ${\hat{\tau }}_{\textrm{vol}}$ is interpreted as a Lagrange multiplier enforcing this constraint and no longer is computed from a constitutive relation. The stress–strain constitutive equations for incompressible hyper-elasticity, summarized in Fig. 9, can be expressed in general as

$$\begin{aligned} \hat{\tau }=&{\hat{\tau }}_{\textrm{vol}}I_n + \frac{\partial {\hat{\psi }}_{\textrm{dev}}}{\partial {{\hat{\zeta }}}_{\textrm{dev}}},\\ \partial _t {{\hat{\zeta }}}_{\textrm{dev}}=&{\hat{\varepsilon }}_{\textrm{dev}},\\ 0 =&{\hat{\varepsilon }}_{\textrm{vol}}. \end{aligned}$$

One can also model near-incompressibility effects trivially by choosing ${\hat{\tau }}_{\textrm{vol}}= \frac{\partial {f}_{\textrm{pen}}}{\partial {{\hat{\zeta }}}_{\textrm{vol}}}({{\hat{\zeta }}}_{\textrm{vol}})$ where ${{f}_{\textrm{pen}}}:{C^\infty (\mathcal {B})}\rightarrow {C^\infty (\mathcal {B})}$ denotes the penalty function of choice.

The aforementioned constitutive equation for incompressible hyper-elasticity can be recovered from the compressible models discussed earlier as the volumetric parameters of the model (given by the bulk modulus $\kappa$ in the isotropic case) tend to zero. In addition, from the trace property $\textrm{tr}(\textrm{ln}(\hat{C})) = \textrm{ln}(\det (\hat{C})),$ one can see that the incompressibility constraint expressed in terms of $\hat{C}$ is given by

$$\textrm{ln}(\det (\hat{C})) = 0 \implies \det (\hat{C}) = 1.$$

Using the linear strain measure (87) then implies that the rotational invariant $I_3=\det ({\hat{\zeta }}) = 0.$ Thus, incompressible versions of isotropic nonlinear models, such as the Mooney–Rivlin model, use (87); the invariant $\det ({\hat{\zeta }})$ is the variable that vanishes and not $\textrm{tr}({\hat{\zeta }})$.

5.5 Viscous Fluid Flow

The last type of constitutive relations that we discuss is viscous fluid flow, in particular Newtonian fluids. Newtonian fluids have the simplest mathematical form as they are characterized by an isotropic linear constitutive relation similar to (86). As illustrated in Fig. 10, the constitutive model for Newtonian fluids, neglecting thermal effects, is characterized by

• 1. a linear resistive relation between $({\hat{\varepsilon }}_{\textrm{dev}},{\hat{\tau }}_{\textrm{dev}})$ representing the viscous response to distortion,

• 2. a combined resistive and storage relation between $({\hat{\varepsilon }}_{\textrm{vol}},{\hat{\tau }}_{\textrm{vol}})$ representing a visco-elastic response to volumetric deformation.

Fig. 10

Convective constitutive relations of Newtonian fluids

Let the convective internal energy functional be given by

$$\begin{aligned} \hat{U}[{{\hat{\zeta }}}_{\textrm{vol}};\hat{\mu }] = \int _{\mathcal {B}}\hat{\mathcal {U}}({{\hat{\zeta }}}_{\textrm{vol}}) \hat{\mu }, \end{aligned}$$

(89)

with ${\hat{\mathcal {U}}}:{\varOmega ^{0}(\mathcal {B})}\rightarrow {\varOmega ^{0}(\mathcal {B})}$ denoting the specific internal energy density function that describes the equation of state of the fluid. The internal energy’s state is given by the volumetric strain ${{\hat{\zeta }}}_{\textrm{vol}}$ and depends parametrically on $\hat{\mu }$. The constitutive stress law for Newtonian fluids is then expressed in the convective representation as

$$\begin{aligned} {\hat{\tau }}_{\textrm{dev}}= 2 \theta {\hat{\varepsilon }}_{\textrm{dev}}, \qquad \qquad {\hat{\tau }}_{\textrm{vol}}= \underbrace{\kappa {\hat{\varepsilon }}_{\textrm{vol}}}_{{\hat{\tau }}_{\textrm{b}}} + \underbrace{\frac{\partial \hat{\mathcal {U}}}{\partial {{\hat{\zeta }}}_{\textrm{vol}}}({{\hat{\zeta }}}_{\textrm{vol}})}_{{\hat{\tau }}_{\textrm{int}}}, \qquad \qquad \partial {{\hat{\zeta }}}_{\textrm{vol}}= {\hat{\varepsilon }}_{\textrm{vol}}, \end{aligned}$$

(90)

with $\kappa ,\theta \in C^\infty (\mathcal {B})$ denoting the bulk and shear viscosity coefficients, respectively, which are in general state dependent. The variational derivative of the internal energy functional is given by $\delta _{{{\hat{\zeta }}}_{\textrm{vol}}}{\hat{U}}({{\hat{\zeta }}}_{\textrm{vol}}) = \hat{\star } \frac{\partial \hat{\mathcal {U}}}{\partial {{\hat{\zeta }}}_{\textrm{vol}}}({{\hat{\zeta }}}_{\textrm{vol}})= \frac{\partial \hat{\mathcal {U}}}{\partial {{\hat{\zeta }}}_{\textrm{vol}}}({{\hat{\zeta }}}_{\textrm{vol}}) \hat{\mu }\in \varOmega ^{n}(\mathcal {S}).$

Using $\dot{\hat{U}} = \langle \delta _{{{\hat{\zeta }}}_{\textrm{vol}}}{\hat{U}} | \partial _t {{\hat{\zeta }}}_{\textrm{vol}} \rangle _{\mathcal {B}} = \langle \hat{\star }{\hat{\tau }}_{\textrm{int}} | {\hat{\varepsilon }}_{\textrm{vol}} \rangle _{\mathcal {B}}$, the above constitutive law has a power balance expressed as

$$\begin{aligned} \langle \hat{\star }_c\hat{\tau } | \hat{\varepsilon } \rangle _{\mathcal {B}} =&\langle \hat{\star } {\hat{\tau }}_{\textrm{vol}} | {\hat{\varepsilon }}_{\textrm{vol}} \rangle _{\mathcal {B}} + \langle \hat{\star }_c{\hat{\tau }}_{\textrm{dev}} | {\hat{\varepsilon }}_{\textrm{dev}} \rangle _{\mathcal {B}} \\ =&\dot{\hat{U}} + \int _{\mathcal {B}}\kappa {\hat{\varepsilon }}_{\textrm{vol}}\wedge \hat{\star } {\hat{\varepsilon }}_{\textrm{vol}}+ 2 \theta {\hat{\varepsilon }}_{\textrm{dev}}\ \dot{\wedge }\ \hat{\star }_c{\hat{\varepsilon }}_{\textrm{dev}}, \end{aligned}$$

which characterizes the storage of internal energy $\hat{U}$ and two Rayleigh dissipation functions.

Finally, we conclude this section by deriving the spatial counterpart of the constitutive law (90) for $({\varepsilon }_{\textrm{dev}},{\tau }_{\textrm{dev}})$ and $({\varepsilon }_{\textrm{vol}},{\tau }_{\textrm{vol}})$ which is more common for fluids as presented in Sect. 4. We start with re-expressing the internal energy (89) in terms of other state variables that are more common, namely the mass density and specific volume.

Proposition 7

Let $\hat{\rho },\hat{\nu } \in \varOmega ^{0}(\mathcal {B})$ denote the convective mass density and specific volume functions, defined such that the mass form $\hat{\mu }\in \varOmega ^{n}(\mathcal {B})$ can be expressed as

$$\hat{\mu }= \hat{\rho }\hat{\omega }_{\hat{g}}, \qquad \text {or} \qquad \hat{\omega }_{\hat{g}}= \hat{\nu } \hat{\mu },$$

where $\hat{\omega }_{\hat{g}}\in \varOmega ^{n}(\mathcal {B})$ is the volume form associated to the convective metric $\hat{g}\in \mathcal {M}(\mathcal {B})$. The volumetric strain is related to these two variables by

$${{\hat{\zeta }}}_{\textrm{vol}}= \ln (\hat{\nu }) = - \ln (\hat{\rho }).$$

Proof

Recall that $\hat{\omega }_{\hat{g}}= \varphi ^*(\omega _{g})$ where $\omega _{g}\in \varOmega ^{n}(\mathcal {S})$ is the spatial volume form associated with $g$. From the Lie derivative definition, one has that

$$\partial _t \hat{\omega }_{\hat{g}}= \varphi ^*(\mathcal {L}_{v}{\omega _{g}}) = \mathcal {L}_{\varphi ^*(v)}{\varphi ^*(\omega _{g})} = \mathcal {L}_{{\hat{v}}}{\hat{\omega }_{\hat{g}}} = \widehat{\textrm{div}}({\hat{v}}) \hat{\omega }_{\hat{g}}.$$

Since $\hat{\mu }$ is constant, then $\partial _t \hat{\omega }_{\hat{g}}= \partial _t({\hat{v}})\hat{\mu }= \frac{\partial _t{\hat{v}}}{{\hat{v}}}$ which implies that $\partial _t{\hat{v}}= {\hat{v}}\widehat{\textrm{div}}({\hat{v}})$. From Th. 5 and the fact that $\partial _t {{\hat{\zeta }}}_{\textrm{vol}}= \frac{\partial _t{\hat{v}}}{{\hat{v}}}$, then one has that ${{\hat{\zeta }}}_{\textrm{vol}}=\ln (\hat{\nu })$. Using the properties of the natural logarithm, it is straightforward that ${{\hat{\zeta }}}_{\textrm{vol}}=\text {ln}(\frac{1}{\hat{\rho }}) = -\ln (\hat{\rho })$. $\square$

With the convective mass density as the state variable, one can rewrite the internal energy (89) as

$$\begin{aligned} \hat{U}[\hat{\rho };\hat{g}] = \int _{\mathcal {B}}\hat{\mathcal {U}}(\hat{\rho }) \hat{\rho }\hat{\omega }_{\hat{g}}, \end{aligned}$$

(91)

with $\hat{\mathcal {U}}$ now being a function of $\hat{\rho }$ with a slight abuse of notation. Using $\partial _t \hat{\rho }= - \hat{\rho }\partial _t {{\hat{\zeta }}}_{\textrm{vol}}= - \hat{\rho }{\hat{\varepsilon }}_{\textrm{vol}}$, the rate of change of $\hat{U}$ now takes the form

$$\dot{\hat{U}} = \int _{\mathcal {B}}\partial _t \hat{\mathcal {U}}(\hat{\rho }) \hat{\rho }\hat{\omega }_{\hat{g}}=\int _{\mathcal {B}}\partial _t \hat{\rho }\frac{\partial \hat{\mathcal {U}}}{\partial \hat{\rho }}(\hat{\rho }) \hat{\mu }= \int _{\mathcal {B}}-{\hat{\varepsilon }}_{\textrm{vol}}\hat{\rho }\frac{\partial \hat{\mathcal {U}}}{\partial \hat{\rho }}(\hat{\rho }) \hat{\mu }.$$

Since $\dot{\hat{U}}= \langle \hat{\star }{\hat{\tau }}_{\textrm{int}} | {\hat{\varepsilon }}_{\textrm{vol}} \rangle _{\mathcal {B}}$, then ${\hat{\tau }}_{\textrm{int}}$ in (90) can be expressed as ${\hat{\tau }}_{\textrm{int}} = - \hat{\rho }\frac{\partial \hat{\mathcal {U}}}{\partial \hat{\rho }}(\hat{\rho }).$

Now let $U:= \hat{U}\circ \varphi$ denote the spatial representation of the internal energy functional defined such that $U[\rho ;g] = \hat{U}[\varphi ^*(\rho );\varphi ^*(g)]$ and expressed as

$$\begin{aligned} U[\rho ;g] = \int _{\mathcal {S}}\mathcal {U}(\rho ) \rho \omega _{g}, \end{aligned}$$

(92)

with the spatial mass density $\rho \in \varOmega ^{0}(\mathcal {S})$ as the state variable and ${\mathcal {U}}:{\varOmega ^{0}(\mathcal {S})}\rightarrow {\varOmega ^{0}(\mathcal {S})}$ denoting the spatial specific internal energy. Using Reynold’s transport theorem, the rate of change of U is given by

$${\dot{U}} = \frac{\textrm{d}}{\textrm{d}t} \int _{\mathcal {S}}\mathcal {U}(\rho ) \rho \omega _{g}=\int _{\mathcal {S}}\textrm{D}_t (\mathcal {U}(\rho ) \rho \omega _{g})=\int _{\mathcal {S}}\textrm{D}_t (\mathcal {U}(\rho )) \rho \omega _{g},$$

which follows from the mass balance $\textrm{D}_t(\rho \omega _{g}) = \textrm{D}_t \mu = 0$ with $\textrm{D}_t:= \partial _t + \mathcal {L}_{v}$ being a shorthand notation. From the non-equilibrium thermodynamics principle $\textrm{D}_t (\mathcal {U}(\rho )) = \frac{\partial \mathcal {U}}{\partial \rho }(\rho )\textrm{D}_t\rho$ and mass balance $\textrm{D}_t\rho = -\rho \textrm{div}(v)$Málek and Pruša (2018), one can express ${\dot{U}}$ as

$${\dot{U}} = \int _{\mathcal {S}}\rho \frac{\partial \mathcal {U}}{\partial \rho }(\rho ) \textrm{D}_t\rho \omega _{g}= - \int _{\mathcal {S}}p \textrm{div}(v)\omega _{g},$$

where $p:= \rho ^2 \frac{\partial \mathcal {U}}{\partial \rho }(\rho ) \in C^\infty (\mathcal {S})$ is the thermodynamic pressure function. Given that ${\varepsilon }_{\textrm{vol}}=\textrm{div}(v)$, then for ${\dot{U}} = \langle \star {\tau }_{\textrm{int}} | {\varepsilon }_{\textrm{vol}} \rangle _{\mathcal {S}}$ to hold then one has that ${\tau }_{\textrm{int}} = - \frac{p}{\rho }$.

The spatial counterpart of the constitutive stress law for Newtonian fluids (90) is then expressed as

$$\begin{aligned} {\tau }_{\textrm{dev}}= 2 \theta {\varepsilon }_{\textrm{dev}}, \qquad \qquad {\tau }_{\textrm{vol}}= \kappa {\varepsilon }_{\textrm{vol}}{\tau }_{\textrm{b}} - p, \qquad \qquad \textrm{D}_t\rho = -\rho {\varepsilon }_{\textrm{vol}}, \end{aligned}$$

(93)

with $\kappa ,\theta \in C^\infty (\mathcal {S})$ denoting the bulk and shear viscosity coefficients considered now functions on $\mathcal {S}$. The above constitutive law has a power balance expressed as

$$\begin{aligned} \langle \star _c\tau | \varepsilon \rangle _{\mathcal {S}}&=\ \langle \star {\tau }_{\textrm{vol}} | {\varepsilon }_{\textrm{vol}} \rangle _{\mathcal {S}} + \langle \star _c{\tau }_{\textrm{dev}} | {\varepsilon }_{\textrm{dev}} \rangle _{\mathcal {S}} \\&=\ {\dot{U}} + \int _{\mathcal {S}}\kappa {\varepsilon }_{\textrm{vol}}\wedge \star {\varepsilon }_{\textrm{vol}}+ 2 \theta {\varepsilon }_{\textrm{dev}}\ \dot{\wedge }\ \star _c{\varepsilon }_{\textrm{dev}}. \end{aligned}$$

By adding the constitutive relation (90) to the port-Hamiltonian system (74-75), one has a port-Hamiltonian representation of the Navier–Stokes equations for Newtonian fluids. The total energy balance after adding the stress constitutive law can be expressed as

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t} ({{H}}_{\textrm{k}}+ U) = \int _{{\partial \mathcal {S}}} v|_{{\partial \mathcal {S}}} \ \dot{\wedge }\ (\mathcal {T}-\frac{1}{2}\iota _v\mathcal {M})|_{{\partial \mathcal {S}}} - \int _{\mathcal {S}}\kappa {\varepsilon }_{\textrm{vol}}\wedge \star {\varepsilon }_{\textrm{vol}}+ 2 \theta {\varepsilon }_{\textrm{dev}}\ \dot{\wedge }\ \star _c{\varepsilon }_{\textrm{dev}}. \end{aligned}$$

The above power balance states that the rate of change of total energy within the domain $\mathcal {S}$ is equal to the difference between the power flow through the boundary ${\partial \mathcal {S}}$, due to stress and momentum flux, and the power dissipated in the domain, due to bulk and shear viscosity.

Remark 11

i) The stress $\tau \in \varOmega ^{1}(\mathcal {S};T\mathcal {S}) \cong \varGamma (T^1_1\mathcal {S})$ is the spatial representation of the intensive Rougee stress discussed in Sect. 5 and is related to the standard Cauchy stress tensor $\sigma \in \varGamma (T^1_1\mathcal {S})$ by $\sigma = \rho \tau$. The final constitutive law of compressible Newtonian fluids in terms of $\sigma$ can then be summarized as

$$\sigma = \rho \tau = \rho ({\tau }_{\textrm{vol}}I_n + {\tau }_{\textrm{dev}}) = (-p + {\kappa }_{\textrm{dyn}} {\varepsilon }_{\textrm{vol}}) I_n + {\theta }_{\textrm{dyn}} {\varepsilon }_{\textrm{dev}},$$

where ${\kappa }_{\textrm{dyn}}:= \rho \kappa \in C^\infty (\mathcal {S})$ denotes the dynamic bulk viscosity and ${\theta }_{\textrm{dyn}}:= \rho \theta \in C^\infty (\mathcal {S})$ denotes the dynamic shear viscosity. Furthermore, by introducing the mechanical pressure function $p_m:= -\frac{1}{n} \textrm{tr}(\sigma )\in C^\infty (\mathcal {S})$, then one recovers the standard relation between mechanical and thermodynamic functions $p_m = p - {\kappa }_{\textrm{dyn}} \textrm{div}(v).$

ii) Incorporating incompressibility in the port-Hamiltonian model is done identically to the procedure explained in Sect. 5.4. The volumetric component of the stress ${\tau }_{\textrm{vol}}$ is then treated as a Lagrange multiplier enforcing the incompressibility constraint ${\varepsilon }_{\textrm{vol}}= \textrm{div}(v) = 0$. Furthermore, the constitutive model (90) can be easily extended to model more complex fluid behavior. For instance, representing the shear viscosity by a fourth-order tensor instead of a function allows modeling anisotropic Newtonian fluids, while representing the bulk and shear viscosity as nonlinear functions of ${\varepsilon }_{\textrm{vol}}$ and ${\varepsilon }_{\textrm{dev}}$, respectively, allows modeling non-Newtonian fluids. As for modeling the transfer of the dissipated energy to the thermal domain, the interested reader is referred to Califano et al. (2022); Málek and Pruša (2018) and references therein for a discussion of the topic.

6 Conclusion

In this paper, we presented decomposed port-Hamiltonian models that represent the governing equations of continuum mechanics, complementing our recent work (Rashad et al. 2023). These models provide a geometric coordinate-free insight by using bundle-valued differential forms to mathematically represent physical variables. Thanks to the elegant machinery of exterior calculus, the models were derived using Hamiltonian reduction theory from first principles. A distinguishing feature of our work is the port-based philosophy of tearing the overall system into its constituent energetic units. An immediate benefit of this approach is that one does not have to consider all physical variables of the whole system at once. Consequently, this simplifies greatly the application of Hamiltonian reduction techniques.

The derived port-Hamiltonian models expose the rich geometric, topological and energetic structure underlying the theory of continuum mechanics. Our work lays the foundations for exploiting this rich structure in analysis and control of distributed parameter systems, in addition to the development of efficient numerical algorithms for discretization and model order reduction of these dynamic equations that preserve this structure at the discrete level.

A promising future extension of our theoretical framework could be to derive simplified port-Hamiltonian models tailored for engineering purposes, for example, classical linear beam formulations (of Euler–Bernoulli or Timoshenko) or geometrically exact beam formulations (of Von Karman). Deriving such reduced port-Hamiltonian models from our generic ones for continuum mechanics will target a broader audience and provide more intuition for the benefits of the port-Hamiltonian framework.

Appendix

1.1 Proof of Prop. 2

(i) Let $\delta \tilde{\varphi }\in T_\varphi \mathscr {C}$ be the tangent vector to the curve $s\mapsto \varphi _s$ at $\varphi \in \mathscr {C}$ at $s=0$. Then $\delta \mu \in \varOmega ^{n}(\mathcal {S})$ is the tangent vector to the induced curve $s\mapsto \mu _s:= \varphi _{s,*}(\tilde{\mu })$ at $\mu \in \varOmega ^{n}(\mathcal {S})$ at $s=0$. One has from the Lie derivative definition that

$$\delta \mu := \left. \frac{d}{ds}\right| _{s=0}\varphi _{s,*}(\tilde{\mu }) = - \mathcal {L}_{\delta \varphi }{\mu },$$

with $\delta \varphi := \varphi _{\textrm{f},*}(\delta \tilde{\varphi }) = \delta \tilde{\varphi }\circ \varphi ^{-1}\in \varOmega ^{0}(\mathcal {S};T\mathcal {S})$.

The variation $\delta \mathcal {M}\in \varOmega ^{n}(\mathcal {S};T^*\mathcal {S})$ consists of two components; one due to the variation of $\varphi$ and another due to the variation of $\tilde{\mathcal {M}}$, such that: $\delta \mathcal {M}:= \delta \mathcal {M}_1 + \delta \mathcal {M}_2$. For the later one due to varying $\tilde{\mathcal {M}}$, it is simply given by $\delta \mathcal {M}_2:= \left. \frac{d}{ds}\right| _{s=0}\varphi _{\textrm{f},*}(\tilde{\mathcal {M}}+s\delta \tilde{\mathcal {M}}) =\varphi _{\textrm{f},*}(\delta \tilde{\mathcal {M}}).$ As for $\delta \mathcal {M}_1$, it is more involved due to the intuition mentioned preceding (22). From Remark 3, one can identify $\mathcal {M}$ with $\mu \otimes v^\flat$ or equivalently $\omega _{g}\otimes m^\flat$, where $\omega _{g}\in \varOmega ^{n}(\mathcal {S})$ and $m^\flat \in \varOmega ^{1}(\mathcal {S})$ are the volume top form and intensive momentum one-form, respectively. Since $\omega _{g}$ is constant, one can express $\delta \mathcal {M}_1 = ^{g}\star _c(\delta m) = \omega _{g}\otimes \delta m^\flat$, with $\delta m^\flat \in \varOmega ^{1}(\mathcal {S})$ and $^{g}\star _c$ denotes the (constant) Hodge star operator constructed from the volume form $\omega _{g}$. It has been shown in Mazer and Ratiu (1989) that the vector field $\delta m:= g^{-1}\cdot \delta m^\flat \in \varGamma (T\mathcal {S})$ can be related to $\delta \varphi$ by

$$\delta m = - \textrm{div}(\delta \varphi m) - \nabla _{\delta \varphi } m = - \textrm{div}(m\otimes \delta \varphi ),$$

where $m\otimes \delta \varphi$ represents a momentum flux. Furthermore, it has been shown in (Rashad et al. 2023, Th. 1) that the divergence of a two-covariant tensor is equivalent in exterior calculus to the exterior covariant derivative of a related covector-valued $n-1$ form. For the case of $m\otimes \delta \varphi$, one can show that $\delta \mathcal {M}_1 = - \textrm{d}_{\nabla }( ^{g}\star _c(m^\flat \otimes \delta \varphi ))$. Finally, one can show that the covector-valued $n-1$ form $^{g}\star _c(m^\flat \otimes \delta \varphi )$ is equivalent to the extensive momentum flux expressed as $\iota _{\delta \varphi } \mathcal {M}:= \iota _{\delta \varphi }\omega _{g}\otimes m^\flat \in \varOmega ^{n-1}(\mathcal {S};T^*\mathcal {S})$, which concludes the first part of the proof.

(ii) From the chain rule of functionals and (21,33), one has that

$$\begin{aligned} \textrm{D}_{\tilde{\mathcal {M}}} \tilde{\mathscr {F}}(\varphi ,\tilde{\mathcal {M}}) \cdot \delta \tilde{\mathcal {M}}=&\left. \frac{d}{ds}\right| _{s=0}\tilde{\mathscr {F}}[\varphi ,\tilde{\mathcal {M}}+ s \delta \tilde{\mathcal {M}}]\\ =&\left. \frac{d}{ds}\right| _{s=0}\mathscr {F}[\varphi _*( \tilde{\mu }),\varphi _{\textrm{f},*}(\tilde{\mathcal {M}}+ s \delta \tilde{\mathcal {M}})] = \textrm{D}_{\mathcal {M}}\mathscr {F}(\mu ,\mathcal {M}) \cdot \delta \mathcal {M}_2. \end{aligned}$$

Using the change of variables formula and the expression for $\delta \mathcal {M}_2$ in (i) above, one has that

$$\begin{aligned} \textrm{D}_{\mathcal {M}}\mathscr {F}(\mu ,\mathcal {M}) \cdot \delta \mathcal {M}_2&= \int _{\mathcal {S}}\delta \mathcal {M}_2 \ \dot{\wedge }\ \delta _{\mathcal {M}}{\mathscr {F}} = \int _{\mathcal {B}}\varphi ^*(\delta \mathcal {M}_2 \ \dot{\wedge }\ \delta _{\mathcal {M}}{\mathscr {F}} ) \\&= \int _{\mathcal {B}}\varphi _{\textrm{f}}^*(\delta \mathcal {M}_2) \ \dot{\wedge }\ \varphi _{\textrm{f}}^*(\delta _{\mathcal {M}}{\mathscr {F}}) = \int _{\mathcal {B}}\delta \tilde{\mathcal {M}}\ \dot{\wedge }\ \varphi _{\textrm{f}}^*(\delta _{\mathcal {M}}{\mathscr {F}}) . \end{aligned}$$

Comparing the above two expressions and using (20) yields that $\delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}} = \varphi _{\textrm{f}}^*(\delta _{\mathcal {M}}{\mathscr {F}})$.

Similarly, from the chain rule of functionals and (22,33), one has that

$$\begin{aligned} \textrm{D}_{\varphi } \tilde{\mathscr {F}}(\varphi ,\tilde{\mathcal {M}}) \cdot \delta \tilde{\varphi }=&\left. \frac{d}{ds}\right| _{s=0}\tilde{\mathscr {F}}[\varphi _s,\tilde{\mathcal {M}}_s] = \left. \frac{d}{ds}\right| _{s=0}\mathscr {F}[\varphi _{s,*}( \tilde{\mu }),(\varphi _s)_{\textrm{f},*}(\tilde{\mathcal {M}}_s)] \nonumber \\ =&\textrm{D}_{\mu }\mathscr {F}(\mu ,\mathcal {M}) \cdot \delta \mu + \textrm{D}_{\mathcal {M}}\mathscr {F}(\mu ,\mathcal {M}) \cdot \delta \mathcal {M}_1. \end{aligned}$$

(94)

Using (34) and the expressions for $\delta \mu ,\delta \mathcal {M}_1$ in (i) above, one can rewrite (94) as

$$\begin{aligned} \textrm{D}_{\mu }\mathscr {F}(\mu ,\mathcal {M}) \cdot \delta \mu + \textrm{D}_{\mathcal {M}}\mathscr {F}(\mu ,\mathcal {M}) \cdot \delta \mathcal {M}_1 = - \int _{\mathcal {S}}\delta _{\mu }{\mathscr {F}} \wedge \mathcal {L}_{\delta \varphi }{\mu } + \delta _{\mathcal {M}}{\mathscr {F}} \ \dot{\wedge }\ \textrm{d}_{\nabla }\iota _{\delta \varphi } \mathcal {M}. \nonumber \\ \end{aligned}$$

(95)

Consider the following identities between scalar-valued and bundle-valued forms

$$\iota _{\delta \varphi }(\textrm{d}f) \wedge \mu = \delta \varphi \ \dot{\wedge }\ (\mu \otimes \textrm{d}f) \in \varOmega ^{n}(\mathcal {S}), \qquad \iota _{\delta \varphi } (f \wedge \mu ) = \delta \varphi \ \dot{\wedge }\ (f \wedge \mu ) \in \varOmega ^{n-1}(\mathcal {S}),$$

for any $f\in C^\infty (\mathcal {S})$, with $\mu$ in the second LHS is considered as a top-form and in the RHS as a covector-valued $n-1$ form. Using the identities above, the Leibniz rule and Cartan’s formula for the Lie derivative, and Stokes theorem, one has that

$$\begin{aligned} \int _{\mathcal {S}}\delta _{\mu }{\mathscr {F}} \wedge \mathcal {L}_{\delta \varphi }{\mu }&= - \int _{\mathcal {S}}\mathcal {L}_{\delta \varphi }{(\delta _{\mu }{\mathscr {F}})} \wedge \mu + \int _{\mathcal {S}}\mathcal {L}_{\delta \varphi }{(\delta _{\mu }{\mathscr {F}}\wedge \mu )} \nonumber \\&= - \int _{\mathcal {S}}\iota _{\delta \varphi } \textrm{d}{(\delta _{\mu }{\mathscr {F}})} \wedge \mu + \int _{{\partial \mathcal {S}}} \iota _{\delta \varphi }{(\delta _{\mu }{\mathscr {F}}\wedge \mu )}|_{{\partial \mathcal {S}}} \nonumber \\&= - \int _{\mathcal {S}}\delta \varphi \ \dot{\wedge }\ (\mu \otimes \textrm{d}\delta _{\mu }{\mathscr {F}}) + \int _{{\partial \mathcal {S}}} \delta \varphi |_{{\partial \mathcal {S}}} \ \dot{\wedge }\ \delta _{\mu }{\mathscr {F}}\wedge \mu |_{{\partial \mathcal {S}}}. \end{aligned}$$

(96)

Now consider the following identities for any vector field $\eta \in \varGamma (T\mathcal {S})$

$$\eta \ \dot{\wedge }\ \iota _{\delta \varphi }\mathcal {M}= \delta \varphi \ \dot{\wedge }\ \iota _{\eta }\mathcal {M}\in \varOmega ^{n}(\mathcal {S}), \qquad \nabla \eta \ \dot{\wedge }\ \iota _{\delta \varphi } \mathcal {M}= \delta \varphi \ \dot{\wedge }\ \mu \otimes (\nabla \eta \ \dot{\wedge }\ v^\flat ) \in \varOmega ^{n}(\mathcal {S}), $$

with $\mathcal {M}= \mu \otimes v^\flat$. Using the above identities, the Leibniz rule for the exterior covariant derivative, and Stokes theorem, one has that

$$\begin{aligned} \int _{\mathcal {S}}\delta _{\mathcal {M}}{\mathscr {F}} \ \dot{\wedge }\ \textrm{d}_{\nabla }\iota _{\delta \varphi } \mathcal {M}=&- \int _{\mathcal {S}}\nabla (\delta _{\mathcal {M}}{\mathscr {F}}) \ \dot{\wedge }\ \iota _{\delta \varphi } \mathcal {M}+ \int _{\mathcal {S}}\textrm{d}( \delta _{\mathcal {M}}{\mathscr {F}}\ \dot{\wedge }\ \iota _{\delta \varphi } \mathcal {M}) \nonumber \\ =&- \int _{\mathcal {S}}\delta \varphi \ \dot{\wedge }\ \mu \otimes (\nabla (\delta _{\mathcal {M}}{\mathscr {F}})\ \dot{\wedge }\ v^\flat ) + \int _{{\partial \mathcal {S}}} \delta \varphi |_{{\partial \mathcal {S}}} \ \dot{\wedge }\ \iota _{\delta _{\mathcal {M}}{\mathscr {F}}} \mathcal {M}|_{{\partial \mathcal {S}}}. \end{aligned}$$

(97)

Finally, by substituting (96,97) in (95), combining (95) with (94) and using the change of variables formula, the expressions of $\delta _{\varphi }{\tilde{\mathscr {F}}}$ and $\delta ^\cup _{\varphi }{\tilde{\mathscr {F}}}$ are concluded by comparing the result with (19).

1.2 Proof of Prop. 4

(i) Recall the curve $s\mapsto \varphi _s$ in $\mathscr {C}$ introduced in Prop. 2. One has that $\delta \hat{g}:= \left. \frac{d}{ds}\right| _{s=0}\hat{g}_s$, where $\hat{g}_s:= \varphi _s^*(g)$ is the induced curve on $\mathcal {M}(\mathcal {B})$. From the Lie derivative properties, it follows that $\delta \hat{g}= \left. \frac{d}{ds}\right| _{s=0}\hat{g}_s = \varphi ^*(\mathcal {L}_{\delta \varphi }{g}) = \mathcal {L}_{\delta \hat{\varphi }}{\hat{g}}$, where $\delta \hat{\varphi }= T \varphi ^{-1} \circ \delta \tilde{\varphi }$.

The variation $\delta \hat{\mathcal {M}}= \delta \hat{\mathcal {M}}_1 + \delta \hat{\mathcal {M}}_2 \in \varOmega ^{n}(\mathcal {B};T^*\mathcal {B})$ also consists of two components; one due to the variation of $\varphi$ and another due to $\tilde{\mathcal {M}}$, where the latter is given by $\delta \hat{\mathcal {M}}_2:= \varphi _{\textrm{v}}^*(\delta \tilde{\mathcal {M}})$. The first part is defined as $\delta \hat{\mathcal {M}}_1:= \left. \frac{d}{ds}\right| _{s=0}\hat{\mathcal {M}}_s$, where $s\mapsto \hat{\mathcal {M}}_s$ is the induced curve in $\varOmega ^{n}(\mathcal {B};T^*\mathcal {B})$ defined in local components as $\hat{\mathcal {M}}_s:= \hat{\mu } ^sF^i_I \delta \tilde{\varphi }_i \otimes E^I$, where $\delta \tilde{\varphi }_i\in \varOmega ^{0}(\mathcal {B})$ are the components of $\delta \tilde{\varphi }$ and $^sF^i_I \in \varOmega ^{0}(\mathcal {B})$ are the components of $^sF:= T\varphi _s$. Using the identities $\left. \frac{d}{ds}\right| _{s=0} ^sF^i_I:= (\tilde{\nabla }\delta \tilde{\varphi })^i_I$ and $(\tilde{\nabla }\delta \tilde{\varphi })^i_I = F^i_J (\hat{\nabla }\delta \hat{\varphi })^J_I$ with $\tilde{\nabla }$ denoting the material covariant derivative (cf. Simo et al. (1988); Rashad et al. (2023)), then one has that

$$\begin{aligned} \delta \hat{\mathcal {M}}_1:= \left. \frac{d}{ds}\right| _{s=0}\hat{\mathcal {M}}_s =&\hat{\mu }\delta \tilde{\varphi }_i \left. \frac{d}{ds}\right| _{s=0} ^sF^i_I \otimes E^I = \hat{\mu }\delta \tilde{\varphi }_i(\tilde{\nabla }\delta \tilde{\varphi })^i_I\otimes E^I\\ =&\hat{\mu }\delta \tilde{\varphi }_i F^i_J (\hat{\nabla }\delta \hat{\varphi })^J_I \otimes E^I = \hat{\mu }\delta \hat{\varphi }_J (\hat{\nabla }\delta \hat{\varphi })^J_I \otimes E^I, \end{aligned}$$

which can be written as $\delta \hat{\mathcal {M}}_1 = \hat{\mu }\otimes (\hat{\nabla }\delta \hat{\varphi }\ \dot{\wedge }\ \hat{v}^\flat )$.

(ii) The derivation of $\delta _{\tilde{\mathcal {M}}}{\tilde{\mathscr {F}}} = \varphi _{\textrm{v},*}(\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}})$ follows identically its spatial counterpart in Prop. 2 (ii). For the other variational derivatives, one has from the chain rule similar to (94) that

$$\textrm{D}_{\varphi } \tilde{\mathscr {F}}(\varphi ,\tilde{\mathcal {M}}) \cdot \delta \tilde{\varphi }= \textrm{D}_{\hat{g}}\hat{\mathscr {F}}(\hat{g},\hat{\mathcal {M}}) \cdot \delta \hat{g}+ \textrm{D}_{\hat{\mathcal {M}}}\hat{\mathscr {F}}(\hat{g},\hat{\mathcal {M}}) \cdot \delta \hat{\mathcal {M}}_1,$$

where the RHS can be further expanded using identities $\frac{1}{2}\mathcal {L}_{{\hat{v}}}{\hat{g}} = \textrm{sym}(\hat{\textrm{d}}_{\hat{\nabla }}\hat{v}^\flat )$ and $\hat{\nabla }\hat{\eta } \ \dot{\wedge }\ \iota _{\delta \hat{\varphi }} \hat{\mathcal {M}}= \delta \hat{\varphi }\ \dot{\wedge }\ \hat{\mu }\otimes (\hat{\nabla }\hat{\eta }\ \dot{\wedge }\ \hat{v}^\flat )$ as

$$\begin{aligned} \textrm{D}_{\varphi } \tilde{\mathscr {F}}(\varphi ,\tilde{\mathcal {M}}) \cdot \delta \tilde{\varphi }&= \int _{\mathcal {B}}\delta _{\hat{\mathscr {F}}}{\hat{g}} \ \dot{\wedge }\ \delta \hat{g}+ \delta _{\hat{\mathscr {F}}}{\hat{\mathcal {M}}} \ \dot{\wedge }\ \delta \hat{\mathcal {M}}_1\\&= \int _{\mathcal {B}}\delta _{\hat{\mathscr {F}}}{\hat{g}} \ \dot{\wedge }\ \mathcal {L}_{\delta \hat{\varphi }}{\hat{g}} + \delta _{\hat{\mathscr {F}}}{\hat{\mathcal {M}}} \ \dot{\wedge }\ \hat{\mu }\otimes (\hat{\nabla }\delta \hat{\varphi }\ \dot{\wedge }\ \hat{v}^\flat )\\&= \int _{\mathcal {B}}\delta _{\hat{\mathscr {F}}}{\hat{g}} \ \dot{\wedge }\ 2 \textrm{sym}(\hat{\textrm{d}}_{\hat{\nabla }}\hat{v}^\flat ) + \hat{\nabla }\delta \hat{\varphi }\ \dot{\wedge }\ \iota _{\delta _{\hat{\mathscr {F}}}{\hat{\mathcal {M}}}} \hat{\mathcal {M}}= \int _{\mathcal {B}}\hat{\nabla }\delta \hat{\varphi }\ \dot{\wedge }\ \mathcal {E}_{\hat{F}}, \end{aligned}$$

with $\mathcal {E}_{\hat{\mathscr {F}}}:= (2 (\delta _{\hat{g}}{\hat{\mathscr {F}}}) ^\flat + \iota _{\delta _{\hat{\mathcal {M}}}{\hat{\mathscr {F}}}}\hat{\mathcal {M}}\in \varOmega ^{n-1}(\mathcal {B};T^*\mathcal {B})$. Using integration by parts, the definition of $\delta \hat{\varphi }$ and the duality of $\varphi _{\textrm{v},*}$ and $\varphi _{\textrm{v}}^*$, one has that

$$\begin{aligned} \int _{\mathcal {B}}\hat{\nabla }\delta \hat{\varphi }\ \dot{\wedge }\ \mathcal {E}_{\hat{F}} =&- \int _{\mathcal {B}}\delta \hat{\varphi }\ \dot{\wedge }\ \hat{\textrm{d}}_{\hat{\nabla }}\mathcal {E}_{\hat{F}} + \int _{\partial \mathcal {B}} \delta \hat{\varphi }|_{{\partial \mathcal {B}}} \ \dot{\wedge }\ \mathcal {E}_{\hat{F}}|_{{\partial \mathcal {B}}}\\ =&- \int _{\mathcal {B}}\delta \tilde{\varphi }\ \dot{\wedge }\ \varphi _{\textrm{v},*}(\hat{\textrm{d}}_{\hat{\nabla }}\mathcal {E}_{\hat{F}}) + \int _{\partial \mathcal {B}} \delta \tilde{\varphi }|_{{\partial \mathcal {B}}} \ \dot{\wedge }\ \varphi _{\textrm{v},*}(\mathcal {E}_{\hat{F}})|_{{\partial \mathcal {B}}}. \end{aligned}$$

Finally, the results follow by comparing with (19).

1.3 Vector Space Decomposition

Let V, W be two vector spaces over the field $\mathbb {R}$. Let $\text {Hom}(V,W)$ denote te set of linear maps from V to W. An endomorphism of V is an element of $\text {End}(V):= \text {Hom}(V,V)$. The dual vector space to V is $V^*:= \text {Hom}(V,\mathbb {R})$ with the duality product between any $v\in V,\tilde{v}\in V^*$ denoted by $\langle \tilde{v} | v \rangle _{V}:= \tilde{v}(v) \in \mathbb {R}$. The orthogonal complement of a subspace $W\subset V$ is the subspace $W^\perp \subset V^*$ defined as

$$W^\perp := \{\tilde{w}\in V^* | \langle \tilde{w} | w \rangle _{V} = 0,\ \forall w\in W\}.$$

Let M be a smooth manifold of dimension n, and let $T_p M,T_p^*M$ denote the tangent and cotangent spaces of M at a point $p\in M$, respectively. The space of all (r, s) tensor at $p\in M$ is denoted by $T^r_{s,p}M$ with $T^0_{0,p}M = \mathbb {R}, T^1_{0,p}M =T_p M,$ and $T^0_{1,p}M = T_p^*M$. Furthermore, $T^0_{2,p}M = \text {Hom}(T_p M,T_p^*M)$ and $T^1_{1,p}M = \text {End}(T_p M)$.

Given a choice of basis for $T_p M$ and its dual basis for $T_p^*M$, one has that

$$T^0_{2,p}M \cong \mathfrak {bl}_n, \qquad \qquad T^1_{1,p}M \cong \mathfrak {gl}_n,$$

where $\mathfrak {bl}_n$ and $\mathfrak {gl}_n$ denote the vector spaces of bilinear forms and linear transformations on $\mathbb {R}^n$, respectively. Elements of $\mathfrak {bl}_n$ and $\mathfrak {gl}_n$ are expressed in components as

$$A = A_{ij} e^i\otimes e^j \in \mathfrak {bl}_n, \qquad \qquad X = X^i_j e_i\otimes e^j \in \mathfrak {gl}_n,$$

where $\{e_i\}, \{e^j\}$ denote the standard basis and its dual on $\mathbb {R}^n$. Furthermore, elements of the dual spaces $\mathfrak {bl}_n^*$ and $\mathfrak {gl}_n^*$ are expressed in components as

$$B = B^{ij} e_i\otimes e_j \in \mathfrak {bl}_n^*, \qquad \qquad Y = Y^i_j e^i\otimes e_j\in \mathfrak {gl}_n^*.$$

The duality product ${\langle \cdot | \cdot \rangle _{\mathfrak {bl}_n}}:{\mathfrak {bl}_n^*\times \mathfrak {bl}_n}\rightarrow {\mathbb {R}}$ is defined by $\langle B | A \rangle _{\mathfrak {bl}_n}:= B^{ij}A_{ij},$ whereas the duality product ${\langle \cdot | \cdot \rangle _{\mathfrak {gl}_n}}:{\mathfrak {gl}_n^*\times \mathfrak {gl}_n}\rightarrow {\mathbb {R}}$ is defined by $\langle Y | X \rangle _{\mathfrak {bl}_n}:= Y_i^j X^i_j.$ If $\varvec{A},\varvec{B} \in \mathbb {R}^{n\times n}$ denote the matrix representations of $A\in \mathfrak {bl}_n$ and $B\in \mathfrak {bl}_n^*$, respectively, then the duality products is more commonly expressed as $\langle B | A \rangle _{\mathfrak {bl}_n} = \textrm{tr}(\varvec{B}\varvec{A})$. The same expression holds also for $\mathfrak {gl}_n$.

1.3.1 Volumetric–Deviatoric Decomposition of $\mathfrak {gl}_n$

Consider the following projection map

$$\begin{aligned} {{\pi }_{\textrm{vol}}}:{\mathfrak {gl}_n}\rightarrow {\mathfrak {gl}_n}, \qquad \qquad {X}\mapsto {\frac{1}{n} \textrm{tr}(X)I_n,} \end{aligned}$$

(98)

where ${\textrm{tr}}:{\mathfrak {gl}_n}\rightarrow {\mathbb {R}}$ denotes the trace map and $I_n \in \mathfrak {gl}_n$ denotes the identity element (with components $\delta _j^i$). It is straightforward to see that ${\pi }_{\textrm{vol}}$ is idempotent (i.e. ${\pi }_{\textrm{vol}}\circ {\pi }_{\textrm{vol}}= {\pi }_{\textrm{vol}}$) which necessitates the division by n. Let $\mathfrak {k}:= \text {im}({\pi }_{\textrm{vol}}) \subset \mathfrak {gl}_n$ and $\mathfrak {sl}_n:= \text {ker}({\pi }_{\textrm{vol}}) \subset \mathfrak {gl}_n$ denote the vector subspaces of scalar operators and traceless operators, respectively. Any $X\in \mathfrak {gl}_n$ can be decomposed then as $X = {X}_{\textrm{vol}}+ {X}_{\textrm{dev}}$ with

$$ {X}_{\textrm{vol}}:= {\pi }_{\textrm{vol}}(X) \in \mathfrak {k}, \qquad {X}_{\textrm{dev}}:= {\pi }_{\textrm{dev}}(X):= X -{\pi }_{\textrm{vol}}(X) \in \mathfrak {sl}_n,$$

where ${X}_{\textrm{vol}}$ and ${X}_{\textrm{dev}}$ will be referred to as the volumetric and deviatoric components of X, respectively. By construction, one has that ${\pi }_{\textrm{vol}}({\pi }_{\textrm{dev}}(X)) = 0$. We refer to the decomposition $\mathfrak {gl}_n= \mathfrak {k}\oplus \mathfrak {sl}_n,$ as the volumetric–deviatoric decomposition of $\mathfrak {gl}_n$. While $\text {dim}(\mathfrak {gl}_n) = n^2$, one has that $\text {dim}(\mathfrak {k}) = 1$ and $\text {dim}(\mathfrak {sl}_n) = n^2 -1$.

The dual space of $\mathfrak {gl}_n$ is decomposed as $\mathfrak {gl}_n^*= \mathfrak {k}^*\oplus \mathfrak {sl}_n^*$ where

$$\begin{aligned} \mathfrak {k}^*&= \mathfrak {sl}_n^\perp := \{Y \in \mathfrak {gl}_n^*| \langle Y | {X}_{\textrm{dev}} \rangle _{\mathfrak {gl}_n} = 0,\ \forall {X}_{\textrm{dev}}\in \mathfrak {sl}_n\},\\ \mathfrak {sl}_n^*&= \mathfrak {k}^\perp := \{Y \in \mathfrak {gl}_n^*| \langle Y | {X}_{\textrm{vol}} \rangle _{\mathfrak {gl}_n} = 0,\ \forall {X}_{\textrm{vol}}\in \mathfrak {k}\}. \end{aligned}$$

If ${{\pi }_{\textrm{vol}}}:{\mathfrak {gl}_n}\rightarrow {\mathfrak {k}}$ and ${{\pi }_{\textrm{dev}}}:{\mathfrak {gl}_n}\rightarrow {\mathfrak {sl}_n}$ are defined as above, then ${{\tilde{\pi }}_{\textrm{vol}}}:{\mathfrak {gl}_n^*}\rightarrow {\mathfrak {k}^*}$ and ${{\tilde{\pi }}_{\textrm{dev}}}:{\mathfrak {gl}_n^*}\rightarrow {\mathfrak {sl}_n^*}$ are defined such that $\langle {\tilde{\pi }}_{\textrm{vol}}(Y) | {\pi }_{\textrm{dev}}(X) \rangle _{\mathfrak {gl}_n} = 0$ and $\langle {\tilde{\pi }}_{\textrm{dev}}(Y) | {\pi }_{\textrm{vol}}(X) \rangle _{\mathfrak {bl}_n} = 0.$ One can show that ${\tilde{\pi }}_{\textrm{vol}}(Y) = \frac{1}{n} \textrm{tr}(X)I_n=: {Y}_{\textrm{vol}}$ and ${\tilde{\pi }}_{\textrm{dev}}(Y) = Y -{Y}_{\textrm{vol}}=: {Y}_{\textrm{dev}}.$ Consequently, one has that the duality pairing between any $X\in \mathfrak {gl}_n$ and $Y\in \mathfrak {gl}_n^*$ to be decomposed as

$$\langle Y | X \rangle _{\mathfrak {gl}_n} = \langle {Y}_{\textrm{vol}} | {X}_{\textrm{vol}} \rangle _{\mathfrak {gl}_n} + \langle {Y}_{\textrm{dev}} | {X}_{\textrm{dev}} \rangle _{\mathfrak {gl}_n}.$$

Since the dimension of the subspace $\mathfrak {k}$ is 1, it is isomorphic to $\mathbb {R}$. However, this isomorphism is not canonical. In this work, we choose to identify $\mathfrak {k}$ with $\mathbb {R}$ with the following isomorphism

$$\begin{aligned} {\textrm{tr}}:{\mathfrak {k}}\rightarrow {\mathbb {R}},\qquad \qquad {{X}_{\textrm{vol}}}\mapsto {\textrm{tr}({X}_{\textrm{vol}})=: {x}_{\textrm{vol}}.} \end{aligned}$$

The dual map ${\textrm{tr}^*}:{\mathbb {R}}\rightarrow {\mathfrak {k}^*}$ is given by $\textrm{tr}^*({y}_{\textrm{vol}}) = {y}_{\textrm{vol}} I_n$. Consequently, one can express the duality pairing on $\mathfrak {k}$ as a product of scalars such that

$$\langle {Y}_{\textrm{vol}} | {X}_{\textrm{vol}} \rangle _{\mathfrak {gl}_n} = {y}_{\textrm{vol}} {x}_{\textrm{vol}}.$$

One can see that ${x}_{\textrm{vol}} = \textrm{tr}({X}_{\textrm{vol}}) = \textrm{tr}(X)$ while ${y}_{\textrm{vol}} = \frac{1}{n} \textrm{tr}(Y)$. In principle, one can modify the above isomorphism between $\mathfrak {k}$ and $\mathbb {R}$ such that the factor $\frac{1}{n}$ is multiplicatively distributed between ${x}_{\textrm{vol}}$ and ${y}_{\textrm{vol}}$ in any way.

The above construction of the volumetric–deviatoric decomposition of $\mathfrak {gl}_n$ can be summarized in the following map ${J_{v/d}}:{\mathfrak {k}^*\times \mathfrak {sl}_n^*\times \mathfrak {gl}_n}\rightarrow {\mathfrak {k}\times \mathfrak {sl}_n\times \mathfrak {gl}_n^*}$

$$\begin{aligned} \begin{pmatrix} {x}_{\textrm{vol}}\\ {X}_{\textrm{dev}}\\ Y \end{pmatrix} = \begin{pmatrix} 0 & 0 & \textrm{tr}\circ {\pi }_{\textrm{vol}}\\ 0 & 0 & {\pi }_{\textrm{dev}}\\ \textrm{tr}^* & \mathbbm {1} & 0 \end{pmatrix} \begin{pmatrix} {y}_{\textrm{vol}}\\ {Y}_{\textrm{dev}}\\ X \end{pmatrix}, \end{aligned}$$

such that the duality pairing between $X\in \mathfrak {gl}_n$ and $Y\in \mathfrak {gl}_n^*$ is given by

$$\langle Y | X \rangle _{\mathfrak {gl}_n} = {y}_{\textrm{vol}} {x}_{\textrm{vol}} + \langle {Y}_{\textrm{dev}} | {X}_{\textrm{dev}} \rangle _{\mathfrak {gl}_n}.$$

1.3.2 Decomposition of Second-Rank Tensor Fields on M

One can easily apply the same line of thought in decomposing $\mathfrak {gl}_n$ above to $\mathfrak {bl}_n$ by decomposing it into symmetric and asymmetric subspaces. Consequently, we can extend the volumetric–deviatoric decomposition of $\mathfrak {gl}_n\cong T^1_{1,p}M$ to the space of (1,1) tensor fields $\varGamma (T^1_1M)$ and the symmetric–asymmetric decomposition of $\mathfrak {bl}_n\cong T^0_{2,p}M$ to the space of (0,2) tensor fields $\varGamma (T^0_2M)$ such that

$$\varGamma (T^0_2M) = \varGamma _{\text {s}}(T^0_2M) \oplus \varGamma _{\text {a}}(T^0_2M), \qquad \varGamma (T^1_1M) = \varGamma _{\text {v}}(T^1_1M) \oplus \varGamma _{\text {d}}(T^1_1M).$$

The duality product between $A\in \varGamma (T^0_2M)$ and $B\in \varGamma (T^2_0M)$ can then be written as

$$\left\langle \left\langle B | A \right\rangle \right\rangle _{M} = \left\langle \left\langle {B}_{\textrm{sym}} | {A}_{\textrm{sym}} \right\rangle \right\rangle _{M} + \left\langle \left\langle {B}_{\textrm{asy}} | {A}_{\textrm{asy}} \right\rangle \right\rangle _{M},$$

with $\left\langle \left\langle B | A \right\rangle \right\rangle _{M}:= \int _M B:A\ \omega$, such that $B:A\in C^\infty (M)$ denotes the double contraction of the tensor fields A and B and $\omega$ denotes the volume form of M. Similarly, the duality product between $X\in \varGamma (T^1_1M)$ and $Y\in \varGamma (T^1_1M)$ can then be written as

$$\left\langle \left\langle X | Y \right\rangle \right\rangle _{M} = \langle {y}_{\textrm{vol}} | {x}_{\textrm{vol}} \rangle _{M} + \left\langle \left\langle {Y}_{\textrm{dev}} | {X}_{\textrm{dev}} \right\rangle \right\rangle _{M},$$

while ${\langle \cdot | \cdot \rangle _{M} }:{C^\infty (M)\times C^\infty (M)}\rightarrow {\mathbb {R}}$ is defined as $\langle {y}_{\textrm{vol}} | {x}_{\textrm{vol}} \rangle _{M}:= \int _M {y}_{\textrm{vol}} {x}_{\textrm{vol}}\ \omega$.

It is important to note that the above decomposition of second rank tensor fields did not require the manifold M to have a Riemannian g. In the presence of this extra metric structure, then one could apply the volumetric–deviatoric decomposition to any $A\in \varGamma (T^0_2M)$ using its tensorial variants by raising and lowering indices. This is achieved by extending the projection map (98) to

$${{\bar{\pi }}_{\textrm{vol}}:= g \circ {\pi }_{\textrm{vol}}\circ g^{-1}}:{\varGamma (T^0_2M)}\rightarrow {\varGamma (T^0_2M)}.$$

In the same manner, one could apply the symmetric–asymmetric decomposition to any $X\in \varGamma (T^1_1M)$ using

$${{\bar{\pi }}_{\textrm{sym}}:= g^{-1} \circ {\pi }_{\textrm{sym}}\circ g}:{\varGamma (T^1_1M)}\rightarrow {\varGamma (T^1_1M)}.$$