11institutetext:  I. Steinbrecher, N. Hagmeyer, A. Popp22institutetext:
Institute for Mathematics and Computer-Based Simulation,
University of the Bundeswehr Munich,
Werner-Heisenberg-Weg 39, D-85577 Neubiberg, Germany
22email<EMAIL_ADDRESS>33institutetext: C. Meier44institutetext:
Institute for Computational Mechanics,
Technical University of Munich,
Boltzmannstrasse 15, D-85748 Garching b. München, Germany

# A consistent mixed-dimensional coupling approach for 1D Cosserat beams and
2D solid surfaces

Ivo Steinbrecher  Nora Hagmeyer  Christoph Meier  Alexander Popp

(Received: date / Accepted: date)

###### Abstract

The present article proposes a novel computational method for coupling 1D
fibers with the 2D surface of a 3D solid. The fibers are modeled as 1D
Cosserat continua (beams) with six local degrees of freedom, three positional
and three rotational ones. A kinematically consistent 1D-2D coupling scheme
for this problem type is proposed considering the positional and rotational
degrees of freedom along the beams. The positional degrees of freedom are
coupled by enforcing a constant normal distance between a point on the beam
centerline and a corresponding point on the solid surface. This strategy
requires a consistent description of the solid surface normal vector field to
guarantee fundamental mechanical properties such as conservation of angular
momentum. Coupling of the rotational degrees of freedom of the beams and a
suitable rotation tensor representative for the local orientation within a
Boltzmann continuum has been considered in a previous contribution. In the
present work, this coupling approach will be extended by constructing rotation
tensors that are representative for local surface orientations. Several
numerical examples demonstrate the consistency, robustness and accuracy of the
proposed method. To showcase its applicability to multi-physics systems of
practical relevance, the fluid-structure interaction example of a vascular
stent is presented.

###### Keywords:

Beam-to-surface coupling 1D-2D position and rotation coupling Mixed-
dimensional coupling Finite element method Geometrically exact beam theory
Mortar methods

††journal: Computational Mechanics

##  1 Introduction

Compound structures consisting of one-dimensional (1D) beam-like structures,
i.e., structures with one spatial dimension being much larger than the other
two, tied to three-dimensional (3D) solid continua can be found in a variety
of different fields. From a geometrical point of view, the problem considered
in the present work consists of 1D beams coupled to the two-dimensional (2D)
surface of a 3D continuum (solid). This will be denoted as beam-to-solid
surface (BTSS) coupling problem throughout this contribution. Applications for
this class of problems can be found in, e.g., civil engineering, where steel
girders are used to support concrete slabs, or in mechanical engineering,
where lightweight structures are realized by stabilizing thin shells with
struts. Leaving the realm of classical engineering applications, BTSS coupling
approaches can also be used to model biomechanical systems, for example to
capture the interaction between stent and graft (encasing fabric) as used for
endovascular aneurysm repair. Numerical simulation of such applications is of
high importance during the development and design phase to accurately predict
and control the desired system behavior.

The present modeling approach for the BTSS coupling problem employs accurate
and efficient 1D models for the beam-like structures based on geometrically
exact beam theory Reissner1972 ; Simo1986a ; Simo1986b ; Cardona1988 ;
Ibrahimbegovic1995 ; Crisfield1999 ; Romero2004 ; Meier2019 ; Betsch2002 ;
Sonneville2014 . The beams are represented by 1D curves in 3D space, i.e., the
beam centerline that connects the centroids of the beam cross-sections. Each
point along the beam centerline has six degrees of freedom (three positional
and three rotational ones), i.e., the beam model can be identified as a 1D
Cosserat continuum. The solid body is modeled as a classical 3D Boltzmann
continuum. The 1D beams are coupled to the 2D surfaces of the 3D continua,
thus resulting in a mixed-dimensional 1D-2D coupling problem. The resulting
BTSS coupling problem has two desirable features, cf. Steinbrecher2020 ;
Steinbrecher2022 : (i) Both the solid and the beam can be modeled and
discretized individually. Therefore, well-established discretization schemes
for the solid and the beam can be used without modifications. (ii) Employing
1D beam models results in computationally efficient finite element
discretizations, which reduces the number of unknowns required to represent
the beam-like structures by several orders of magnitude as compared to a
modeling approach based on 3D continuum theory. In the literature, mixed-
dimensional coupling between structural beam theories and solid continua is
often addressed to model fiber-reinforced materials, e.g., Phillips1976 ;
Chang1987 ; Elwi1989 ; Ranjbaran1996 ; Gomes2001 ; Kerfriden2020 . However, in
all of the aforementioned works, the coupling is an embedded 1D-3D coupling
since the 1D fiber reinforcements are placed inside the solid domain and are
directly coupled to the 3D solid volume. Furthermore, 1D string-like models
with a limited representation of the relevant modes of deformation, i.e., only
axial deformation, were used in these contributions to represent the 1D
curves. Coupling approaches for full beam theories in 1D-3D beam-to-solid
volume (BTSV) coupling problems have been developed more recently, e.g., in
Durville2007 ; Steinbrecher2020 ; Steinbrecher2022 ; Khristenko2021 . Compared
to the previously mentioned string models, beam theories contain additional
deformations modes, i.e., bending, torsion and shear, which allows for a more
realistic representation of the nonlinear force-displacement relations caused
by the reinforcements. In Durville2007  a collocation method is used to couple
the beams to the solid. A mortar-type approach to couple the positional
degrees of freedom of the beam centerline to the solid is presented in
Steinbrecher2020  and a mortar-type approach for _full_ coupling, i.e.,
positional _and_ rotational coupling, is presented in Steinbrecher2022 . In
Khristenko2021  the coupling constraints of a BTSV problem are formulated on
the surface of the beam and are subsequently projected onto the beam
centerline considering a Taylor series expansion of the solid displacement
field. Apart from the 1D-3D coupling problems discussed so far, a truly 1D-2D
coupling is presented in Konyukhov2015 . However, the solid surfaces are
assumed to be rigid in that contribution, which heavily limits the
applicability to real life engineering problems.

In the present work, we propose the first truly mixed-dimensional 1D-2D
mortar-type approach for BTSS coupling problems. This is an extension of the
authors’ previous contributions on positional beam-to-solid volume (BTSV-POS)
and rotational beam-to-solid volume (BTSV-ROT) coupling problems
Steinbrecher2020 ; Steinbrecher2022  to BTSS coupling problems. The transition
from a 1D-3D to a 1D-2D mixed-dimensional coupling introduces two additional
challenges: (i) The positional BTSS coupling (BTSS-POS) constraints between
the beam and the solid surface depend on the surface normal vector. In this
work, the term _consistent_ implies that no simplifications regarding the
surface normal vector are introduced in the further derivation of the
positional coupling constraints. A consistent treatment of the surface normal
vector, especially in the discretized problem setting can become cumbersome.
Therefore, different possible simplifications of the consistent positional
coupling constraints (BTSS-POS-$(\cdot)$) are presented. As a main scientific
contribution of this work, it is demonstrated that in the general case of non-
matching 1D-2D interfaces only a fully consistent handling of the surface
normal vector within the coupling constraints allows to fulfill fundamental
mechanical properties and to give accurate results. In particular, exact
conservation of linear and angular momentum is shown for the resulting 1D-2D
coupling scheme. To the authors’ best knowledge, this is the first time that
exact conservation of angular momentum is shown for a surface coupling scheme
with non-vanishing surface normal distance. (ii) For rotational BTSS (BTSS-
ROT) coupling a suitable solid orientation field is required on the solid
surface. A detailed discussion on rotation tensors that are representative for
the local orientation of a solid continuum is given in Steinbrecher2022 .
However, a direct application of these approaches to solid surfaces leads to
undesirable effects, i.e., the solid surface orientation would depend on the
deformations inside the solid volume. Therefore, the second main scientific
contribution of this work is the construction of a suitable solid surface
orientation field. In the remainder of this work, BTSS-FULL refers to
positional _and_ rotational BTSS coupling.

Within this contribution the solid surface is exclusively considered as
boundary of a three-dimensional solid volume. For standard Lagrangian finite
element interpolations of the solid, the solid surface normal field is non-
continuous across element edges. To recover $C^{0}$-continuity, a re-
interpolation scheme for the surface normal vectors is employed, comparable to
the evaluation of the surface normal field in surface-to-surface contact
schemes, e.g., Yang2005 ; Popp2009 .

Eventually, it is emphasized that a modeling approach based on mixed-
dimensional coupling influences the nature of the underlying mechanical
problem. In the context of embedded 1D-3D coupling this issue has been
thoroughly discussed and analyzed for the cases of positional coupling
Steinbrecher2020  and rotational coupling Steinbrecher2022 . One of the main
consequences is that the analytical solution of the mixed-dimensional coupling
problem exhibits a local singularity at the position of the beam centerline.
In the embedded 1D-3D positional coupling case this can be interpreted as a
generalization of the well-known Kelvin problem PodioGuidugli2014 ; Kelvin1848
, cf. Figure 1, i.e., a line load acting on an infinite solid. The same issue
arises in the considered case of 1D-2D BTSS coupling, which corresponds to the
Flamant problem of a line load acting on an infinite half space
PodioGuidugli2005 , cf. Figure 1. However, similar to the BTSV case this
aspect does not impact the applicability of the proposed BTSS coupling method
for the envisioned range of practically relevant discretization resolutions,
i.e., solid element sizes in the range of the beam cross-section diameter or
above. As discussed in detail in Steinbrecher2020 ; Steinbrecher2022 , the
aforementioned singularity does not occur for this range of solid mesh sizes.

Figure 1: 1 The Kelvin problem of an embedded line load acting on an infinite
solid and 1 the Flamant problem of a line load acting on an infinite solid
half space.

The remainder of this work is organized as follows: In Section 2, we state the
governing equations for solid and beam formulations as well as for the BTSS-
FULL method. In LABEL:sec:surface_triad, a suitable procedure for constructing
the solid surface triad is presented to couple the rotations of the solid
surface to the beam cross-section orientations. The finite element
discretization of the BTSS-FULL method is presented in
LABEL:sec:discretization. Furthermore, the construction of a
$C^{0}$-continuous surface normal field is elaborated. Finally, in
LABEL:sec:examples, we present numerical examples to demonstrate the
consistency of the presented BTSS-FULL method and the applicability to real-
life engineering and biomedical problems. Detailed comparisons with full 3D
continuum approaches, i.e., beam and solid are modeled as 3D continua and
discretized by 3D solid finite elements, for the BTSS coupling problem are
presented. Furthermore, the importance of coupling both positions _and_
rotations for beam-to-solid surface coupling problems is shown.

##  2 Problem formulation

We consider a quasi-static 3D finite deformation BTSS-FULL coupling problem as
shown in Figure 2. It is emphasized that the presented BTSS-FULL method is not
restricted to quasi-static problems, but can directly be applied to time-
dependent problems as well. A Cartesian frame $\underline{\bm{e}}_{1}$,
$\underline{\bm{e}}_{2}$ and $\underline{\bm{e}}_{3}$ serves as fixed frame of
reference. The principle of virtual work (PVW) serves as basis for a
subsequent finite element discretization and reads

$\delta W^{S}+\delta W^{B}+\delta\Pi_{\lambda}=0.$ (1)

Here, $\delta W^{S}$ is the total virtual work of the pure solid problem,
$\delta W^{B}$ is the total virtual work of the pure beam problem and
$\delta\Pi_{\lambda}$ is the virtual work due to coupling forces and moments.
Contributions to the total virtual work of the pure solid and beam problem are
independent of the coupling constraints, i.e., well-established modeling and
discretization techniques can be used for the solid and the beam, cf.
Steinbrecher2020 ; Steinbrecher2022 .

Figure 2: Notation for the finite deformation BTSS-FULL coupling problem.

###  2.1 Finite rotations

Before stating the governing equations for the solid, beam and BTSS-FULL
coupling problem, a short recap on finite rotations is given here, as a
consistent treatment of large rotations is required for the rotational
coupling conditions BTSS-ROT. In geometrically exact beam theory, the term
triad is commonly used to describe the set of three orthonormal vectors
defining a beam cross-section orientation, i.e.,

$\underline{\bm{\Lambda}}=\begin{bmatrix}\underline{\bm{g}}_{1},\underline{\bm{g}}_{2},\underline{\bm{g}}_{3}\end{bmatrix}\in
SO^{3}.$ (2)

Here, $SO^{3}$ is the special orthogonal group and $\underline{\bm{g}}_{i}$
are the base vectors of the triad. The triad is equivalent to a rotation
tensor, mapping the Cartesian basis vectors $\underline{\bm{e}}_{i}$ onto
$\underline{\bm{g}}_{i}$. Among others, a triad can be parameterized by the
rotation (pseudo-)vector $\underline{\bm{\psi}}$, i.e.,
$\underline{\bm{\Lambda}}=\underline{\bm{\Lambda}}(\underline{\bm{\psi}})$.
The rotation vector describes a rotation by an angle
$\psi=\left\|\underline{\bm{\psi}}\right\|$ around the rotation axis
$\underline{\bm{e}}_{\psi}=\underline{\bm{\psi}}/\left\|\underline{\bm{\psi}}\right\|$.
The parametrization can be evaluated with the well-known Rodrigues formula
Argyris1982

$\begin{split}\underline{\bm{\Lambda}}(\underline{\bm{\psi}})&=\exp(\underline{\bm{S}}(\underline{\bm{\psi}}))\\\
&=\underline{\bm{I}}+\sin\psi\underline{\bm{S}}(\underline{\bm{e}}_{\psi})+(1-\cos\psi)\underline{\bm{S}}()^{2}(\underline{\bm{e}}_{\psi}),\end{split}$
(3)

where $\exp(\cdot)$ is the exponential map and $\underline{\bm{S}}()$ is a
skew-symmetric tensor such that
$\underline{\bm{S}}(\underline{\bm{a}})\underline{\bm{b}}=\underline{\bm{a}}\times\underline{\bm{b}}\
\forall\ \underline{\bm{a}},\underline{\bm{b}}\in\mathbb{R}^{3}$. The
calculation of the inverse of the Rodrigues formula is not straight forward.
For simplicity, it is abbreviated by the expression
$\underline{\bm{\psi}}(\underline{\bm{\Lambda}})=\operatorname{rv}(\underline{\bm{\Lambda}})$
in the following. In practice, Spurrier’s algorithm Spurrier1978  can be used
for the extraction of the rotation vector. Let us consider two triads
$\underline{\bm{\Lambda}}_{1}(\underline{\bm{\psi}}_{1})$ and
$\underline{\bm{\Lambda}}_{2}(\underline{\bm{\psi}}_{2})$ with their
respective rotation vectors $\underline{\bm{\psi}}_{1}$ and
$\underline{\bm{\psi}}_{2}$. They are related to each other by the relative
rotation $\underline{\bm{\Lambda}}_{21}(\underline{\bm{\psi}}_{21})$. The
relative rotation is given by

$\begin{split}\underline{\bm{\Lambda}}_{2}(\underline{\bm{\psi}}_{2})&=\underline{\bm{\Lambda}}_{21}(\underline{\bm{\psi}}_{21})\underline{\bm{\Lambda}}_{1}(\underline{\bm{\psi}}_{1})\\\
&\Updownarrow\\\
\underline{\bm{\Lambda}}_{21}(\underline{\bm{\psi}}_{21})&=\underline{\bm{\Lambda}}_{2}(\underline{\bm{\psi}}_{2})\underline{\bm{\Lambda}}_{1}(\underline{\bm{\psi}}_{1})^{\mathrm{T}},\end{split}$
(4)

with the identity
$\underline{\bm{\Lambda}}^{\mathrm{T}}=\underline{\bm{\Lambda}}^{-1}$ for all
elements of $SO^{3}$. The relative rotation vector
$\underline{\bm{\psi}}_{21}=\operatorname{rv}(\underline{\bm{\Lambda}}_{21})$
describes the relative rotation between $\underline{\bm{\Lambda}}_{1}$ and
$\underline{\bm{\Lambda}}_{2}$. Rotation vectors are non-additive, i.e.,
$\underline{\bm{\psi}}_{21}\neq\underline{\bm{\psi}}_{2}-\underline{\bm{\psi}}_{1}$.
For a more comprehensive treatment of this topic, the interested reader is
referred to Simo1986b ; Cardona1988 ; Ibrahimbegovic1995 ; Romero2004 ;
Meier2019 ; Betsch1998 . In the following sections, both symbols
$\underline{\bm{\Lambda}}$ and $\underline{\bm{R}}$ will be used to represent
rotation tensors.

###  2.2 Solid formulation

The solid is modeled as a 3D Boltzmann continuum. The solid domain in the
reference configuration is $\Omega_{S,0}\subset\mathbb{R}^{3}$ and
$\partial\Omega_{S,0}$ is the boundary of the solid domain, i.e., the solid
surface. Throughout this work, the subscript $(\cdot)_{0}$ indicates a
quantity in the reference configuration. Accordingly, $\Omega_{S}$ and
$\partial\Omega_{S}$ are the current solid domain and current solid surface,
respectively. Furthermore, the current position
$\underline{\bm{x}}_{S}\in\mathbb{R}^{3}$ of a solid material point relates to
the reference position $\underline{\bm{X}}_{S}\in\mathbb{R}^{3}$ via the solid
displacement field $\underline{\bm{u}}_{S}\in\mathbb{R}^{3}$, i.e.,

$\underline{\bm{x}}_{S}(\underline{\bm{X}}_{S})=\underline{\bm{X}}_{S}+\underline{\bm{u}}_{S}(\underline{\bm{X}}_{S}).$
(5)

The virtual work contributions of the solid domain read

$\begin{split}\delta
W^{S}=&\int_{\Omega_{S,0}}{\underline{\bm{S}}:\delta\underline{\bm{E}}\,\mathrm{d}V_{0}}\,\\\
&-\int_{\Omega_{S,0}}{\hat{\underline{\bm{b}}}\cdot\delta\underline{\bm{u}}_{S}\,\mathrm{d}V_{0}}\,-\int_{\Gamma_{\sigma}}{\hat{\underline{\bm{t}}}\cdot\delta\underline{\bm{u}}_{S}\,\mathrm{d}A_{0}}\,.\end{split}$
(6)

Here, $\delta$ denotes the (total) variation of a quantity,
$\underline{\bm{S}}\in\mathbb{R}^{3\times 3}$ is the second Piola–Kirchhoff
stress tensor,
$\underline{\bm{E}}=\frac{1}{2}(\underline{\bm{F}}^{\mathrm{T}}\underline{\bm{F}}-\underline{\bm{I}})\in\mathbb{R}^{3\times
3}$ is the work-conjugated Green–Lagrange strain tensor,
$\hat{\underline{\bm{b}}}\in\mathbb{R}^{3}$ is the body load vector and
$\hat{\underline{\bm{t}}}\in\mathbb{R}^{3}$ are surface tractions on the
Neumann boundary $\Gamma_{\sigma}\subset\partial\Omega_{S,0}$. Furthermore,
$\underline{\bm{F}}\in\mathbb{R}^{3\times 3}$ is the solid deformation
gradient, which is defined according to

$\underline{\bm{F}}=\frac{\partial\underline{\bm{x}}_{S}}{\partial\underline{\bm{X}}_{S}}.$
(7)

If e.g., a hyperelastic strain energy function $\Psi(\underline{\bm{E}})$
exists, the constitutive relations between stresses and strains can be stated
as
$\underline{\bm{S}}=\frac{\partial\Psi(\underline{\bm{E}})}{\partial\underline{\bm{E}}}$.

###  2.3 Geometrically exact beam theory

In this work the geometrically exact Simo–Reissner beam theory is used to
describe the embedded beams as 1D Cosserat continua, e.g., Reissner1972 ;
Meier2019 ; Simo1986a ; Simo1986b . Each beam cross-section along the beam
centerline is described by six degrees of freedom, three positional and three
rotational ones, thus resulting in six deformation modes of the beam: axial
tension, bending (2$\times$), shear (2$\times$) and torsion.

The complete beam kinematics can be defined by a centerline curve
$\underline{\bm{r}}(s)\in\mathbb{R}^{3}$, connecting the cross-section
centroids, and a field of triads
$\underline{\bm{\Lambda}}_{B}(s)=\underline{\bm{\Lambda}}_{B}(\underline{\bm{\psi}}_{B}(s))$
defining the orientation of the beam cross-sections. Here,
$s\in[0,L]=:\Omega_{B,0}$ is the arc-length coordinate along the beam
centerline in the reference configuration and $L$ is the reference length of
the beam. The triad $\underline{\bm{\Lambda}}_{B}$ is chosen such that the
second and third basis vectors, $\underline{\bm{g}}_{B2}$ and
$\underline{\bm{g}}_{B3}$, span the beam cross-section, i.e., the first triad
basis vector $\underline{\bm{g}}_{B1}$ is normal to the beam cross-section. A
total hyperelastic stored-energy function of the Simo–Reissner beam can be
stated as

$\Pi_{\text{int},B}=\int_{\Omega_{B,0}}{\tilde{\Pi}_{\text{int},B}\,\mathrm{d}s}\,,$
(8)

with

$\tilde{\Pi}_{\text{int},B}=\frac{1}{2}(\underline{\bm{\Gamma}}^{\mathrm{T}}\underline{\bm{C}}_{F}\underline{\bm{\Gamma}}+\underline{\bm{\Omega}}^{\mathrm{T}}\underline{\bm{C}}_{M}\underline{\bm{\Omega}}).$
(9)

Here, $\underline{\bm{\Gamma}}\in\mathbb{R}^{3}$ is a material deformation
measure representing axial tension and shear,
$\underline{\bm{\Omega}}\in\mathbb{R}^{3}$ is a material deformation measure
representing torsion and bending, and
$\underline{\bm{C}}_{F}\in\mathbb{R}^{3\times 3}$ and
$\underline{\bm{C}}_{M}\in\mathbb{R}^{3\times 3}$ are cross-section
constitutive matrices. The material force stress resultants
$\underline{\bm{F}}=\partial\tilde{\Pi}_{\text{int},B}/\partial\underline{\bm{\Gamma}}$
and moment stress resultants
$\underline{\bm{M}}=\partial\tilde{\Pi}_{\text{int},B}/\partial\underline{\bm{\Omega}}$
can be derived from the hyperelastic stored-energy function. Finally, the beam
contributions to the weak form are given by

$\delta W^{B}=\delta\Pi_{\text{int},B}+\delta W^{B}_{\text{ext}},$ (10)

where $\delta W^{B}_{\text{ext}}$ is the virtual work of external forces and
moments acting on the beam. For a more comprehensive presentation of the weak
form of the geometrically exact Simo–Reissner beam theory, the interested
reader is referred to Meier2019 .

###  2.4 Beam-to-solid surface coupling (BTSS-FULL)

The BTSS-FULL method proposed in this work couples all six cross-section
degrees of freedom of the beam to the solid surface. This is realized by
coupling the positions of the beam centerline as well as the orientation of
the beam cross-section to the solid surface. One advantage of a 1D-2D coupling
approach solely enforced at the beam centerline is the decoupling of the
positional and rotational coupling conditions, i.e., both of them can be
formulated independently. For embedded 1D geometrically exact beams in 3D
solid volumes such an approach has recently been presented in Steinbrecher2022
. The same general strategy is followed here for BTSS coupling problems, where
we define two sets of coupling constraints, the positional coupling
constraints (BTSS-POS) and the rotational coupling constraints (BTSS-ROT).
With this split, the total BTSS-FULL coupling contribution to the weak form
reads,

$\delta\Pi_{\lambda}=\delta\Pi_{\lambda_{r}}+\delta\Pi_{\lambda_{\theta}},$
(11)

where $\delta\Pi_{\lambda_{r}}$ and $\delta\Pi_{\lambda_{\theta}}$ are the
virtual work contributions from the positional and rotational coupling
conditions, respectively.

####  2.4.1 Closest point projection

Figure 3: Illustration of possible BTSS coupling problems – 3 a curved beam on
a matching curved solid surface, 3 a curved beam centerline offset by the
cross-section radius in surface normal direction and 3 a general non-matching
case.

In the considered BTSS-FULL coupling problem, cf. Figure 2, no requirements on
the initial beam position relative to the solid surface exist. This is
illustrated in Figure 3. Obviously, the coupling scheme has to be applicable
to cases where the beam centerline curve lies on the solid surface, cf. Figure
3, and cases where the beam centerline is offset by the cross-section radius
in surface normal direction, cf. Figure 3. However, also general cases, where
no strict requirements on the reference placement of the beam centerline
relative to the solid surface are made, are considered in the presented
coupling schemes, cf. Figure 3. The only requirement considered in this work
is a unique closest point projection of each beam centerline point onto the
solid surface. For the envisioned application cases, it can be assumed that a
unique solution of the closest point projection exists in the vicinity of each
beam centerline point $\underline{\bm{r}}_{0}$, cf. Konyukhov2008 . In order
to formulate the closest point projection, the solid surface is parameterized
with the two surface parameter coordinates $\xi^{S}\in\mathbb{R}$ and
$\eta^{S}\in\mathbb{R}$. In the reference configuration each point
$\underline{\bm{r}}_{0}(s)$ on the beam centerline is assigned to a
corresponding closest point $\underline{\bm{X}}_{S}(\xi^{S}_{c},\eta^{S}_{c})$
on the solid surface, where $\xi^{S}_{c}=\xi^{S}_{c}(s)$ and
$\eta^{S}_{c}=\eta^{S}_{c}(s)$ are the surface parameter coordinates of the
closest point. The closest point can be found by formulating a unilateral
minimal distance problem in the reference configuration:

$d_{c,0}(s):=\min_{\xi^{S},\eta^{S}}d(s,\xi^{S},\eta^{S})=d(s,\xi^{S}_{c},\eta^{S}_{c})$
(12)

with

$d(s,\xi^{S},\eta^{S})=\left\|\underline{\bm{r}}_{0}(s)-\underline{\bm{X}}_{S}(\xi^{S},\eta^{S})\right\|.$
(13)

The two orthogonality conditions obtained from the minimal distance problem
(12) read

$\begin{split}\underline{\bm{X}}_{S,\xi^{S}}(\xi^{S},\eta^{S})^{\mathrm{T}}(\underline{\bm{r}}_{0}(s)-\underline{\bm{X}}_{S}(\xi^{S},\eta^{S}))=0,\\\
\underline{\bm{X}}_{S,\eta^{S}}(\xi^{S},\eta^{S})^{\mathrm{T}}(\underline{\bm{r}}_{0}(s)-\underline{\bm{X}}_{S}(\xi^{S},\eta^{S}))=0.\end{split}$
(14)

For a given beam coordinate $s$, these conditions can be solved for the
unknown surface coordinates $\xi^{S}$ and $\eta^{S}$. The non-trivial solution
of (14) requires the surface directors
$\underline{\bm{X}}_{S,\xi^{S}}=\partial\underline{\bm{X}}_{S}/\partial\xi^{S}$
and
$\underline{\bm{X}}_{S,\eta^{S}}=\partial\underline{\bm{X}}_{S}/\partial\eta^{S}$
to be orthogonal to the relative vector between the surface point and the beam
centerline point, i.e., this relative vector is parallel to the outward
pointing surface normal vector $\underline{\bm{N}}\in\mathbb{R}^{3}$,

$\underline{\bm{r}}_{0}(s)-\underline{\bm{X}}_{S}(\xi^{S}_{c},\eta^{S}_{c})=d_{c,0}(s)\underline{\bm{N}}(\xi^{S}_{c},\eta^{S}_{c})$
(15)

with

$\underline{\bm{N}}(\xi^{S},\eta^{S})=\frac{\underline{\bm{X}}_{S,\xi^{S}}(\xi^{S},\eta^{S})\times\underline{\bm{X}}_{S,\eta^{S}}(\xi^{S},\eta^{S})}{\left\|\underline{\bm{X}}_{S,\xi^{S}}(\xi^{S},\eta^{S})\times\underline{\bm{X}}_{S,\eta^{S}}(\xi^{S},\eta^{S})\right\|}.$
(16)

###  2.5 Positional beam-to-solid surface coupling (BTSS-POS)

In this section, three different variants of the BTSS-POS coupling constraints
are presented. They will be compared with each other in more detail in
LABEL:sec:examples. The first presented variant is consistent with the
kinematic relations between beam centerline and solid surface. The resulting
coupling terms contain the surface normal vector, i.e., the coupling terms
become non-linear. Furthermore, the second derivative of the surface normal
vector is required for a consistent linearization of the problem in tangent-
based nonlinear solution schemes (such as the Newton–Raphson algorithm). To
avoid this computationally expensive linearization, two additional variants to
formulate the positional coupling constraints, commonly used in classical
surface-to-surface mesh tying problems Puso2004 , will be investigated. Both
of them do not require an evaluation of the current surface normal vector or
its derivatives, and the resulting coupling operators only depend on the
reference configuration, i.e., they are constant. The different coupling
variants are visualized in Figure 4.

Figure 4: Illustration of the different positional beam-to-solid surface
(BTSS-POS) coupling variants. 4 Consistent positional coupling (BTSS-POS-CONS)
via the surface normal vector, 4 forced reference configuration coupling
(BTSS-POS-REF) by forcing beam centerline points to lie on the solid surface
and 4 displacement coupling (BTSS-POS-DISP), where the displacement of beam
centerline and solid surface are coupled. The surface-to-surface equivalents
of the BTSS-POS-REF and BTSS-POS-DISP variants are commonly used in classical
surface-to-surface mesh tying problems Puso2004 .

####  2.5.1 Consistent positional coupling (BTSS-POS-CONS)

The BTSS-POS coupling constraints are exclusively formulated along the beam
centerline and couple the beam and solid material points associated by (15) to
each other. For the considered consistent variant, the surface normal distance
$d_{c}$ at each beam centerline point shall be constant over the simulation
(pseudo-)time, i.e., $d_{c}\equiv d_{c,0}$. Therefore, the coupling equations
in the current configuration can be formulated as

$\underline{\bm{r}}(s)-\underline{\bm{x}}_{S}(\xi^{S}_{c},\eta^{S}_{c})-d_{c,0}(s)\,\underline{\bm{n}}(\xi^{S}_{c},\eta^{S}_{c})=\underline{\bm{0}}\quad\text{on}\quad\Gamma_{c}.$
(17)

The current normal vector is defined in analogy to the reference normal vector
(4), i.e.,

$\underline{\bm{n}}(\xi^{S}_{c},\eta^{S}_{c})=\frac{\underline{\bm{x}}_{S,\xi^{S}}(\xi^{S}_{c},\eta^{S}_{c})\times\underline{\bm{x}}_{S,\eta^{S}}(\xi^{S}_{c},\eta^{S}_{c})}{\left\|\underline{\bm{x}}_{S,\xi^{S}}(\xi^{S}_{c},\eta^{S}_{c})\times\underline{\bm{x}}_{S,\eta^{S}}(\xi^{S}_{c},\eta^{S}_{c})\right\|},$
(18)

with the current surface directors
$\underline{\bm{x}}_{S,\xi^{S}}=\partial\underline{\bm{x}}_{S}/\partial\xi^{S}$
and
$\underline{\bm{x}}_{S,\eta^{S}}=\partial\underline{\bm{x}}_{S}/\partial\eta^{S}$.
The constraints (17) are enforced along the one-dimensional coupling domain
$\Gamma_{c}\subseteq\Omega_{B,0}$ between the beam centerline and the solid
surface, i.e., the part of the beam that is coupled to the solid surface. In
the following considerations, the explicit dependency on the beam and solid
parameter coordinates will mostly be omitted for improved readability.

In the remainder of this work, the positional coupling constraints (17) will
be referred to as the _consistent_ (BTSS-POS-CONS) surface coupling variant.
The name refers to the fact that the coupling definition is consistent with
the kinematic relations between solid surface and beam centerline, cf. Figure
4. Furthermore, it will be shown that this variant leads to vanishing
constraint forces in the (undeformed) reference configuration and exact
conservation of linear and angular momentum in the discretized coupled system,
cf. LABEL:sec:discretization:conservation_laws.

The Lagrange multiplier method is used to weakly enforce the coupling
constraints (17). Therefore, a Lagrange multiplier vector field
$\underline{\bm{\lambda}}_{r}(s)\in\mathbb{R}^{3}$, defined along the beam
centerline, is introduced. The total Lagrange multiplier potential reads:

$\Pi_{\lambda_{r}}=\int_{\Gamma_{c}}{\underline{\bm{\lambda}}_{r}^{\mathrm{T}}(\underline{\bm{r}}-\underline{\bm{x}}_{S}-d_{c,0}\underline{\bm{n}})\,\mathrm{d}s}\,.$
(19)

Variation of the Lagrange multiplier potential leads to the constraint
contribution to the weak form,

$\begin{split}\delta\Pi_{\lambda_{r}}=&\underbrace{\int_{\Gamma_{c}}{\delta\underline{\bm{\lambda}}_{r}^{\mathrm{T}}(\underline{\bm{r}}-\underline{\bm{x}}_{S}-d_{c,0}\underline{\bm{n}})\,\mathrm{d}s}\,}_{\delta
W_{\lambda_{r}}}\\\
&+\underbrace{\int_{\Gamma_{c}}{\underline{\bm{\lambda}}_{r}^{\mathrm{T}}(\delta\underline{\bm{r}}-\delta\underline{\bm{x}}_{S}-d_{c,0}\delta\underline{\bm{n}})\,\mathrm{d}s}\,}_{-\delta
W_{C_{r}}}.\end{split}$ (20)

Therein, $\delta W_{\lambda_{r}}$ and $\delta W_{C_{r}}$ are the variational
form of the coupling constraints and the virtual work of the Lagrange
multiplier field $\underline{\bm{\lambda}}_{r}$, respectively. It is well-
known from geometrically exact beam theory that the variation of the
centerline position $\delta\underline{\bm{r}}$ is work-conjugated with the
point-wise forces acting on the beam centerline, i.e.,
$\underline{\bm{\lambda}}_{r}\mathrm{d}s$. Therefore, the Lagrange multiplier
field $\underline{\bm{\lambda}}_{r}$ can be directly interpreted as the
coupling line load acting on the beam centerline. On the solid side, the
variation of the solid displacement $\delta\underline{\bm{x}}_{S}$ is work
conjugated with the point-wise force
$-\underline{\bm{\lambda}}_{r}\mathrm{d}s$ acting on the solid, i.e., the
(negative) Lagrange multiplier field also acts as a line load on the solid.
Additionally, the term
$d_{c,0}\underline{\bm{\lambda}}_{r}^{\mathrm{T}}\delta\underline{\bm{n}}\mathrm{d}s$
arises, which represents a point-wise moment contribution of the coupling line
load on the solid. If the beam centerline lies exactly on the beam surface,
i.e., $d_{c,0}\equiv 0$, the BTSS-POS-CONS method (in the space continuous
form) is equivalent to the BTSV-POS method, cf. Steinbrecher2020 . The
drawback of the BTSS-POS-CONS variant is that for general scenarios the weak
form contains the surface normal vector variation, thus requiring the second
derivatives of the surface normal vector for a consistent linearization of
$\delta\underline{\bm{n}}$ as required for tangent-based nonlinear solution
schemes. Furthermore, the positional coupling operators become non-linear due
to this contribution, i.e., they depend on the current configuration.

####  2.5.2 Forced reference configuration coupling (BTSS-POS-REF)

The considered 1D-2D line-to-surface coupling constraints are very similar to
the ones in classical 2D-2D surface-to-surface coupling problems, cf. Puso2004
; Park2002 ; Dohrmann2000 . In such problems, the space continuous interfaces
are usually matching, i.e., the normal distance vanishes. Even if the surfaces
do not match exactly, e.g., due to incompatible CAD files, the influence of
the surface normal vector can usually be neglected since it is in the range of
the discretization error. Therefore, the coupling constraints (17) can be
simplified to

$\underline{\bm{r}}-\underline{\bm{x}}_{S}=\underline{\bm{0}}\quad\text{on}\quad\Gamma_{c}.$
(21)

This type of positional coupling constraint will be referred to as the _forced
reference configuration_ surface coupling (BTSS-POS-REF). The Lagrange
multiplier coupling contributions to the global weak form read

In this case, the surface normal vector is not contained in the resulting
coupling equations, thus simplifying the numerical evaluation of the coupling
terms. However, the coupling constraints (21) in the reference configuration
are only fulfilled if the beam centerline lies exactly on the solid surface,
i.e., $d_{c,0}\equiv 0$. If the beam centerline is not a subset of the solid
surface, the coupling constraints (21) will lead to non-vanishing virtual work
contributions in the reference configuration, i.e., initial stresses and
deformations in the unloaded coupled system. In other words, the BTSS-POS-REF
coupling conditions force the beam centerline to exactly lie on the solid
surface, which is illustrated in Figure 4.

####  2.5.3 Displacement coupling (BTSS-POS-DISP)

Another alternative coupling approach in surface-to-surface mesh tying is to
directly couple the displacements instead of the positions in (21). This
variant will be referred to as the _displacement_ surface coupling (BTSS-POS-
DISP). The BTSS-POS-DISP coupling constraints read,

$\underline{\bm{u}}_{B}-\underline{\bm{u}}_{S}=\underline{\bm{0}}\quad\text{on}\quad\Gamma_{c},$
(22)

with the beam centerline displacement
$\underline{\bm{u}}_{B}=\underline{\bm{r}}-\underline{\bm{r}}_{0}$. The
Lagrange multiplier coupling contributions to the global weak form are
$0\mathrm{e}{0}$$1.5\mathrm{e}{0}$$\left\|\underline{\bm{u}}\right\|$ Figure
21: Supported plate – deformed configurations for various modeling techniques.
LABEL:fig:examples:btssc_rotation_plate:results_deformed:full_large Full 3D
model, LABEL:fig:examples:btssc_rotation_plate:results_deformed:bts_full BTSS-
FULL-CONS (with rotational coupling) and
LABEL:fig:examples:btssc_rotation_plate:results_deformed:bts_trans BTSS-POS-
CONS (without rotational coupling). The contour plots visualize the
displacement magnitude.

$-0.5$$0$$0.5$$0.6$$0.8$$1$$1.2$$1.4$$\underline{\bm{e}}_{2}$$\underline{\bm{e}}_{3}$full
3DBTSS-FULL-CONSBTSS-POS-CONS Figure 22: Supported plate – deformed
configurations of the material line indicated in
LABEL:fig:examples:btssc_rotation_plate:problem for various modeling
techniques.

###  5.4 Towards biomedical applications

The last example is designed to give an outlook towards real-life applications
and the suitability of the proposed BTSS-FULL approach for more complex
coupling scenarios. Specifically, we want to analyze the applicability of our
approach in the context of vascular angioplasty. To this end, we set up a
variant of the well-known fluid-structure interaction (FSI) benchmark problem
of a pressure wave traveling through an elastic tube, that was originally
proposed in Gerbeau2003  to validate the suitability of FSI algorithms for
blood flow simulations. In addition to the original problem, we will use our
BTSS-FULL coupling approach to capture the effect of a diamond-shaped stent
structure on the behavior of the overall system. In particular, we expect to
capture the large change in compliance between the stented and unstented
regions of the pipe, thus leading to stress peaks in these transitional
regions as well as an altered fluid flow. Such effects have been linked to the
occurrence of in-stent restenosis and are of high interest when analyzing the
suitability of endovascular devices and their effect on the patient Kim2013 ;
Colombo2021 ; Kohler1992 .

Figure 23: Towards biomedical applications – problem setup. 23 Geometric configuration of the stent and artery, and 23 unwrapped stent geometry. Table 3: Table containing the parameters for the stented elastic pipe problem. Geometry | $r_{i}$ | $0.0125\,\mathrm{m}$
---|---|---
| $r_{a}$ | $0.01375\,\mathrm{m}$
| $r_{s}$ | $0.01246\,\mathrm{m}$
| $l_{1}$ | $0.15\,\mathrm{m}$
| $l_{2}$ | $0.06\,\mathrm{m}$
Beam | $E_{B}$ | $9\cdot 10^{8}\,\mathrm{N/m^{2}}$
| $\rho_{B}$ | $7800\,\mathrm{kg/m^{3}}$
| $\nu_{B}$ | $0.3$
| $R$ | $0.0004\,\mathrm{m}$
Solid | $E_{S}$ | $3\cdot 10^{5}\,\mathrm{N/m^{2}}$
| $\rho_{S}$ | $1200\,\mathrm{kg/m^{3}}$
| $\nu_{S}$ | $0.3$
Fluid | $p_{\text{in}}$ | $500\,\mathrm{N/m^{2}}$
| $\rho_{F}$ | $1000\,\mathrm{kg/m^{3}}$
| $\eta_{F}$ | $0.003\,\mathrm{kg/(ms)}$

As in the original benchmark problem, a constant pulse $p_{\text{in}}$ is
applied for $3\cdot 10^{-3}\,\mathrm{s}$ at the fluid inlet. Besides the
pulse, zero traction conditions are applied to the fluid inflow as well as
outflow boundary on the left and right end of the pipe, respectively, while
both ends of the pipe are assumed to be clamped. In addition to the BTSS-FULL
problem introduced in Section 2, this example contains a fluid, modeled as
Newtonian with a constant dynamic viscosity $\eta_{F}$ and a density
$\rho_{F}$, using the incompressible Navier-Stokes equations. Figure 23
illustrates the problem setup. The fluid is coupled to the solid via classical
surface-coupled FSI Kloeppel2011  in a partitioned manner aided by a matrix-
free Newton Krylov method Kuettler2008  to accelerate convergence. Classical
no-slip conditions are enforced on the FSI boundary. The beam centerline
geometry depicted in Figure 23 is wrapped around a cylinder with a radius of
$r_{s}=r_{i}-R$ to create the used diamond-shaped stent geometry such that the
stent perfectly fits into the pipe structure up to an offset the size of the
beam radius. Since FSI problems are necessarily transient, the BTSS-FULL
problem is enhanced by a Generalized-$\alpha$ Lie group time integration
method for all structural degrees of freedom Bruels2010 ; Bruels2012 . Here,
the parameters are chosen to obtain a fully implicit scheme, and a time step
size $\Delta t=10^{-4}\,\mathrm{s}$ is used. To the fluid field, a classical
second-order accurate Generalized-$\alpha$ time integration scheme, with the
same time step size as for the structure field, is applied Jansen2000 . The
mortar-type BTSS-FULL-CONS method is used with linear shape functions for the
Lagrange multiplier fields and the penalty parameters
$\epsilon_{r}=10^{9}\,\mathrm{N/m^{2}}$ and
$\epsilon_{\theta}=10^{-1}\,\mathrm{Nm/m}$. To discretize the problem, $264$
Reissner beam elements, $2{,}880$ solid shell elements and $22{,}800$
PSPG/SUPG stabilized Q1-Q1 fluid elements with an additional div-grad
stabilization term schott2015  are employed. All dimensions and material
parameters of the problem setup are summarized in Table 3.

Figures 24(a) to 24(d) depict the structural displacement scaled with a factor
of 15 and the fluid pressure after $0.01\,\mathrm{s}$, $0.016\,\mathrm{s}$,
$0.024\,\mathrm{s}$ and $0.030\,\mathrm{s}$. It is evident, that the wall
displacement caused by the pressure wave in the stiffer stented region in
Figure 24 is smaller than in the unstented region. This, in turn, affects the
fluid since a constant flow throughput requires increased velocities within
the stented region compared to the more compliant unstented region. Figure 25
illustrates the fluid velocity $v_{2}$ in channel direction along the pipe’s
centerline. The fluid velocity plot demonstrates the previously mentioned
phenomenon as the maximum fluid velocity increases slightly and the wave
broadens while traveling through the stented region. This effect on the fluid
flow is still visible even after the pressure wave leaves the stented region.

(a)(b)(c)(d)

$-1.9\mathrm{e}{2}$$4.3\mathrm{e}{2}$$p$$0\mathrm{e}{0}$$1.8\mathrm{e}{-4}$$\left\|\underline{\bm{u}}\right\|$

Figure 24: Towards biomedical applications – deformed configuration of the
stented elastic pipe problem at various simulation times. The snapshots are
taken at (a) $t=0.01\,\mathrm{s}$, (b) $t=0.016\,\mathrm{s}$, (c)
$t=0.024\,\mathrm{s}$ and (d) $t=0.030\,\mathrm{s}$ respectively. The norm of
the displacements is shown in the solid and the pressure is shown in the
fluid. The displacements are scaled with a factor of 15. $0$$5\cdot
10^{-2}$$0.1$$0.15$$-5$$0$$5$$\cdot
10^{-2}$$\underline{\bm{e}}_{2}$$v_{2}$0.010 s0.016 s0.024 s0.030 s Figure 25:
Plots of the fluid velocity $v_{2}$ in channel direction along the pipe’s
centerline.

While the change of compliance in the artery, and thus also its effect on the
fluid flow, could also be modeled by a simpler homogenized approach, the
proposed approach allows to quantify the forces interchanged on the coupling
interface. Figures 26 and 26 illustrate the coupling interactions, i.e., the
line loads excerted on the beam system by the surface. In general, it can be
observed that the interaction is highest at the ends of the stent, i.e., at
the transition between a compliant and a very stiff region. This is
particularly notable in Figure 26, where the pressure wave is right at the
transition between the unstented and stented region. Furthermore, dividing the
1D coupling loads by the beam diameter results in an approximation of the
interaction stresses between the beam and the artery. The maximum absolute
values of normal and shear stresses can be estimated for this example as
$1.8179\cdot 10^{3}\,\mathrm{N/m^{2}}$ and $1.28899\cdot
10^{3}\,\mathrm{N/m^{2}}$ (not visualized in the figures), respectively.

$2\mathrm{e}{-2}$$1.8\mathrm{e}{0}$$\left\|\underline{\bm{\lambda}}_{r}\right\|$

Figure 26: Negative coupling line loads for beam-to-solid surface coupling at
various simulation times – the snapshots are taken at (a)
$t=0.01\,\mathrm{s}$, (b) $t=0.016\,\mathrm{s}$, (c) $t=0.024\,\mathrm{s}$ and
(d) $t=0.030\,\mathrm{s}$ respectively. Only the positional Lagrange
multiplier field is shown and five values are visualized along each beam
element. The displacements are scaled with a factor of 15.

The demonstrated example certainly represents a simplified model. In
particular, the use of beam-to-solid surface coupling, as presented here,
instead of frictional beam-to-solid contact, prevents the observation of some
real-life phenomena such as stent migration. Nevertheless, because of growth
and remodeling of the artery and successive protrusion of the stent struts,
coupling, i.e., mesh tying, is a valid assumption in many patient-specific
cases. A further interesting novel computational method to incorporate was
recently reported in Hagmeyer2022 . It enables capturing the effect of the
stent struts on the fluid flow, which is linked to altered wall shear stresses
that may lead to in-stent restenosis Johari2020 ; Pant2010 . In any case, the
presented simulation results serve as a proof of concept to show that the
proposed BTSS-FULL coupling approach can generally be used for geometrically
complex beam systems such as stent geometries. The ability to capture
important phenomena, such as changes in compliance and its effect on the blood
flow as well as the distribution of the interaction forces, which may provide
insight into the long-term success of vascular angioplasty, has been
demonstrated.

##  6 Conclusion

In this work, we have proposed a 1D-2D mixed-dimensional coupling method to
consistently couple 1D Cosserat beams to 2D surfaces of 3D Boltzmann continua
(solid). Therein, the six coupling constraints act along the beam centerline,
i.e., three positional constraints and three rotational constraints. Three
different variants of the positional coupling constraints have been
investigated. One of them, the _consistent_ variant, requires the expensive
evaluation of the current surface normal field. The other two variants are
commonly used in surface-to-surface mesh tying problems. Numerical examples
have shown that only the consistent positional coupling constraints, i.e.,
with inclusion of the surface normal vector, lead to physically correct
results and fulfill basic mechanical consistency properties, such as
conservation of angular momentum. Existing coupling methods for the rotational
degrees of freedom are extended by constructing a suitable surface triad field
on the 2D surface of the 3D Boltzmann continuum. The Lagrange multiplier
method is used to enforce the positional and rotational coupling constraints.
The coupling equations are discretized using a mortar-type approach and the
resulting discrete constraint equations are regularized via a weighted penalty
approach. Furthermore, the numerical examples illustrate the importance of
combining both positional _and_ rotational coupling via a practically
motivated example. Finally, a multi-physics simulation, inspired by models of
stented arteries, has demonstrated the method’s suitability for complex beam
geometries and its ability to capture global effects on the solid as well as
the fluid field.

Future work will focus on the extension of the presented beam-to-surface
coupling approach to beam-to-surface contact and finite sliding problems,
i.e., replacing the coupling constants with unilateral and frictional contact
constraints.

###### Acknowledgements.

Sketches in this work have been created using the Adobe Illustrator plug-in
LaTeX2AI (https://github.com/isteinbrecher/LaTeX2AI).

## References

  * (1) Argyris, J.: An excursion into large rotations.  Computer Methods in Applied Mechanics and Engineering 32(1), 85–155 (1982)
  * (2) BACI: A Comprehensive Multi-Physics Simulation Framework.  https://baci.pages.gitlab.lrz.de/website (2021)
  * (3) Ben Belgacem, F.: The mortar finite element method with Lagrange multipliers.  Numerische Mathematik 84(2), 173–197 (1999)
  * (4) Betsch, P., Menzel, A., Stein, E.: On the parametrization of finite rotations in computational mechanics: A classification of concepts with application to smooth shells.  Computer Methods in Applied Mechanics and Engineering 155(3), 273–305 (1998)
  * (5) Betsch, P., Steinmann, P.: Frame-indifferent beam finite elements based upon the geometrically exact beam theory.  International Journal for Numerical Methods in Engineering 54(12), 1775–1788 (2002)
  * (6) Bischoff, M., Ramm, E.: Shear deformable shell elements for large strains and rotations.  International Journal for Numerical Methods in Engineering 40(23), 4427–4449 (1997)
  * (7) Brüls, O., Cardona, A.: On the use of Lie group time integrators in multibody dynamics.  Journal of Computational and Nonlinear Dynamics 5(3) (2010)
  * (8) Brüls, O., Cardona, A., Arnold, M.: Lie group generalized-$\alpha$ time integration of constrained flexible multibody systems.  Mechanism and Machine Theory 48, 121–137 (2012)
  * (9) Cardona, A., Geradin, M.: A beam finite element non-linear theory with finite rotations.  International Journal for Numerical Methods in Engineering 26(11), 2403–2438 (1988)
  * (10) Chang, T.Y., Taniguchi, H., Chen, W.F.: Nonlinear finite element analysis of reinforced concrete panels.  Journal of Structural Engineering 113(1), 122–140 (1987)
  * (11) Colombo, M., He, Y., Corti, A., Gallo, D., Ninno, F., Casarin, S., Rozowsky, J., Migliavacca, F., Berceli, S., Chiastra, C.: In-stent restenosis progression in human superficial femoral arteries: Dynamics of lumen remodeling and impact of local hemodynamics.  Annals of Biomedical Engineering 49, 2349–2364 (2021)
  * (12) Crisfield, M.A., Jelenić, G.: Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation.  Proceedings of the Royal Society of London A 455(1983), 1125–1147 (1999)
  * (13) Dohrmann, C.R., Key, S.W., Heinstein, M.W.: Methods for connecting dissimilar three-dimensional finite element meshes.  International Journal for Numerical Methods in Engineering 47(5), 1057–1080 (2000)
  * (14) Durville, D.: Finite element simulation of textile materials at mesoscopic scale.  In: Finite element modelling of textiles and textile composites. Saint-Petersbourg, Russian Federation (2007)
  * (15) Elwi, A.E., Hrudey, T.M.: Finite element model for curved embedded reinforcement.  Journal of Engineering Mechanics 115(4), 740–754 (1989)
  * (16) Farah, P., Popp, A., Wall, W.A.: Segment-based vs. element-based integration for mortar methods in computational contact mechanics.  Computational Mechanics 55(1), 209–228 (2015)
  * (17) Gerbeau, J.F., Vidrascu, M.: A quasi-Newton algorithm based on a reduced model for fluid-structure interaction problems in blood flows.  ESAIM: Mathematical Modelling and Numerical Analysis 37(4), 631–647 (2003)
  * (18) Gomes, H.M., Awruch, A.M.: Some aspects on three-dimensional numerical modelling of reinforced concrete structures using the finite element method.  Advances in Engineering Software 32(4), 257–277 (2001)
  * (19) Hagmeyer, N., Mayr, M., Steinbrecher, I., Popp, A.: One-way coupled fluid-beam interaction: Capturing the effect of embedded slender bodies on global fluid flow and vice versa.  Advanced Modeling and Simulation in Engineering Sciences 9(1), 9 (2022)
  * (20) Ibrahimbegović, A., Frey, F., Kožar, I.: Computational aspects of vector-like parametrization of three-dimensional finite rotations.  International Journal for Numerical Methods in Engineering 38(21), 3653–3673 (1995)
  * (21) Jansen, K., Whiting, C., Hulbert, G.: Generalized-$\alpha$ method for integrating the filtered Navier-Stokes equations with a stabilized finite element method.  Computer Methods in Applied Mechanics and Engineering 190, 305–319 (2000)
  * (22) Johari, N., Hamady, M., Xu, X.: A computational study of the effect of stent design on local hemodynamic factors at the carotid artery bifurcation.  Artery Research 26, 161 – 169 (2020)
  * (23) Kerfriden, P., Claus, S., Mihai, I.: A mixed-dimensional CutFEM methodology for the simulation of fibre-reinforced composites.  Advanced Modeling and Simulation in Engineering Sciences 7(1), 18 (2020)
  * (24) Khristenko, U., Schuß, S., Krüger, M., Schmidt, F., Wohlmuth, B., Hesch, C.: Multidimensional coupling: A variationally consistent approach to fiber-reinforced materials.  Computer Methods in Applied Mechanics and Engineering 382, 113869 (2021)
  * (25) Kim, Y.G., Oh, I.Y., Kwon, Y.W., Han, J.K., Yang, H.M., Park, K.W., Lee, H.Y., Kang, H.J., Koo, B.K., Kim, H.S.: Mechanism of edge restenosis after drug-eluting stent implantation.  Circulation Journal 77, 2928–2935 (2013)
  * (26) Klöppel, T., Popp, A., Küttler, U., Wall, W.A.: Fluid–structure interaction for non-conforming interfaces based on a dual mortar formulation.  Computer Methods in Applied Mechanics and Engineering 200(45), 3111–3126 (2011)
  * (27) Kohler, T.R., Jawień, A.: Flow affects development of intimal hyperplasia after arterial injury in rats.  Arteriosclerosis and thrombosis: A journal of vascular biology 12, 963–71 (1992)
  * (28) Konyukhov, A., Schweizerhof, K.: On the solvability of closest point projection procedures in contact analysis: Analysis and solution strategy for surfaces of arbitrary geometry.  Computer Methods in Applied Mechanics and Engineering 197(33), 3045–3056 (2008)
  * (29) Konyukhov, A., Schweizerhof, K.: On some aspects for contact with rigid surfaces: Surface-to-rigid surface and curves-to-rigid surface algorithms.  Computer Methods in Applied Mechanics and Engineering 283, 74–105 (2015)
  * (30) Korelc, J., Wriggers, P.: Automation of finite element methods.  Springer International Publishing (2016)
  * (31) Küttler, U., Wall, W.A.: Fixed-point fluid-structure interaction solvers with dynamic relaxation.  Computational Mechanics 43(1), 61–72 (2008)
  * (32) Meier, C., Grill, M.J., Wall, W.A.: Generalized section-section interaction potentials in the geometrically exact beam theory (2021).  Preprint, https://arxiv.org/abs/2105.10032
  * (33) Meier, C., Popp, A., Wall, W.A.: An objective 3D large deformation finite element formulation for geometrically exact curved Kirchhoff rods.  Computer Methods in Applied Mechanics and Engineering 278, 445–478 (2014)
  * (34) Meier, C., Popp, A., Wall, W.A.: A finite element approach for the line-to-line contact interaction of thin beams with arbitrary orientation.  Computer Methods in Applied Mechanics and Engineering 308, 377–413 (2016)
  * (35) Meier, C., Popp, A., Wall, W.A.: Geometrically exact finite element formulations for slender beams: Kirchhoff–Love theory versus Simo–Reissner theory.  Archives of Computational Methods in Engineering 26(1), 163–243 (2019)
  * (36) Pant, S., Bressloff, N.W., Forrester, A.I.J., Curzen, N.: The influence of strut-connectors in stented vessels: A comparison of pulsatile flow through five coronary stents.  Annals of Biomedical Engineering 38, 1893–1907 (2010)
  * (37) Park, K.C., Felippa, C.A., Rebel, G.: A simple algorithm for localized construction of non-matching structural interfaces.  Int. J. Numer. Meth. Engng. 53(9), 2117–2142 (2002)
  * (38) Phillips, D.V., Zienkiewicz, O.C.: Finite element non-linear analysis of concrete structures.  Proceedings of the Institution of Civil Engineers 61(1), 59–88 (1976)
  * (39) Podio-Guidugli, P.: Examples of concentrated contact interactions in simple bodies.  Journal of Elasticity 75(2), 167–186 (2005)
  * (40) Podio-Guidugli, P., Favata, A.: Elasticity for geotechnicians: A modern exposition of Kelvin, Boussinesq, Flamant, Cerruti, Melan, and Mindlin problems, vol. 204.  Springer, Cham (2014)
  * (41) Popp, A., Gee, M.W., Wall, W.A.: A finite deformation mortar contact formulation using a primal–dual active set strategy.  International Journal for Numerical Methods in Engineering 79(11), 1354–1391 (2009)
  * (42) Popp, A., Gitterle, M., Gee, M.W., Wall, W.A.: A dual mortar approach for 3D finite deformation contact with consistent linearization.  International Journal for Numerical Methods in Engineering 83(11), 1428–1465 (2010)
  * (43) Puso, M.A.: A 3D mortar method for solid mechanics.  International Journal for Numerical Methods in Engineering 59(3), 315–336 (2004)
  * (44) Puso, M.A., Laursen, T.A.: A mortar segment-to-segment contact method for large deformation solid mechanics.  Computer Methods in Applied Mechanics and Engineering 193(6), 601–629 (2004)
  * (45) Puso, M.A., Laursen, T.A.: A mortar segment-to-segment frictional contact method for large deformations.  Computer Methods in Applied Mechanics and Engineering 193(45), 4891–4913 (2004)
  * (46) Ranjbaran, A.: Mathematical formulation of embedded reinforcements in 3D brick elements.  Communications in Numerical Methods in Engineering 12(12), 897–903 (1996)
  * (47) Reissner, E.: On one-dimensional finite-strain beam theory: The plane problem.  Zeitschrift für angewandte Mathematik und Physik ZAMP 23(5), 795–804 (1972)
  * (48) Romero, I.: The interpolation of rotations and its application to finite element models of geometrically exact rods.  Computational Mechanics 34(2), 121–133 (2004)
  * (49) The Sacado Project Website.  https://trilinos.github.io/sacado.html (2021)
  * (50) Schott, B., Rasthofer, U., Gravemeier, V., Wall, W.A.: A face-oriented stabilized Nitsche-type extended variational multiscale method for incompressible two-phase flow.  International Journal for Numerical Methods in Engineering 104(7), 721–748 (2015)
  * (51) Simo, J.C., Vu-Quoc, L.: On the dynamics of flexible beams under large overall motions – The plane case: Part I.  Journal of Applied Mechanics 53(4), 849–854 (1986)
  * (52) Simo, J.C., Vu-Quoc, L.: On the dynamics of flexible beams under large overall motions – The plane case: Part II.  Journal of Applied Mechanics 53(4), 855–863 (1986)
  * (53) Sonneville, V., Cardona, A., Brüls, O.: Geometrically exact beam finite element formulated on the special Euclidean group SE(3).  Computer Methods in Applied Mechanics and Engineering 268, 451–474 (2014)
  * (54) Spurrier, R.A.: Comment on ”Singularity-free extraction of a quaternion from a direction-cosine matrix”.  Journal of Spacecraft and Rockets 15(4), 255–255 (1978)
  * (55) Steinbrecher, I., Mayr, M., Grill, M.J., Kremheller, J., Meier, C., Popp, A.: A mortar-type finite element approach for embedding 1D beams into 3D solid volumes.  Computational Mechanics 66(6), 1377–1398 (2020)
  * (56) Steinbrecher, I., Popp, A.: MeshPy – A general purpose 3D beam finite element input generator.  https://compsim.gitlab.io/codes/meshpy (2021)
  * (57) Steinbrecher, I., Popp, A., Meier, C.: Consistent coupling of positions and rotations for embedding 1D Cosserat beams into 3D solid volumes.  Computational Mechanics 69(3), 701–732 (2022)
  * (58) Taylor, R.L., Simo, J.C., Zienkiewicz, O.C., Chan, A.C.H.: The patch test – A condition for assessing FEM convergence.  International Journal for Numerical Methods in Engineering 22(1), 39–62 (1986)
  * (59) Thomson, W.: Note on the integration of the equations of equilibrium of an elastic solid.  The Cambridge and Dublin mathematical journal 3, 87–89 (1848)
  * (60) The Trilinos Project Website.  https://trilinos.github.io (2021)
  * (61) Vetyukov, Y.: Nonlinear mechanics of thin-walled structures: Asymptotics, direct approach and numerical analysis.  Foundations of Engineering Mechanics. Springer (2014)
  * (62) Vu-Quoc, L., Tan, X.G.: Optimal solid shells for non-linear analyses of multilayer composites. I. Statics.  Computer Methods in Applied Mechanics and Engineering 192(9), 975–1016 (2003)
  * (63) Wohlmuth, B.I.: A mortar finite element method using dual spaces for the lagrange multiplier.  SIAM Journal on Numerical Analysis 38(3), 989–1012 (2000)
  * (64) Yang, B., Laursen, T.A., Meng, X.: Two dimensional mortar contact methods for large deformation frictional sliding.  International Journal for Numerical Methods in Engineering 62(9), 1183–1225 (2005)


# Enhancing Affective Representations of Music-Induced EEG through Multimodal
Supervision and latent Domain Adaptation

###### Abstract

The study of Music Cognition and neural responses to music has been invaluable
in understanding human emotions. Brain signals, though, manifest a highly
complex structure that makes processing and retrieving meaningful features
challenging, particularly of abstract constructs like affect. Moreover, the
performance of learning models is undermined by the limited amount of
available neuronal data and their severe inter-subject variability. In this
paper we extract efficient, personalized affective representations from EEG
signals during music listening. To this end, we employ music signals as a
supervisory modality to EEG, aiming to project their semantic correspondence
onto a common representation space. We utilize a bi-modal framework by
combining an LSTM-based attention model to process EEG and a pre-trained model
for music tagging, along with a reverse domain discriminator to align the
distributions of the two modalities, further constraining the learning process
with emotion tags. The resulting framework can be utilized for emotion
recognition both directly, by performing supervised predictions from either
modality, and indirectly, by providing relevant music samples to EEG input
queries. The experimental findings show the potential of enhancing neuronal
data through stimulus information for recognition purposes and yield insights
into the distribution and temporal variance of music-induced affective
features.

22footnotetext: This research work was performed at NTUA.

Index Terms—  Music Cognition, Emotion Recognition, Electroencephalography,
Cross-Modal Learning

##  1 Introduction

Music is an abstract, yet densely emotional form of art. It is universally
enjoyed, due to its ability to induce powerful emotions irrespective of the
underlying mood [1] and has been characterized to greatly affect the function
of the human brain [2, 3]. Hence it is widely used to study emotion
recognition, both by analyzing the mood produced by several musical features
[4, 5] and by studying its effects on human neural and physiological responses
[6]. However, the task of extracting emotional information from brain activity
poses severe challenges, due to the inherently abstract nature of the induced
emotions, the variability in emotional and physiological responses between
different individuals and the lack of large-scale databases of emotionally
coordinated neural activity. In this paper, we propose a deep multimodal
approach [7], using musical stimuli.

The scope of this study lies at the intersection of Music Cognition, Emotion
Recognition from neuronal signals and Multimodal Learning. We choose to study
brain responses to music by employing a cross-modal system to identify the
correspondence between these modalities. We use the electroencephalogram (EEG)
to model brain responses for this task and we constrain the learning process
with emotion labels. Therefore, we aim to derive important insights regarding
the affective role that music can play on humans and the extent to which it
can help us build richer neuronal representations of affect. To conduct the
experiment, we exploit multimodal optimization and domain adaptation
strategies to project EEG and music features onto a common latent space, from
which we could assess their similarity. By conditioning the learning process
with emotion tags, the constructed space represents affect, enabling thus
emotion recognition both directly, by performing supervised predictions, and
indirectly, by ranking music tracks to EEG inputs, based on their distance. To
the best of our knowledge, this is the first study to propose such a
framework, thus it could be utilized as a baseline reference. We also perform
an extensive qualitative study across 32 subjects of the DEAP dataset [8] to
derive inter-subject affective patterns.

The remainder of this paper is organized as follows: Section 2 reviews the
related work in cross-modal learning and research on EEG processing and music
cognition. In Section 3 we introduce our problem and present the proposed
framework along with the optimization methods we utilized. Section 4 includes
information about the data, their pre-processing and implementation details,
while in Section 5 we provide the experimental results. In that Section and in
Section 6 we provide extensive quantitative and qualitative analysis of the
outcomes, while Section 7 concludes our study.

##  2 Related Work

Music Cognition: Studying the human brain’s responses to music stimuli has
always been a lively field of research in neuroscience and signal processing
[9] aiming to answer fundamental questions regarding our enjoyment of music.
The field has gained a lot of attention in recent years, with the upsurge in
available neuronal data. Many studies in the field rely on EEG recordings as
they provide better temporal resolution than other techniques, such as
functional magnetic resonance imaging (fMRI). In addition to the traditional,
well-controlled auditory experiments, modern approaches gather physiological
data from music listeners as they enjoy or imagine naturalistic music [10], in
order for instance to examine correlations in temporal structure [11] or the
perceived tempo [12]. One of the core findings on music cognition is the
correlation between characteristics of the neural oscillation patterns and
rhythmical patterns in music [13]. Additionally, Event-Related Potentials
(ERP) have been utilized to extract brain activity patterns that can relate to
note onsets or pitch [14, 15]. In parallel to the above, there has also been a
shift towards deep learning approaches for information retrieval from music
stimuli [16], in which we focus on the present study.

Fig. 1: The proposed bi-stream network. The output embedding layer of each
stream is fed to the common 64D dense layer (common space).

Emotion Recognition: Undoubtedly, the most powerful impact of music on humans
concerns the induced emotions. Emotion Recognition is a widely researched
field of contemporary Machine Learning and Behavioral Signal Processing [17]
and several studies have focused on the musical features [5, 18] that
determine affective attributes of music listening in a wide range of emotions.

Recently, several studies have examined physiological signals to analyze
humans’ felt emotions [19], with music emerging as an efficient method to
elicit them. Due to its temporal resolution, the electroencephalogram is the
most widely researched signal of this type and various statistical, spectral
or time-frequency features have been proposed for Emotion Recognition [20,
21]. Due to the noisy structure of EEG signals, many studies incorporate
entropy [22] and fractal [23] algorithms to extract emotion-related features.
Of course, variations of deep neural networks have been proposed and exceeded
the performance of traditional feature extraction methods [24, 25], however
the limited data availability and inter-subject variability present serious
barriers for this kind of modeling.

Cross-Modal Learning: The task of learning a shared embedding space from
different datasets or modalities is being studied through various approaches,
which are predominantly applied to image and text modalities. A widely used
baseline is Canonical Correlation Analysis (CCA). CCA is non-probabilistic and
enables the extraction of linear components to optimize the correlation of
pairs of vectors. One can find in the literature various non-linear CCA-based
frameworks and architectures utilized to learn inter-modal similarities, such
as Deep CCA [26]. Besides CCA, other methods that have been used include an
HGR-based maximal correlation metric [27] and adversarial training [28],
focusing mainly on the optimization function of the respective model, and on
adaptive hidden layers [29]. Another study [30] incorporated music to co-train
a shared space with images using a contrastive loss. Further, in [31] a state-
of-the-art framework exploits label supervision to better manipulate the
latent space, a key concept that we also follow in our study.

##  3 Methodology

In this study, we extract the semantic relationship between music tracks and
corresponding EEG recordings, so that an EEG could be mapped to an efficient
affective representation and retrieve emotionally consistent music samples.
Let us assume a collection of $n$ instances of EEG-music pairs, denoted as
$T=\\{(x_{i}^{a},x_{i}^{b})\\}_{i=1}^{n}$ where $x_{i}^{a}$ is the input EEG
sample of the ith instance and $x_{i}^{b}$ the input music stimulus
corresponding to that sample. Each instance has been assigned an affective
annotation $y_{i}\in\mathbb{R}^{2}$ for valence and arousal dimensions. For
each instance $i$ we aim to learn an EEG embedding
$u_{i}=f(x_{i}^{a},Y^{a})\in\mathbb{R}^{d}$ and a musical audio embedding
$v_{i}=g(x_{i}^{b},Y^{b})\in\mathbb{R}^{d}$, where $d$ is the dimensionality
of the common representation space and $Y^{a},Y^{b}$ the trainable parameters.

###  3.1 The Proposed Framework

We use a bi-stream Neural Network with one branch for each modality. The EEG
branch is a recurrent network, comprised of two LSTM modules and a softmax
attention layer. The model takes as input an EEG trial of shape (channels,
features) and attempts to capture its inter-channel correlations. Next, the
output features of all channels are flattened and passed through an attention
module to identify the most important components, that will lead the embedding
vector in the common space. We utilize a lightweight network in order to avoid
overfitting to the limited range of the available data, however any state-of-
the-art model in the task could be applied. For the music branch we use the
MusiCNN model [32] to extract high-level embeddings from the available audio
stimuli. MusiCNN is a robust network, pre-trained on large audio databases,
and produces high-quality music embeddings that compensate for the limited
size of our track set and further assist the learning process. The extracted
embeddings are then fed into a feed-forward neural network.

To construct the final bi-modal framework, we connect the last layer of each
of the previous networks to a dense layer (Fig. 1) constituting the common
representation space, from which we output emotion predictions. Inspired by
[25], we further apply a Gradient Reversal Layer (GRL) [33], aiming to reduce
the distribution shift between EEG and music modalities. In specific, both 64D
embeddings are fed to this layer, from where we output a prediction regarding
the modality type. From each batch, we randomly permute half of the EEG
samples and their respective music samples, forming a new equal-sized batch
that we shuffle and input to the GRL module, along with a binary label vector
to denote the modality. Subsequently, these embeddings are passed through
dense layers to predict the modality of each sample. By reversing the
gradients corresponding to these predictions during back-propagation, we help
the feature extractor produce modality-invariant features.

###  3.2 Objective Function

Our goal is to learn a common space where the samples from the same semantic
category should be similar, even though they come from different modalities.
To learn discriminative features we want to minimize the discrimination loss
in both the label space and the representation space, by reducing the cross-
modal discrepancy. With regard to the label space, we use a linear classifier
to predict the emotion labels of the samples projected in the common space.
Outputs of each modality are passed through a sigmoid activation and a binary
cross-entropy (BCE) loss $\ell$ is computed. For the cross-modal task we apply
a weighted linear combination of those losses:
$\mathcal{J}_{1}=\lambda_{11}\ell_{a}+\lambda_{12}\ell_{b}$. To reduce the
cross-modal discrepancy between EEG and music representations, we also compute
the BCE loss of the modality prediction after GRL:
$\mathcal{J}_{2}=\ell_{dd}$. By combining terms
$\mathcal{J}_{1},\mathcal{J}_{2}$ we obtain the proposed objective, in which
the hyper-parameters $\lambda_{i}$ control the contribution of each separate
component and are determined through trial and error:
$\mathcal{J}=\lambda_{1}\mathcal{J}_{1}+\lambda_{2}\mathcal{J}_{2}$.

##  4 Experimental Setup

###  4.1 The DEAP Dataset

DEAP [8] is a comprehensive dataset that includes EEG signals of music
listening, collected from 32 subjects. Each subject watches forty 1-min long
music videos while having their EEG recorded. After each video trial, the
subject was instructed to rate the emotion that was elicited upon the entire
trial in 5 dimensions: valence, arousal, dominance, liking and familiarity to
the track. In this paper we solely experiment with the 2D emotion space,
determined by valence and arousal, whose ratings range from 1 (weakest) to 9
(strongest). We use the EEG signals in their already preprocessed form:
recorded at a sampling rate of 512 Hz and denoised by bandpass filtering,
after downsampling to 128 Hz. Eye-related artefacts were removed whereas the
10-20 electrode placement system was followed.

Specifying Music Tracks: The 40 one-minute music stimuli of DEAP are not
included in the dataset, so we located the video clips of the corresponding
tracks and isolated the minute of interest for each one, according to the
metadata provided. The task of deriving the common latent space poses a
crucial challenge: the semantic gap between the “subjective” affective
responses of participants and the emotion tags of the songs. Ideally, we need
musical stimuli that are tagged in accordance with the participants’
annotations. DEAP stimuli have been selected for this purpose and have been
independently annotated by the experimenters at track level. Nearly all songs
received average ratings from the participants that were in accordance with
those annotations. We found that only $6/40$ songs had such an inconsistency
and discarded them. The resulting track set is used to extract MusiCNN
embeddings which we make
available111https://github.com/klean2050/EEG_CrossModal as well.

###  4.2 Input Feature Extraction

EEG and music signals are processed differently in order to produce an
embedding form suitable for multimodal training. DEAP signals are first cut to
3-second segments, while having their preparatory phase removed. For feature
extraction purposes we consider differential entropy features, reported to
achieve superior performance in the task [34]. Differential entropy (DE) $h$
is defined as:

$h(X)=-\int_{X}f(x)\log(f(x))dx,$ (1)

where $X$ is an EEG segment and $f(x)$ its distribution. Assuming further that
the utilized signals can be modeled as Gaussian distributions, i.e.
$f(x)=N(\mu,\sigma)$, then $h(X)$ can be determined by the logarithm energy
spectrum of $X$ as follows:

$\begin{gathered}h(X)=-\int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\frac{(x-\mu)^{2}}{2\sigma^{2}}\\\
\log\left(\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\frac{(x-\mu)^{2}}{2\sigma^{2}}\right)dx=\frac{1}{2}\log
2\pi e\sigma^{2}.\end{gathered}$ (2)

Thus, for each EEG segment we use the Short-Time Fourier Transform with an
1-sec non-overlapping Hanning window to compute the variance $\sigma^{2}$ for
each of the three windows in the frequency domain and subsequently we compute
$h$ in each channel of the four available EEG rhythms: $\theta$ (4-7Hz),
$\alpha$ (7-13Hz), $\beta$ (13-30Hz) and $\gamma$ (31-50Hz). The features for
all four bands are concatenated and the resulting feature vector is then used
as channel-wise input. On the other hand, music tracks are cut in segments
aligned with the EEG throughout their whole (3 sec) duration and fed directly
into the pre-trained model, from which we extract its “pool5” embeddings.

###  4.3 Evaluation Protocol

We evaluate our proposed method using accuracy to assess the supervised
predictions for each modality and the Precision@10 (P@10) and mean Average
Precision (mAP) metrics for the retrieval of music tracks given EEG queries.
Those two metrics have been widely used to assess retrieval tasks in the
literature [35, 31] as they evaluate the response’s distance-based ranking to
each query. In particular, P@10 considers the top 10 ranked tracks whereas mAP
evaluates the whole ranking. Results are also presented after trial
aggregation: For the accuracy, we simply denote a prediction as correct by
majority voting on the segment-wise predictions. For the retrieval metrics,
since no such voting can be made, we consider the median of the segment-wise
distance scores as the overall query score.

##  5 Experiments

The training procedure considers personalized models, each one trained on data
of a specific subject and the respective audio stimuli. To compensate for
possible annotation noise, we binarize the emotion labels by setting the
threshold to the median score 5, as in [8]. Following the same paradigm, we
consider separate experiments for valence and arousal dimensions. We apply
5-fold stratified cross-validation to train each network, where each fold
holds 20% of the total trials (7 tracks). Additionally, we apply class weights
to alleviate any subject-specific data imbalance. All networks are optimized
using Adam at a $10^{-4}$ learning rate and patience of 15 epochs of non-
decreasing validation loss.

###  5.1 Predicting Emotion Tags

We first evaluate the models’ performance on Emotion Recognition for both EEG
and Music modalities (Table 1). EEG scores show high variance per subject, on
average reaching up to $70.4\%$ on valence and $68.9\%$ on arousal after trial
aggregation. The obtained scores are competitive for the specific dataset,
despite the simple utilized architecture, something we attribute to the impact
of music co-training and the adaptation of the common latent space. This
contribution is further quantified in Section 5.3. Additionally, aggregating
predictions on a per-track basis provides substantially enhanced results
compared to non-aggregated ones, with the EEG accuracy increasing by above 5%
in arousal recognition and about 8% in valence, implying that there is strong
correlation (e.g., in the form of clusters) between same-track samples,
especially in valence. On the other hand, despite the small number of tracks
in our music set, the recognition performance is substantially high for the
music branch, 78.8% average on valence and 91.9% on arousal, something that
indicates the robustness of our transfer learning module.

Dimension | Non-Aggregated | Aggregated
---|---|---
Valence | 62.9% – 71.5% |  70.4% – 78.7%
Arousal | 63.3% – 88.0% |  68.9% – 91.9%
Table 1: Emotion Accuracy Scores for (EEG – Music) modalities, reporting mean values over 32 subject-specific models. Dimension | Precision@10 | mAvg. Precision
---|---|---
Valence |  19.4% – 63.8% | 18.8% – 59.1%
Arousal | 18.4% – 65.0% |  19.9% – 67.8%
Table 2: (Track – Emotion) Retrieval Scores on EEG input queries, reporting
mean aggregated scores over 32 subjects.

Fig. 2: t-SNE visualisation of the common space for subjects –from left to
right– 8, 15 (Valence) and 18, 20 (Arousal). 0
$\rightarrow\text{Low}\,|\,1\rightarrow$ High

###  5.2 Retrieving Tracks from EEG Queries

Table 2 summarizes the retrieval scores from the personalized models, acquired
by querying the common representation space of each network with a test EEG
sample and then evaluating the ranking of music samples based on their
distance to the query. Retrieval metrics provide robust results in both cases,
indicating that the EEG samples are well-situated in the common space and the
majority of them are capable of retrieving tracks that are emotionally
consistent. Specifically, in the case of induced valence, a P@10 value of
63.8% is achieved. We note that this percentage is higher than the reported
mAP (59.1%), strongly implying that the learned valence space is fragmented
into local subspaces of high similarity. Arousal on the other hand seems to be
more consistently represented, as both mAP and P@10 median retrieval scores
indicate that the majority of tested tracks can derive consistent music
rankings, in contrast to valence where the emotional response similarity seems
concentrated to the top-ranked elements. As a result, the correct retrieval
percentage conditioned on arousal approaches 68% on average across subjects.
We also note some preliminary results in approaching retrieval of the exact
stimulus of an EEG sample. The derived scores, around 20%, are clearly above
random selection, however we believe that further experimentation is required
on the temporal resolution of the input samples, yielding an interesting
direction of future study.

###  5.3 Ablation Study

In our study we incorporated a complex objective function, combining 3 BCE
terms to minimize the discrimination loss in both the label space and the
common latent space. To further investigate the impact of our proposals on the
models’ performance, we trained separate sessions, first by considering sole
EEG samples without music supervision, and second by avoiding the domain
discrimination module. From the results in Table 3 we deduce that our full
objective $\mathcal{J}$ leads to higher overall performance, indicating that
all utilized terms contribute to richer EEG affective representations.
Specifically, we can see that the absence of multimodal training sharply
impacts the validity of the common space and reduces the classification
performance, 2.6% in valence and 0.9% in arousal. On the other side, the
absence of domain adaptation causes slighter modifications to the correlation
of samples and stimuli, as measured by precision metrics. Through this module
we manage though to better distribute samples in the common space, break
modality-specific clusters and reduce the overall sample distances (Section
6.1), which is reflected in the improved classification performance in both
experiments.

Metric | $\mathcal{J}$ |  $\ell_{a}$ only | $\neg\,\ell_{dd}$
---|---|---|---
Acc${}_{\text{EEG}}$ |  70.4% – 68.9% | 67.8% – 68.0% | 67.9% – 63.4%
P@10 |  63.8% – 65.0% | 57.3% – 53.1% |  63.4% – 66.7%
mAP | 59.1% – 67.8% | 51.9% – 55.8% |  59.8% – 68.1%
Table 3: Ablation on the Objective Function for (Valence – Arousal). Here we
solely consider mean aggregated scores over 32 subjects.

Fig. 3: Arousal mAP scores over the 58 time samples for the numbered tracks,
averaged across all subjects.

##  6 Qualitative Analysis

Studying the Common Space: We visually inspect the produced latent space using
t-SNE to reduce its 64D dimension to 2D. We select one of the 5 trained models
for 2 subjects in Valence and Arousal to display their results in Fig. 2
(similar trends are observed for most subjects). It is evident that latent
domain adaptation has alleviated the cross-modal discrepancy and the
modalities homogenize their embeddings to a certain degree. Cohesive sub-
clusters are though visible, especially in the case of valence. This provides
an explanation towards the discrepancy we observed between P@10 and mAP
metrics, since the top-ranked track retrievals originate from the
corresponding local subspace, but there is no coarse bisection between high-
and low- valence samples, in contrast to the case of arousal.

Temporal Variation of Recognition: Since each track is segmented into 58
overlapping samples of 3 sec, it is expected that the emotion is not elicited
at the same pace throughout its duration. Hence, the temporal variation of our
scores could indicate important moments in the track. In Fig. 3, we present
the temporal evolution of the mAP scores for selected music tracks, averaged
across all subjects. While the raw plots are noisy, each song individually
exhibits a pattern of variation, which we depict by applying a 7-sample moving
average filter. Scores typically reveal an oscillating pattern on the time
axis and emotions are highly induced at certain peaks of the graph. These
patterns reveal a characteristic picture of emotional induction in songs and
could be subject of further experimentation.

##  7 Conclusion

In this paper we presented a novel approach to analyze emotion induction from
EEG recordings of music listening. We proposed a cross-modal framework to
learn rich affective representations for EEG data through music supervision
and adaptation of a common latent space, from which one could retrieve
consistent music rankings from EEG queries. Our approach indicates that
distilling information from processed musical stimuli to the respective EEG
signals can improve performance and provide insights on personalized emotion
analysis. To the best of our knowledge, this is the first study to propose a
complete framework for the specific task and dataset, thus our results can be
viewed as a concrete baseline. This framework can be used to model the EEG-
Music relationship by using different condition mechanisms, e.g., musical
beat. Another interesting direction would be to explore improvements in exact
stimulus retrieval.

## References

  * [1] K. Jakubowski and A. Ghosh,  “Music-evoked autobiographical memories in everyday life,”  Psychology of music, vol. 49, no. 3, pp. 649–666, 2021.

  * [2] O. W. Sacks,  Musicophilia : Tales of Music and the Brain,  Alfred A. Knopf, New York, 1st edition, 2007.
  * [3] D. Bashwiner,  “Brain and Music. By Stefan Koelsch,”  Music Theory Spectrum, vol. 39, no. 2, pp. 269–274, 2017.
  * [4] R. Panda, R. Malheiro, and R. P. Paiva,  “Novel Audio Features for Music Emotion Recognition,”  IEEE Trans. Affective Computing, 2020.
  * [5] Y. Song, S. Dixon, and M. Pearce,  “Evaluation of Musical Features for Emotion Classification,”  in Proc. ISMIR 2012, Porto, Portugal, 2012.
  * [6] T. Greer, B. Ma, M. Sachs, A. Habibi, and S. Narayanan,  “A Multimodal View into Music’s Effect on Human Neural, Physiological, and Emotional Experience,”  in Proc. ACM Int’l Multimedia Conf. 2019, Nice, France, 2019.
  * [7] J. Ngiam, A. Khosla, M. Kim, J. Nam, H. Lee, and A. Y. Ng,  “Multimodal deep learning,”  in Proc. ICML 2011, Bellevue, WA, USA, 2011.
  * [8] S. Koelstra et al.,  “DEAP: A Database for Emotion Analysis Using Physiological Signals,”  IEEE Trans. Affective Computing, vol. 3, 2011.
  * [9] T. Schäfer, P. Sedlmeier, C. Städtler, and D. Huron,  “The psychological functions of music listening,”  Frontiers in Psychology, vol. 4, pp. 511, 2013.
  * [10] S. Losorelli, D. T. Nguyen, J. Dmochowski, and B. Kaneshiro,  “NMED-T: A Tempo-Focused Dataset of Cortical and Behavioral Responses to Naturalistic Music,”  in Proc. ISMIR 2017, Suzhou, China, 2017.
  * [11] A. Vinay, A. Lerch, and G. Leslie,  “Mind the Beat: Detecting Audio Onsets from EEG Recordings of Music Listening,”  ArXiv, vol. abs/2102.06393, 2021.
  * [12] S. Stober, T. Prätzlich, and M. Müller,  “Brain Beats: Tempo Extraction from EEG Data,”  in Proc. ISMIR 2016, New York, NY, USA, 2016.
  * [13] S. Nozaradan, I. Peretz, M. Missal, and A. Mouraux,  “Tagging the Neuronal Entrainment to Beat and Meter,”  Journal of Neuroscience, pp. 10234–10240, 2011.
  * [14] R. S. Schaefer, P. Desain, and P. Suppes,  “Structural Decomposition of EEG Signatures of Melodic Processing,”  Biological Psychology, pp. 253–259, 2009.
  * [15] H. Poikonen, et al.,  “Event-related Brain Responses while Listening to Entire Pieces of Music,”  Neuroscience, vol. 312, pp. 58–73, 2016.
  * [16] E. J. Humphrey, J. P. Bello, and Y. LeCun,  “Moving Beyond Feature Design: Deep Architectures and Automatic Feature Learning in Music Informatics.,”  in Proc. ISMIR 2012, Porto, Portugal, 2012.
  * [17] M. Sarprasatham,  “Emotion Recognition: A Survey,”  Int’l Journal of Advanced Research in Computer Science, vol. 3, pp. 14–19, 01 2015.
  * [18] R. Panda, B. Rocha, and R. P. Paiva,  “Music Emotion Recognition with Standard and Melodic Audio Features,”  Applied Artificial Intelligence, pp. 313–334, 2015.
  * [19] S. Chen, Z. Gao, and S. Wang,  “Emotion recognition from peripheral physiological signals enhanced by EEG,”  in Proc. ICASSP 2016, New Orleans, LA, USA, 2016.
  * [20] X. Wang, D. Nie, and B. Lu,  “EEG-Based Emotion Recognition Using Frequency Domain Features and Support Vector Machines,”  in Proc. ICONIP 2011, Shanghai, China, 2011.
  * [21] P. C. Petrantonakis and L. J. Hadjileontiadis,  “Emotion Recognition from Brain Signals Using Hybrid Adaptive Filtering and Higher Order Crossings Analysis,”  IEEE Trans. Affective Computing, vol. 1, no. 2, pp. 81–97, 2010\.
  * [22] R. Duan, J. Zhu, and B. Lu,  “Differential Entropy Feature for EEG-based Emotion Classification,”  in Proc. Int’l IEEE/EMBS Conf. on Neural Engineering (NER), San Diego, CA, USA, 2013.
  * [23] K. Avramidis, A. Zlatintsi, C. Garoufis, and P. Maragos,  “Multiscale Fractal Analysis on EEG Signals for Music-Induced Emotion Recognition,”  in Proc. EUSIPCO 2021, Amsterdam, the Netherlands, 2021.
  * [24] Y. Wang, Z. Huang, B. McCane, and P. Neo,  “EmotioNet: A 3-D Convolutional Neural Network for EEG-based Emotion Recognition,”  in Proc. IJCNN 2018, Rio de Janeiro, Brazil, 2018.
  * [25] X. Du et al.,  “An Efficient LSTM Network for Emotion Recognition from Multichannel EEG Signals,”  IEEE Trans. Affective Computing, pp. 1–1, 2020.
  * [26] G. Andrew, R. Arora, J. Bilmes, and K. Livescu,  “Deep Canonical Correlation Analysis,”  in Proc. of Machine Learning Research, Atlanta, GA, USA, 2013, pp. 1247–1255.
  * [27] M. Li, Y. Li, S.-L. Huang, and L. Zhang,  “Semantically Supervised Maximal Correlation For Cross-Modal Retrieval,”  in Proc. ICIP 2020, Abu Dhabi, UAE, 2020.
  * [28] B. Wang, Y. Yang, X. Xu, A. Hanjalic, and H. T. Shen,  “Adversarial Cross-Modal Retrieval,”  in Proc. ACM Int’l Multimedia Conf. 2017, Mountain View, CA, USA, 2017.
  * [29] T.-K. Hu, Y.-Y. Lin, and P.-C. Hsiu,  “Learning Adaptive Hidden Layers for Mobile Gesture Recognition,”  in Proc. AAAI 2018, New Orleans, LA, USA, 2018.
  * [30] B. Li and A. Kumar,  “Query by Video: Cross-modal Music Retrieval,”  in Proc. ISMIR 2019, Delft, the Netherlands, 2019.
  * [31] L. Zhen, P. Hu, X. Wang, and D. Peng,  “Deep Supervised Cross-Modal Retrieval,”  in Proc. CVPR 2019, Salt Lake City, UT, USA, 2019.
  * [32] J. Pons, O. Nieto, M. Prockup, E. M. Schmidt, A. F. Ehmann, and X. Serra,  “End-to-End Learning for Music Audio Tagging at Scale,”  in Proc. ISMIR 2018, Paris, France, 2018.
  * [33] Y. Ganin and V. Lempitsky,  “Unsupervised Domain Adaptation by Backpropagation,”  in Proc. ICML 2015, Lille, France, 2015.
  * [34] W.-L. Zheng and B.-L. Lu,  “Investigating Critical Frequency Bands and Channels for EEG-Based Emotion Recognition with Deep Neural Networks,”  IEEE Trans. on Autonomous Mental Development, vol. 7, pp. 162–175, 2015.
  * [35] M. Won, S. Oramas, O. Nieto, F. Gouyon, and X. Serra,  “Multimodal Metric Learning for Tag-based Music Retrieval,”  arXiv preprint:2010.16030, 2020.