# Elasto-plastic large deformation analysis of multi-patch thin shells by
isogeometric approach

G. D. Huynh  X. Zhuang  H. G. Bui  G. Meschke  H. Nguyen-Xuan  Institute for
Continuum Mechanics, Leibniz University Hannover, Germany  Institute for
Structural Mechanics, Ruhr University Bochum, Germany  CIRTech Institute, Ho
Chi Minh City University of Technology (HUTECH), Ho Chi Minh City, Vietnam
Department of Architectural Engineering, Sejong University, 98 Kunja Dong,
Kwangjin Ku, Seoul, 143-747, South Korea

###### Abstract

This paper studies elasto-plastic large deformation behavior of thin shell
structures using the isogeometric computational approach with the main focus
on the efficiency in modelling the multi-patches and arbitrary material
formulations. In terms of modelling, we employ the bending strip method to
connect the patches in the structure. The incorporation of bending strips
allows to eliminate the strict demand of the $C^{1}$ continuity condition,
which is postulated in the Kirchhoff-Love theory for thin shell, and therefore
it enables us to use the standard multi-patch structure even with $C^{0}$
continuity along the patch boundaries. Furthermore, arbitrary nonlinear
material models such as hyperelasticity and finite strain plasticity are
embedded in the shell formulation, from which a unified thin shell formulation
can be achieved. In terms of analysis, the Bézier decomposition concept is
used to retain the local support of the traditional finite element. The
performance of the presented approach is verified through several numerical
benchmarks.

###### keywords:

Isogeometric analysis, Kirchhoff-Love shell theory, Multi-patch structures,
Finite Strain, Bézier decomposition.

††journal: Finite Elements in Analysis and Design

## Highlights

  * 1.

A unified thin shell formulation is established, allowing arbitrary nonlinear
material models and multi-patch shell structure to be applicable.

  * 2.

The Kirchhoff-Love shell theory is employed and the $C^{1}$ continuity at
patch boundaries is remained by the bending strip method.

  * 3.

The Bézier decomposition concept is used to retain the local support of the
traditional finite element.

  * 4.

The computational model is validated and its numerical results agree well with
ones in literature.

Nomenclatures

u |  Displacement vector
---|---
x |  Position Vector in reference configuration
X |  Position Vector in current configuration
$\bm{\xi}$ |  Curvilinear coordinates
$\bm{G}_{\alpha}$ |  Covariant base vectors in reference configuration
$\bm{G}^{\alpha}$ |  Contravariant base vectors in reference configuration
$G_{\alpha\beta}$ |  Covariant metric coefficients
$B_{\alpha\beta}$ |  Covariant curvature coefficients
$\bm{F}$ |  Deformation gradient tensor
$\bm{F}^{e}$ |  Elastic part of the deformation gradient tensor
$\bm{F}^{p}$ |  Plastic part of the deformation gradient tensor
$\bm{C}$ |  Right Cauchy deformation tensor
$\bm{C}^{p}$ |  Plastic right Cauchy deformation tensor
$\bm{b}^{e}$ |  Elastic left Cauchy deformation tensor
$\lambda_{33}$ |  thickness stretch
$\bm{E}$ |  Green-Lagrange strain tensor
$\bm{S}$ |  2nd Piola-Kirchhoff stress tensor
$\varepsilon_{\alpha\beta}$ |  Membrane strain components
$\kappa_{\alpha\beta}$ |  Bending strain components
$\Psi^{e}$ |  Strain energy density
$\Psi^{vol},\,\Psi^{dev}$ |  Volumetric and deviatoric parts of the strain energy density
$\alpha$ |  Equivalent plastic strain
$k(\alpha)$ |  Hardening function
$\Phi$ |  Yield function
$\bm{\mathbbm{C}}$ |  Material tensor
$\bm{n}$ |  Membrane force vector
$\bm{m}$ |  Bending moment vector

##  1 Introduction

Thin shell structures play a vital role in the automotive and aerospace
industry. Due to the efficiency and the accuracy of the shell formulation, it
is typically used to simulate the manufacturing processes or to determine the
weak spot of the structure during operation. To accommodate for higher
computational fidelity, large strains and rotations are used in the shell
formulation. Because the material composing shell is typically metal, the
elastoplastic constitutive model is used to investigate the failure mode, e.g.
strain localization, under large loading.

There are two typical shell models accounting for the inelastic response of
structures. The first model relies completely on stress resultants in which
the normal vector of the shell is assumed to be inextensible. The description
of the shell surface response is embodied in constitutive equations and the
spread of the plasticity through the shell thickness can not be reached
completely. Accordingly, formulations of stress resultant constitutive
relations are awkward to derive and also sophisticated to implement in the
finite element framework, for further discussion see [1]. In the standard
model, integration points through the thickness direction of the shell are
defined in association with the notion of stress-based elastoplasticity and
the distribution of plasticity through the thickness can be represented.
Therefore, the stress-based approach is more preferred in the finite element
context. The model is originally formulated from the three dimensional theory,
which can be compatible with formulations of shell kinematics in a direct or
indirect manner. Solid-shell element as presented in [2] , where only
displacement degree-of-freedom are used, can match perfectly with 3D
constitutive model without the need of ad hoc assumptions. However, this
element is usually implemented together with methods such as assumed natural
strain (ANS) and enhanced assumed strain (EAS) to alleviate locking effects.
In contrast, elements based on Kirchhoff-Love hypothesis, e.g., [3] are though
rotation-free, need to resort to plane stress enforcement to take through-
thickness behaviours into account. The implementation of three dimensional
constitutive models for finite strain plasticity with the plane stress
assumption can be classified into two common approaches: plane stress-
projected constitutive models and the nested iteration approach for the plane
stress enforcement, which are also applicable to hyperelasticity, see details
in [4] .The former approach is only applicable if involved equations of the
models are adequately simple so that transverse components can be eliminated
from the formulation, while the latter one can be implemented in a
straightforward manner for any models and is adopted in this work, see also
its applicability in recent works [5, 6].

Isogeometric analysis (IGA) emerges as a versatile tool to perform analysis
and modeling simultaneously using the same interpolation bases. IGA offers
numerical properties which are often beneficial in numerical analysis, e.g.
positiveness of the basis functions, higher continuity. IGA is shown to be
superior over standard FEM based on Lagrange shape function for dynamics
analysis . Refinement with IGA is more straightforward than with the standard
FEM counterpart, in which hpk-refinement scheme can be used to offer both
order elevation and addition of degree of freedoms (d.o.fs) conveniently
without introducing internal modes or hanging nodes. It is noted that,
refinement in IGA preserves the geometry, and thus eliminating the geometrical
error stemming from modelling. On the other hand, IGA is particularly suited
for shell modelling. The smoothness of NURBS basis functions and the
similarity between natural coordinates of the isogeometric parameterization
and the shell coordinates lead the displacement based shell formulations to be
implemented in a direct manner without the need for coordinate transformation.
In addition, NURBS is also utilized as surface modelling tool in the CAD
technology. This suggests that, extension of shell formulation for IGA is
practical, and is apparently beneficial for many industrial applications.
Nevertheless, many complex shell structures are not easy to be modelled with a
single NURBS patch. An option is to use a modelling technique that supports
connecting control point, i.e. T-splines . However, this technique is not
broadly implemented in many CAD packages. It is therefore the preferred option
to use multi NURBS patch to model a shell structure. In terms of shell
formulations, the Mindlin-Reissner theory is commonly-used in the finite
element context as only C0-continuity is required. Typical shell elements for
this theory are generations of Mixed Interpolation of Tensorial Components
(MITC) proposed by K. Bathe et al. [7, 8, 9], shell elements based on discrete
Kirchhoff-Love constraints by Areias et al. [10] [11]. These works need adhoc
techniques such as an enhanced assumed strain to alleviate locking effects
that are the most challenged hinge in the application of the finite element
method (FEM) to shell problems. For a comprehensive overview of other FEM
shell elements, readers are referred to an excellent work .In contrast to the
Mindlin-Ressner theory, shell formulations using the Kirchhoff-Love hypothesis
require $C^{1}$ continuity of the basis functions to ensure the smoothness of
bending moment interpolation. It is well known that, if the open knot vector
is used in the modelling, only $C^{0}$ continuity is achieved along the patch
boundary. Thus, additional technique, i.e. bending strip method, mortar
method, shall be used to maintain the $C^{1}$ patch continuity. While the
mortar method or Nitsche method that requires the constraints on the
derivatives of the basis function at the boundary, the bending strip method
only adds fictitious bending stiffness on the overlapping domain in the
vicinity of the boundary and its term does not complicate the weak form and
can be evaluated locally. Therefore, the latter is simple to implement, yet
ensures the reliability in the description of structure behaviours.

In this paper, an efficient computational model for thin shell using
isogeometric analysis is presented. In terms of modelling, the efficiency is
preserved by using multi NURBS patch structure with $C^{0}$ continuity on the
patch boundaries. The bending strip method [12] is employed to maintain the
$C^{1}$ continuity condition of the Kirchhoff-Love formulation. Moreover,
arbitrary 3D nonlinear materials such as hyperelasticity, finite strain
plasticity can be included in the formulation by enforcing the plane stress
condition. This leads to our main focus in this work in which a unified thin
shell formulation is established to not only model shell structures with
complex geometries or multiple patches, but also exhibit a wide range of
material behaviours. Note that, for the plane stress enforcement, we follow
the idea in [4], in which the plane stress constraint is considered an
additional equation that is used to determine the in-plane stress components
and the thickness stretch. Though similar approaches are presented in recent
works [5] [6], these methods are restricted to single-patch shell geometries
and can not be applicable to a large class of shell structures involving
multi-patches. On the aspect of analysis, the Bézier decomposition method [13]
is employed. The advantage of the Bézier decomposition approach is twofold: 1)
maintaining the local characteristic of the finite element and improving the
performance of the global stiffness matrix assembly procedure, 2) enabling the
reuse of standard finite element components, including the kinematics and the
constitutive law. To verify the applicability of the proposed computational
model, it is validated with reference examples using hyperelasticity and
elastoplasticity constitutive laws.

The structure of the paper is as following: in Section 2, geometry description
of the thin shell using Kirchhoff-Love hypothesis is presented, followed by
the shell kinematics and the treatment of plane stress enforcement which
enables the usage of three dimensional constitutive model for the shell
element. Section 3 underlines some important properties of NURBS basis
functions and give the introduction to the bending strip method with an
illustrative example. In Section 2.4, required kinematics modification to
accommodate for elastoplastic constitutive model in finite strain, which
relies on the hypothesis of multiplicative plasticity framework, is briefly
presented. Section 4 presents the discretization of the weak form based on the
directional derivative, and the formation of the Bézier finite element, based
on Bézier decomposition. Section 5 is devoted for selected numerical examples.
The paper is concluded by Section 6.

##  2 Kirchhoff-Love theory for thin shell

###  2.1 Notation

In this paper, the following notation is employed: italic letters ($g$, $G$)
are used to indicate scalar and bold letters ($\mathbf{u}$, $\mathbf{g}$,
$\mathbf{G}$) indicate the vector and second order tensor. Lower case letters
indicate terms in the reference configuration and upper case letters are for
terms in the current configuration. $\xi^{\alpha}$ ($\alpha=1,2$) denote the
convective curvilinear coordinates of the shell and $\xi^{3}$ denotes the
local coordinates in thickness direction. The indices in Greek letters
($\alpha$, $\beta$) take values of $\\{1,2\\}$ and the indices in Latin
letters ($i$, $j$) take values of $\\{1,2,3\\}$.

###  2.2 Geometry definition

The deformation and the rotation of the shell are defined on its mid-surface
as

$\mathbf{u}(\bm{\xi})=\mathbf{x}(\bm{\xi})-\mathbf{X}(\bm{\xi})$ (1)

where $\mathbf{x}$ and $\mathbf{X}$ are the position vectors of a material
point on the mid-surface in the deformed configuration and undeformed
configuration. $\bm{\xi}=\\{\xi^{1},\xi^{2},\xi^{3}\\}$ is the local
coordinates of the material point.

The covariant tangential base vectors are defined by

$\mathbf{G}_{\alpha}=\mathbf{X}_{,\alpha}=\dfrac{\partial\mathbf{X}}{\partial\bm{\xi}^{\alpha}}\qquad\mathbf{g}_{\alpha}=\mathbf{x}_{,\alpha}=\dfrac{\partial\mathbf{x}}{\partial\bm{\xi}^{\alpha}}$
(2)

Figure. 1 reveals the representation of the basis vectors between reference
and current configurations in a curvilinear coordinate system.

Figure 1: Description of body motion between different configurations using
curvilinear coordinate systems

The covariant metric coefficients of the surface are computed by

$G_{\alpha\beta}=\mathbf{G}_{\alpha}\cdot\mathbf{G}_{\beta}\qquad
g_{\alpha\beta}=\mathbf{g}_{\alpha}\cdot\mathbf{g}_{\beta}$ (3)

and the contravariant base vectors by

$\mathbf{G}^{\alpha}=G^{\alpha\beta}\mathbf{G}_{\beta}\qquad\mathbf{g}^{\alpha}=g^{\alpha\beta}\mathbf{g}_{\beta}$
(4)

with ${G^{\alpha\beta}}=({G_{\alpha\beta}})^{-1}$ and
${g^{\alpha\beta}}=({g_{\alpha\beta}})^{-1}$.

The unit normal vectors on the middle surface can be computed from the
covariant tangential base vectors as

$\mathbf{G}_{3}=\dfrac{\mathbf{G}_{1}\times\mathbf{G}_{2}}{\|\mathbf{G}_{1}\times\mathbf{G}_{2}\|}\qquad\mathbf{g}_{3}=\dfrac{\mathbf{g}_{1}\times\mathbf{g}_{2}}{\|\mathbf{g}_{1}\times\mathbf{g}_{2}\|}$
(5)

From Equations 5 and 2, one can compute the curvature tensor $\mathbf{B}$ and
$\mathbf{b}$. They are defined by the second fundamental form of surfaces

$\begin{split}B_{\alpha\beta}&=\dfrac{1}{2}\left(\mathbf{G}_{\alpha}\cdot\mathbf{G}_{3,\beta}+\mathbf{G}_{\beta}\cdot\mathbf{G}_{3,\alpha}\right)\\\
b_{\alpha\beta}&=\dfrac{1}{2}\left(\mathbf{g}_{\alpha}\cdot\mathbf{g}_{3,\beta}+\mathbf{g}_{\beta}\cdot\mathbf{g}_{3,\alpha}\right)\end{split}$
(6)

###  2.3 Kinematics

Strain measurement is represented with respect to the right Cauchy-Green
deformation tensor $\bm{C}=\bm{F}^{T}\bm{F}$, where $\bm{F}$ is the
deformation gradient

$\bm{F}=\frac{\text{d}\mathbf{x}}{\text{d}\mathbf{X}}=\mathbf{g}_{i}\otimes\mathbf{G}^{i}.$
(7)

Then, the deformation tensor can be written in terms of contravariant base
vectors as

$\bm{C}=\bm{F}^{T}\bm{F}=\bm{g}_{i}\cdot\bm{g}_{j}\bm{G}^{i}\otimes\bm{G}^{j}=g_{ij}\bm{G}^{i}\otimes\bm{G}^{j}.$
(8)

In Equation 8, components of the right Cauchy-Green deformation tensor are the
same as the metric coefficients in the deformed configuration. However,
$C_{33}$ that describes the change of the shell thickness during deformations
can not take the same value with $g_{33}$ that is equal $1$. Hence, $C_{33}$
will be determined by an additional constraint, not directly by the solution
of the midsurface variables. The thickness stretch $\lambda_{3}$ is computed
from $C_{33}$ as

$\lambda_{3}=\sqrt{C_{33}}.$ (9)

The Green-Lagrange strain tensor $\bm{E}=\frac{1}{2}(\bm{C}-\bm{I})$ is
typically separated into a constant part representing membrane action and a
linear part representing bending action:

$\displaystyle\bm{E}=E_{ij}\mathbf{G}^{i}\otimes\mathbf{G}^{j}\qquad
E_{\alpha\beta}=\epsilon_{\alpha\beta}+\xi^{3}\kappa_{\alpha\beta}$ (10)

in which

$\displaystyle\varepsilon_{\alpha\beta}$
$\displaystyle=\frac{1}{2}(g_{\alpha\beta}-G_{\alpha\beta})$ (11)
$\displaystyle\kappa_{\alpha\beta}$
$\displaystyle=b_{\alpha\beta}-B_{\alpha\beta}$ (12)

In order to be consistent with material models which are often expressed in
the Cartesian coordinate frame, the strain quantities and the deformation
gradient need to be represented in a local Cartesian coordinate system. A
general expression for such a transformation is given by

$\bar{(\cdot)}_{\gamma\delta}=(\cdot)_{\gamma\delta}(\bm{E}_{\gamma}\cdot\bm{G}^{\alpha})(\bm{G}^{\beta}\cdot\bm{E}_{\delta}),$
(13)

where $\bar{(\cdot)}_{\gamma\delta}$ and $(\cdot)_{\gamma\delta}$ are
components of generic second order tensors obtained in local Cartesian and
curvilinear coordinates respectively. The local Cartesian basis vectors
$\bm{E}_{\gamma}$ and $\bm{E}_{\delta}$ are given by equations

$\bm{E}_{1}=\frac{\bm{G}_{1}}{||\bm{G}_{1}||},\,\,\,\,\,\,\bm{E}_{2}=\frac{\bm{G}_{2}-(\bm{G}_{2}\cdot\bm{E}_{1})\bm{E}_{1}}{||\bm{G}_{2}-(\bm{G}_{2}\cdot\bm{E}_{1})\bm{E}_{1}||},\,\,\,\,\,\,\bm{E}_{3}=\bm{G}_{3}$
(14)

###  2.4 Finite strain kinematics for elastoplasticity

The kinematics presented in Section 2.3 can be applied for the constitutive
laws exhibiting elastic behaviours, such as linear elasticity or
hyperelasticity. For constitutive laws accounting for plasticity effects, the
elastoplastic kinematics model for finite strain [14, 15] shall be employed,
which is based on the concept of the plastically deforming intermediate
configuration and the multiplicative decomposition of the deformation
gradient.

To account for the plasticity behaviour, the deformation gradient is
decomposed as

$\bm{F}=\bm{F}^{e}\bm{F}^{p}$ (15)

with $\bm{F}^{e}$ and $\bm{F}^{p}$ as the elastic and plastic parts of the
deformation gradient tensor.

Following Equation 15, the plastic right Cauchy-Green deformation tensor is
computed as

$\bm{C}^{p}=\left(\bm{F}^{p}\right)^{T}\bm{F}^{p}$ (16)

and the elastic left Cauchy-Green deformation tensor:

$\bm{b}^{e}=\bm{F}^{e}\left(\bm{F}^{e}\right)^{T}$ (17)

From Equations 15, 16 and 17, the elastic left deformation tensor
$\mathbf{b}^{e}$ can be computed from $\mathbf{C}^{p}$ and $\mathbf{F}$:

$\bm{b}^{e}=\bm{F}\left(\bm{C}^{p}\right)^{-1}\bm{F}^{T}$ (18)

The total rate of the elastic deformation tensor is obtained by Lie derivative
as follows

$\mathcal{L}[\bm{b}^{e}]=\bm{F}\left(\dot{\bm{C}}^{p}\right)^{-1}\bm{F}^{T}$
(19)

The elastic behaviour of the material is modeled by strain energy density
function of hyperelastic material, which is further separated into the
deviatoric and volumetric components as

$\Psi^{e}=\Psi^{e}_{vol}+\Psi^{e}_{dev}$ (20)

with

$\displaystyle\Psi^{e}_{vol}$
$\displaystyle=\frac{1}{2}K\left(\frac{1}{2}(\left(J^{e}\right)^{2}-1)-\text{ln}J^{e}\right)$
(21) $\displaystyle\Psi^{e}_{dev}$
$\displaystyle=\frac{1}{2}G\left(\text{tr}[\hat{\bm{b}}^{e}]-3\right)$ (22)

where $J^{e}=\text{det}[\bm{F}^{e}]$,
$\hat{\bm{b}}^{e}=\left(J^{e}\right)^{-2/3}\bm{b}^{e}$, and $G$ and $K$ are
identified as the shear modulus and the bulk modulus in infinitesimal
deformation respectively.

The classical Mises-Huber yield condition established with respect to
deviatoric part of the Kirchhoff stress $\bm{\tau}_{dev}$ and the hardening
function $k(\alpha)$ is given by

$\Phi(\bm{\tau}_{dev},\alpha)=||\bm{\tau}_{dev}||-\sqrt{\frac{2}{3}}k(\alpha)\leq
0$ (23)

where $\Phi$ is the yield function and $\alpha$ the equivalent plastic strain.
Based on the notion of the principle of maximum plastic dissipation, the
associative flow rule can be expressed by

$\left(\dot{\bm{C}}^{p}\right)^{-1}=-\frac{2}{3}\gamma\,\text{tr}[\bm{b}^{e}]\,\bm{F}^{-1}\bm{n}\bm{F}^{-T}$
(24)

with $\bm{n}=\bm{\tau}_{dev}/||\bm{\tau}_{dev}||$ and $\gamma$ is the plastic
multiplier. The evolution of the yield stress known as hardening phenomenon is
driven by the equation

$\dot{\alpha}=\sqrt{\frac{2}{3}}\gamma$ (25)

Finally, the so-called loading/unloading conditions of the elastoplastic model
are governed by the Karush-Kuhn-Tucker condition:

$\gamma\geq 0,\,\,\,\,{\Phi}(\bm{\tau}_{dev},\alpha)\leq
0,\,\,\,\gamma{\Phi}(\bm{\tau}_{dev},\alpha)=0$ (26)

accompanied by the consistency condition:

$\gamma\dot{{\Phi}}(\bm{\tau}_{dev},\alpha)=0.$ (27)

####  2.4.1 Time integration for the finite strain plasticity model

The solution of Equations 24 and 25 requires a time integration scheme, for
which the backward Euler scheme is adopted, such that

$\displaystyle\left(\bm{C}_{n+1}^{p}\right)^{-1}$
$\displaystyle=\left(\bm{C}_{n}^{p}\right)^{-1}-\frac{2}{3}\Delta\gamma[\bm{b}^{e}_{n+1}]\bm{F}_{n+1}^{-1}\bm{n}_{n+1}\bm{F}_{n+1}^{-T}$
(28) $\displaystyle\alpha_{n+1}$
$\displaystyle=\alpha_{n}+\sqrt{\frac{2}{3}}\Delta\gamma$ (29)

$\Delta\gamma$ is the incremental plastic multiplier and satisfies
$\Delta\gamma>0$.

The classical return mapping scheme [16] [14] is utilized to integrate
Equation 29, and summarized in Table.1.

Table 1: Hyperelastic extension of return-mapping algorithm . $\bm{F}_{n+1}$
are given at $t_{n+1}$, $\bm{C}_{n}^{p^{-1}}$ and $\alpha_{n}$ are known at
$t_{n}$.
---
Elastic predictor step stemming from intermediate configuration:
$\bar{\bm{F}}_{n+1}=J_{n+1}^{-1/3}\,\bm{F}_{n+1}$
$J_{n+1}=\text{det}[{\bm{F}}_{n+1}]$
$\bar{\bm{b}}_{n+1}^{e,trial}=\bar{\bm{F}}_{n+1}\,\bm{C}_{n}^{p-1}\,\bar{\bm{F}}_{n+1}^{T}$
$\bm{\tau}_{dev,n+1}^{trial}=\text{G}\,\text{dev}[\bar{\bm{b}}_{n+1}^{e,trial}]$
$\bm{\tau}_{vol,n+1}^{trial}=\frac{\text{1}}{2}K\,(J^{2}_{n+1}-1)\,\bm{I}$
Check for yielding:
IF
$\Phi_{n+1}^{trial}=||\bm{\tau}_{dev,n+1}^{trial}||-\sqrt{\frac{2}{3}}R(\alpha_{n})\leq
0$
THEN set $(\cdot)_{n+1}=(\cdot)_{n+1}^{trial}$ and EXIT
ELSE go to corrector step
Return mapping corrector:
Compute $\Delta\gamma$ by solving the equation
$\Phi_{n+1}^{trial}=||\bm{\tau}_{dev,n+1}^{trial}||-\sqrt{\frac{2}{3}}k(\alpha_{n}+\sqrt{\frac{2}{3}}\Delta\gamma)-\frac{2}{3}\Delta\gamma\text{G}\text{tr}[\bar{\bm{b}}_{n+1}^{e,trial}]=0$
Update the state variables
$\bm{\tau}_{dev,n+1}=\bm{\tau}_{dev,n+1}^{trial}-\frac{2}{3}\Delta\gamma\text{G}\text{tr}[\bar{\bm{b}}_{n+1}^{e,trial}]\bm{n}_{n+1}^{trial}$
$\bm{\tau}_{vol,n+1}=\bm{\tau}_{vol,n+1}^{trial}$
$\bm{\sigma}_{n+1}=\frac{1}{J}(\bm{\tau}_{vol,n+1}+\bm{\tau}_{dev,n+1})$
$\bm{C}_{n+1}^{p^{-1}}=\bm{C}_{n}^{p^{-1}}-\frac{2}{3}\Delta\gamma\text{tr}[\bar{\bm{b}}_{n+1}^{e,trial}]\bm{F}_{n+1}^{-1}\bm{n}_{n+1}^{trial}\bm{F}_{n+1}^{-T}$
$\alpha_{n+1}=\alpha_{n}+\sqrt{\frac{2}{3}}\Delta\gamma$

###  2.5 Zero transverse normal stress enforcement for thin shell

Here, the indigenous three-dimensional constitutive models are utilized within
a Newton-Raphson loop by solving the equation of plane stress constraint at
each Gauss point with thickness stretch $C_{33}$ as the unknown.

Following [4] [17], $C_{33}$ is obtained by satisfying the zero transverse
normal stress condition, which takes the form of

$S^{33}=\mathbbm{C}^{33\alpha\beta}E_{\alpha\beta}+\mathbbm{C}^{3333}E_{33}=0$
(30)

from which the coefficients of the material tensor are modified as

$\hat{\mathbbm{C}}=\mathbbm{C}-\frac{\mathbbm{C}^{\alpha\beta
33}\mathbbm{C}^{33\gamma\delta}}{\mathbbm{C}^{3333}}$ (31)

$C_{33}$ is solved iteratively by the equation

$S^{33}+\frac{\partial S^{33}}{\partial C_{33}}\Delta
C_{33}=S^{33}+\frac{1}{2}\mathbbm{C}^{3333}\Delta C_{33}=0$ (32)

followed by an incremental update

$\Delta C_{33}^{(I)}=-2\frac{S^{33,(I)}}{C^{3333,(I)}}$ (33)

where $I$ denotes the number of iteration step.

The thickness change is corrected by computing $C_{33}$ and $\lambda_{3}$ as
follows

$\displaystyle C_{33}^{(I+1)}$ $\displaystyle=C_{33}^{(I)}+\Delta
C_{33}^{(I)}$ (34) $\displaystyle\lambda_{3}^{(I+1)}$
$\displaystyle=\sqrt{C_{33}^{(I+1)}}$ (35)

$\lambda_{3}$ will be updated until $S^{33}$ converges to a defined tolerance
which fulfills the zero transverse normal stress condition. Subsequently, the
statically material tensor $\hat{\mathbbm{C}}$ is obtained and only its in-
plane components are chosen for the shell structures.

###  2.6 Variational form of the equilibrium equation

The mechanical response of a structure is governed by the principle of virtual
work which requires the sum of internal and external work vanishes in the
equilibrium state

$\delta W=\delta W^{int}-\delta W^{ext}$ (36)

The internal virtual work is defined by

$\displaystyle\delta W^{int}$ $\displaystyle=\int\limits_{\Omega}S_{ij}:\delta
E_{ij}d\,\Omega=\int\limits_{A}\Big{(}n_{\alpha\beta}:\delta\varepsilon_{\alpha\beta}+m_{\alpha\beta}:\delta\kappa_{\alpha\beta}\Big{)}\,dA,$
(37)

and the external virtual work takes the form of

$\delta W^{ext}=\int\limits_{A}\delta u_{i}f_{i}\,dA.$ (38)

##  3 NURBS-based isogeometric analysis for shell structures

In this section, the useful properties of NURBS basis functions which are
suitable for the analysis of Kirchhoff-Love shells are summarized, followed by
an introduction to the bending strip method which is employed in the multi-
patch case.

###  3.1 Interpolation using NURBS basis functions

Due to the appearance of the second derivatives in curvature changes of the
shell kinematics, $C^{1}$-continuity of the employed basis functions must be
required to attain a conforming discretization of the Kirchhoff-Love
hypothesis. This can be simply achieved using NURBS basis functions which are
used to approximate the solution fields and represent the shell geometry
exactly. Following that, given a control net $\bm{P}_{i,j},1\leq i\leq n,1\leq
j\leq m$, polynomial orders $p$ and $q$, and knot vectors
$\Xi_{1}=\\{\xi_{1}^{1},\xi_{1}^{2},...,\xi_{1}^{n+p+1}\\}$ and
$\Xi_{2}=\\{\xi_{2}^{1},\xi_{2}^{2},...,\xi_{2}^{m+q+1}\\}$, a shell surface
can be represented by sum:

$\bm{S}(\xi_{1},\xi_{2})=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}N_{i,j}(\xi_{1},\xi_{2})\bm{P}_{i,j}$
(39)

where $\xi_{1}\in\Xi_{1}$, $\xi_{2}\in\Xi_{2}$ are coordinates in the
parametric space, $N_{i,j}$ are the NURBS basis functions, which are defined
as

$N_{i,j}(\xi_{1},\xi_{2})=\frac{B_{i,p}(\xi_{1})B_{j,q}(\xi_{2})w_{i,j}}{\sum\limits_{k=1}^{n}\sum\limits_{l=1}^{m}B_{k,p}(\xi_{1})B_{l,q}(\xi_{2})w_{k,l}}.$
(40)

In Equation 40, $B$ is the univariate B-spline basis functions and $w$ is the
weight factor. In the parametric space, the elements are defined by
partitioning the knots into knot spans. Inside each knot span, the NURBS basis
functions are $C^{\infty}$-continuous and $C^{p-k}$ at each knot ($k$ denotes
the multiplicity of the knot). When $k\leq p-1$, the $C^{1}$ continuity is
satisfied and is useful for the analysis of Kirchhoff-Love thin shell.
Nevertheless, in the typical scenario, the open knot vector is used [18],
which leads to only $C^{0}$ continuity at the patch boundaries. Without
additional treatment, a kink in bending moment will occur on the patch
boundaries if the Kirchhoff-Love formulation is used. This issue will be
further described in the next section.

Displacements of the shell are approximated with respect to the NURBS basis
functions which are used to model the midsurface of the shell as

$\bm{u}=N_{I}(\xi_{1},\xi_{2})\bm{u}_{I}$ (41)

where $\bm{u}_{I}$ are the displacement at the I-th control point. Following
that, the first and second derivatives of the displacements can be obtained:

$\displaystyle\bm{u}_{,\alpha}$ $\displaystyle=\dfrac{\partial
N_{I}(\xi_{1},\xi_{2})}{\partial\xi_{\alpha}}\bm{u}_{I}$ (42)
$\displaystyle\bm{u}_{,\alpha\beta}$
$\displaystyle=\dfrac{\partial^{2}N_{I}(\xi_{1},\xi_{2})}{\partial\xi_{\alpha}\partial\xi_{\beta}}\bm{u}_{I}$
(43)

###  3.2 The bending strip method

The bending strip method [12] is a mesh connecting method, which fundamentally
adds additional kinematics constraints at patch interface in order that the
bending moment can be transmitted between patches. The strip is constructed by
the triples of control points, one at the shared boundary, one on each side of
two involved patches. With this construction, the parametric domain of the
bending strip surface comprises of one quadratic element in the direction
perpendicular to the interface and linear elements along the boundary. For an
illustration of the bending strip method, the reader is referred to Figure 2
(Left). The principle of virtual work in Equation 36 is rewritten to account
for the bending strips as

$\delta W=\delta W^{int}+\delta W^{int}_{strip}-\delta W^{ext}$ (44)

The new term $\delta W^{int}_{strip}$ is computed as

$\delta
W^{int}_{strip}=\int\limits_{A}\frac{t^{3}}{12}\bm{\kappa}^{T}\bm{\mathbbm{C}}_{strip}\bm{\kappa}\,dA.$
(45)

Figure 2: Left: Illustration of the bending strip domain; Right: Geometry of
the Hex-can example

Adopting the bending material stiffness proposed in [12], we have

$\displaystyle\bm{\mathbbm{C}}_{strip}=\begin{bmatrix}0&0&0\\\
0&E_{strip}&0\\\ 0&0&0\end{bmatrix}$ (46)

in which the value of the bending strip modulus $E_{strip}$ must be chosen to
be high enough and a range of its values varies from $10^{3}\times E$ to
$10^{5}\times E$ with $E$ as Young’s modulus.

---

Figure 3: Hex-can subjected to internal pressure. Deformed configuration (a)
without bending strips (b) with bending strips

.

Here, we illustrate the method by an example representing by a Hex-can
geometry subjected to internal pressure (see Figure 2 (Right) for the geometry
of the Hex-can). The Hex-can geometry comprises of six patches joining at
patch boundaries. Due to the flat geometry of each patch, kinks occur at the
patch interfaces. With this configuration, the bending moment is not
transmitted if C0 continuity between patches is enforced. When bending strips
are added to couple the triples of control points at interfaces of each patch,
the kinks do not arise (see Figure 3 (Left) for the displacement without
bending strips and Figure 3 (Right) for the displacement with bending strips).
As can be observed in Figure 3, kinks do not arise, which means that bending
moment is transferred between patches.

##  4 FE discretization

###  4.1 Discretization of the weak form

To solve Equation 36 using the Newton-Raphson scheme. It shall be linearized
with respect to $\Delta\mathbf{u}$:

$\delta W(\bm{u},\delta\bm{u})+D\delta
W(\bm{u},\delta\bm{u})[\Delta\bm{u}]=\bm{0}$ (47)

which can be simply rewritten in a discretized form as

$\delta v_{a}K_{ab}\Delta u_{b}=\delta v_{a}R_{a}$ (48)

In Equation 48, $a$ and $b$ denote the global degree of freedom (d.o.f) number
of the displacement field. $R_{a}$ and $K_{ab}$ are the derivatives of $W$ and
$R_{a}$, respectively:

$R_{a}=\frac{\partial W}{\partial u_{a}}\qquad\qquad K_{ab}=\frac{\partial
R_{a}}{\partial u_{b}}=\frac{\partial^{2}W}{\partial u_{a}\partial u_{b}}$
(49)

The residual forces vector is the difference of the internal forces vector and
the external forces vector:

$R_{a}=R^{int}_{a}-R^{ext}_{a},$ (50)

in which

$\displaystyle R^{int}_{a}$
$\displaystyle=\int\limits_{A}\left(n_{\alpha\beta}\frac{\partial\varepsilon_{\alpha\beta}}{\partial
u_{a}}+m_{\alpha\beta}\frac{\partial\kappa_{\alpha\beta}}{\partial
u_{a}}\right)\,dA$ (51) $\displaystyle R^{ext}_{a}$
$\displaystyle=\int\limits_{A}\frac{\partial u_{i}}{\partial
u_{a}}t_{i}\,dA=\int\limits_{A}N_{I}t_{i}\,dA$ (52)

The tensors $n_{\alpha\beta}$ and $m_{\alpha\beta}$ in Equation 51 denote the
stress resultants, which are computed by integrating the second Piola-
Kircchoff stress tensor $S^{\alpha\beta}$ over the thickness direction
$\xi^{3}$

$n_{\alpha\beta}=\int\limits_{-t/2}^{t/2}S^{\alpha\beta}\,d\xi^{3}\qquad
m_{\alpha\beta}=\int\limits_{-t/2}^{t/2}S^{\alpha\beta}\,\xi^{3}\,d\xi^{3}$
(53)

The tangential stiffness matrix in Equation 49 can be split into material and
geometric parts as

$K_{ab}=K_{ab}^{geo}+K_{ab}^{mat}$ (54)

in which

$\displaystyle K_{ab}^{mat}$
$\displaystyle=\int\limits_{A}\left(\frac{\partial n_{\alpha\beta}}{\partial
u_{b}}\frac{\partial\varepsilon_{\alpha\beta}}{\partial u_{a}}+\frac{\partial
m_{\alpha\beta}}{\partial u_{b}}\frac{\partial\kappa_{\alpha\beta}}{\partial
u_{a}}\right)dA$ (55) $\displaystyle K_{ab}^{geo}$
$\displaystyle=\int\limits_{A}\left(n_{\alpha\beta}\frac{\partial^{2}\varepsilon_{\alpha\beta}}{\partial
u_{a}\partial
u_{b}}+m_{\alpha\beta}\frac{\partial^{2}\kappa_{\alpha\beta}}{\partial
u_{a}\partial u_{b}}\right)dA$ (56)

It is noted that, the case of displacement-dependent load
$t_{i}=t_{i}(u_{i})$, which leads to the dependency of tangential stiffness to
the internal forces, is not considered. The derivatives of the stress
resultants w.r.t nodal displacements take the form of

$\displaystyle\frac{\partial n_{\alpha\beta}}{\partial u_{b}}$
$\displaystyle=\left(\int\limits_{-t/2}^{t/2}\hat{\mathbbm{C}}^{\alpha\beta\gamma\delta}d\xi^{3}\right)\frac{\partial\varepsilon_{\gamma\delta}}{\partial
u_{b}}+\left(\int\limits_{-t/2}^{t/2}\hat{\mathbbm{C}}^{\alpha\beta\gamma\delta}\xi^{3}d\xi^{3}\right)\frac{\partial\kappa_{\gamma\delta}}{\partial
u_{b}}$ (57) $\displaystyle\frac{\partial m_{\alpha\beta}}{\partial u_{b}}$
$\displaystyle=\left(\int\limits_{-t/2}^{t/2}\hat{\mathbbm{C}}^{\alpha\beta\gamma\delta}\xi^{3}d\xi^{3}\right)\frac{\partial\varepsilon_{\gamma\delta}}{\partial
u_{b}}+\left(\int\limits_{-t/2}^{t/2}\hat{\mathbbm{C}}^{\alpha\beta\gamma\delta}(\xi^{3})^{2}d\xi^{3}\right)\frac{\partial\kappa_{\gamma\delta}}{\partial
u_{b}}$ (58)

For more details on the derivatives of strain tensor with respect to nodal
displacement, the reader is referred to [5, 19]. Note that the four-order
tensor $\hat{\mathbbm{C}}^{\alpha\beta\gamma\delta}$ in the equations above is
expressed in Appendix B which is then modified following Eq.(31).

###  4.2 Bézier decomposition

To evaluate basis function Equation 40, the full knot vector is required. In
addition, to interpolate the displacement within an element, the displacement
at other control points in the patch is necessary. To alleviate this issue and
maintain the local characteristic of the finite element, the Bézier
decomposition concept is proposed, which bridges the gap between isogeometric
method and standard finite element method [13]. Taking into account that not
all B-splines basis function is nonzero on a knot span, the basis functions
$\mathbf{N}$, which is not vanished on a knot span, can be represented by a
linear combination of Bézier basis functions $\mathbf{B}$ constructed on knot
vector
$\\{\underbrace{0,\ldots,0}_{p+1\,\text{times}},\underbrace{1,\ldots,1}_{p+1\,\text{times}}\\}$,
in short:

$\mathbf{N}=\mathcal{C}\mathbf{B}$ (59)

Figure 4: Bézier elements over the isogeometric control mesh;  : control
point;  : control point supported by Bézier element.

The coefficient matrix $\mathcal{C}$ is called Bézier extraction operator and
depends only on the knot vectors, hence it is considered as constant. For more
details on how to compute $\mathcal{C}$, the reader is referred to [13].

Based on the Bézier decomposition concept, the Bézier finite element can be
constructed, which contains only the control points associated with non-
vanished basis function on the knot span. An illustration of the Bézier finite
element for 2D case can be found in Figure 4.

To evaluate the second derivatives terms, i.e. Equation 56, the second
derivatives of the basis functions are required. Appendix A presents a method
to take the second derivatives of NURBS basis functions evaluated through
Bézier extraction operator and Bézier basis functions.

##  5 Numerical examples

In this section, four selected examples are performed to analyze hyperelastic
and elastoplastic behaviours of multi-patch thin shell structures. Comparisons
with the reference results are made. Arc-length control method following the
approach in [20] is adopted to capture the structural response in the last two
examples.

###  5.1 Numerical example 1: Pinched semi-cylindrical shell

For the first example, a semi-cylidrical shell is pinched by an end force at
the middle of one end and fully clamped in the other direction. The point load
is applied incrementally until the maximum load $P=2000.0$ is reached. Figure.
5 shows the problem setup and boundary conditions. The length of the cylinder
is $L=3.048$ and the radius is $R=1.016$. St. Vennant-Kirchhoff material model
is used with $E=2.0685\times 10^{7}$ and $\rho=0.3$.

Figure.6 reveals obtained results on four different mesh sizes with cubic
NURBS elements, indicating that the computational procedure fails at $P=500.0$
with $8\times 8$ elements and accurate results are obtained with the mesh
sizes of $12\times 12$, $16\times 16$ and $20\times 20$. The deformed geometry
of the cylinder under the maximum load is revealed in Fig.7. Table.2 reports
the total number of iterations required for the computation, in which the
number for our approach that requires 3 to 9 iterations at each load step is
less than the one for other shell elements [21] [22] using finite element
methods .

Figure 5: Numerical example 1: Geometry and boundary condition.  Figure 6: Numerical example 1: Investigation of mesh sensitivity with respect to the load-deflection curve.  Figure 7: Numerical example 1: Deformed geometry at the maximum load.  | 16x16 | 32x32 S4R elements | 16x16 | 16x16 elements
---|---|---|---|---
| current elements | [21] | current elements |  ,B. Brank et al. [22]
Incremental load | 40 | 33 | 80 | 80
Number of iterations | 124 | 184 | 76 | 110
Table 2: Numerical example 1: Lists of the numbers of iterations and load
increments.

###  5.2 Numerical example 2: Pinching of cylinder

In the second example, the proposed computational model is validated against
finite strain hyperelasticity model. Following that, a tube with geometry and
boundary conditions as illustrated in Figure 13 is analyzed. The geometric
parameters are chosen as length $L=30cm$, radius $R=9cm$, and thickness
$t=0.2cm$. The bottom of the cylinder is fixed and a uniformly distributed
line load ($p$) is applied on top (see Figure 8). The constitutive law is Neo-
Hookean with the strain energy is defined in Equation 60. The material
parameters are $\mu=60kN/mm^{2}$ and $\lambda=240kN/mm^{2}$. The geometry of
the cylinder consists of 4 patches, and is discretized using 552 cubic Bézier
elements.

$W=\dfrac{\mu}{2}\left(\text{tr}\bm{C}-3\
\right)-\mu\ln{\sqrt{\det{\bm{C}}}}+\dfrac{\lambda}{4}\left(\det{\bm{C}}-1-2\ln{\sqrt{\det{\bm{C}}}}\right)$
(60)

Figure 8: Numerical example 2: Geometry and boundary condition

The load $p$ is increased equally in 8 steps. Table 3 lists computed
displacements at point A corresponding to applied load at each step. From
Figure 9 (Right), the structural load when the vertical deflection is
$u=160mm$ can be approximated as $F\approx 34.965kN$, which is in good
agreement with the results from literature ranging between $34.59kN$ and
$35.47kN$ [23]. A contour plot of the last deformed configuration is depicted
in Figure 9 (Left).

Steps | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---
$F(kN)$ | $4.5$ | $9.0$ | $13.5$ | $18.0$ | $22.5$ | $27.0$ | $31.5$ | $36.0$
$u(mm)$ | $13.817$ | $30.872$ | $51.954$ | $76.473$ | $98.818$ | $121.337$ | $143.815$ | $164.83$
Table 3: Numerical example 2: Resulting displacements corresponding to applied
load at each step.

Figure 9: Numerical example 2: Left) Contour plot of the last deformed
configuration; Right) Load-displacement curve measured on top.

###  5.3 Numerical example 3: Scordelis-Lo roof

In the next example, the failure analysis of the Scordelis-Lo roof, accounting
for fully non-linear thin shell analysis, is considered. Figure 10 illustrates
boundary conditions and the setup of the geometry of the problem which
contains 2 patches connected at the middle line on the roof top.The material
parameters are chosen as elastic modulus $E=2.1\times 10^{4}N/mm^{2}$,
Poisson’s ratio $\rho=0$ and the hardening function
$\kappa(\alpha)=4.2N/mm^{2}$. Because the hardening modulus is vanished, the
perfect plasticity condition is obtained. The structure is applied with
gravity load. A mesh comprises of $3969$ quadratic Bézier elements is used to
discretize the whole domain. Arc-length control method is employed to capture
the structural response with the reference value $f_{0}=4\times
10^{-3}N/mm^{2}$.

Figure 10: Numerical example 3: Geometry and boundary conditions. It is noted
that, the rigid diaphragms are used to support the two ends of the roof.

Figure 11 shows deformed configurations at different loading stages. As can be
observed from Figure 11, the localized failure mode is characterized by the
appearance of the plastic hinge along the axial line on the roof top.

The corresponding load-deflection curve at center point of the side (denoted
A) is plotted in Figure 12. It can be observed that the obtained results agree
well with the reference[6, 24].

|
---|---
|
Figure 11: Numerical example 3: Deformation at different stages.
$0$$500$$1{,}000$$1{,}500$$2{,}000$$2{,}500$$0$$0.5$$1$$1.5$u
(mm)$f/f_{0}$Reference [24]Reference [6]Obtained results Figure 12: Numerical
example 3: Load-deflection curve at point A.

###  5.4 Numerical example 4: Pinching of cylinder with large strains

The last example considers a cylinder supported by two rigid diaphragms and
pinched by two concentrated forces at the middle of the opposite sides. Figure
13 presents the setup of boundary conditions and the geometry of the cylinder
which contains four patches connected at common interfaces A1, A2, A3 and A4.
The material parameters are chosen as $E=3000$, $\rho=0.3$ and the hardening
function $\kappa(\alpha)=24.3+300\alpha$. The geometry contains 8 patches as
depicted in Figure 14 (Top Left). The bending strips are used on each patch
interfaces to enforce the $C^{1}$-continuity condition. A mesh contains 3452
quadratic Bézier elements is used to discretize the whole domain.

Figure 13: Numerical example 4: Geometry and boundary conditions. (A1, A2, A3
and A4 are marked in red.)

Deformed configurations at different stages of loading are illustrated in
Figure 14. The load-deflection curve representing the structural response,
obtained at the loading point, is shown in Figure 15. The obtained results are
in good agreement with the reference [6, 25].

|
---|---
|
Figure 14: Numerical example 4: Deformation at different stages.
$0$$0$$10$$20$$30$u (mm)F (kN)Reference [25]Reference [6]Obtained results
Figure 15: Numerical example 4: Load-deflection curve.

##  6 Conclusion

An efficient computational model for thin shell using isogeometric analysis is
proposed, with the advantage of using bending strip method to maintain the
$C^{1}$ continuity, thus allowing for efficient modelling using multiple NURBS
patch with $C^{0}$ continuity at the patch boundaries. The Bézier
decomposition method is used to retain the local characteristics of the finite
element. This computational model is validated with respect to nonlinear
constitutive model, including hyperelasticity and elastoplasticity. The
benchmark results agree well with the references, showing that the proposed
computational model is promising for efficient shell analysis. The presented
shell formulation is very general and flexible and can be extended to dynamics
and ductile fracture motivated by published papers that are classified into
two categories: discrete crack modeling [11] [26] [27] and smeared crack
modeling [28] [29] . On the other hand, in order to further advance the
efficiency, the parallelization for distributed computing shall be done. The
parallel algorithm shall take into account the distribution of bending strips
over computing processes. That will be addressed in future research.

## Acknowledgements

The first author would like to acknowledge the financial support via RISE
project BESTOFRAC 734370 for this work. The research performed by Hoang-Giang
Bui and Günther Meschke were conducted in the framework of the Collaborative
Research Project SFB 837 ”Interaction Modeling in Mechanized Tunneling”,
financed by the German Research Foundation (DFG). The authors would like to
thank the DFG for the support of this project.

##  Appendix A Second derivatives of NURBS basis functions based on Bézier
decomposition

The NURBS basis function, which is non-vanished on the knot span, is defined
by

$\mathbf{R}=\dfrac{\mathcal{W}\mathcal{C}\mathbf{B}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle}$
(61)

In Equation 61, the operator $\langle\cdot,\cdot\rangle$ denotes the dot
product of 2 vectors. $\mathcal{W}$ is the diagonal matrix with entries as the
corresponding weight of the non-vanished basis function.
$\mathbf{W}=\text{diag}\\{\mathcal{W}\\}$ and
$\mathbf{W}^{b}=\mathcal{C}^{T}\mathbf{W}$.

The first derivatives of Equation 61 can be computed as

$\dfrac{\partial\mathbf{R}}{\partial\xi^{\alpha}}=\mathcal{W}\mathcal{C}\left[\dfrac{\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\alpha}}}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle}-\dfrac{\langle\mathbf{W}^{b},\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\alpha}}}\rangle\mathbf{B}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle^{2}}\right]$
(62)

Taking the derivatives of Equation 62 with respect to $\xi^{\alpha}$ and
$\xi^{\beta}$ gives

$\dfrac{\partial^{2}\mathbf{R}}{\partial(\xi^{\alpha})^{2}}=\mathcal{W}\mathcal{C}\left[\dfrac{\nicefrac{{\partial^{2}\mathbf{B}}}{{\partial(\xi^{\alpha})^{2}}}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle}-2\dfrac{\langle\mathbf{W}^{b},\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\alpha}}}\rangle\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\alpha}}}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle^{2}}+2\dfrac{\langle\mathbf{W}^{b},\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\alpha}}}^{2}\rangle\mathbf{B}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle^{3}}-\dfrac{\langle\mathbf{W}^{b},\nicefrac{{\partial^{2}\mathbf{B}}}{{\partial(\xi^{\alpha})^{2}}}\rangle\mathbf{B}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle^{2}}\right]$
(63)

And

$\begin{split}\dfrac{\partial^{2}\mathbf{R}}{\partial\xi^{\alpha}\partial\xi^{\beta}}=\mathcal{W}\mathcal{C}\bigg{[}&\dfrac{\nicefrac{{\partial^{2}\mathbf{B}}}{{\partial\xi^{\alpha}\partial\xi^{\beta}}}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle}-\dfrac{\langle\mathbf{W}^{b},\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\alpha}}}\rangle\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\beta}}}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle^{2}}-\dfrac{\langle\mathbf{W}^{b},\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\beta}}}\rangle\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\alpha}}}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle^{2}}\\\
&+2\dfrac{\langle\mathbf{W}^{b},\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\alpha}}}\rangle\langle\mathbf{W}^{b},\nicefrac{{\partial\mathbf{B}}}{{\partial\xi^{\beta}}}\rangle\mathbf{B}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle^{3}}-\dfrac{\langle\mathbf{W}^{b},\nicefrac{{\partial^{2}\mathbf{B}}}{{\partial\xi^{\alpha}\partial\xi^{\beta}}}\rangle\mathbf{B}}{\langle\mathbf{W}^{b},\mathbf{B}\rangle^{2}}\bigg{]}\end{split}$
(64)

##  Appendix B Consistent elasto-plastic tangent moduli

Scaling factors:
---
$\bar{G}=\frac{1}{3}G\text{tr}(\bar{\bm{b}}^{e})$
$\Xi_{0}=1+\frac{R^{\prime}}{3\bar{G}}$
$\Xi_{1}=[1-\frac{1}{\Xi_{0}}]\frac{2}{3}\frac{||\bm{\tau}_{dev,n+1}^{trial}||}{\bar{G}}\Delta\gamma$
$\Xi_{2}=\frac{2\bar{G}\Delta\gamma}{||\bm{\tau}_{dev,n+1}^{trial}||}$
$\Xi_{3}=\frac{1}{\Xi_{0}}-\Xi_{2}+\Xi_{1}$
$\Xi_{4}=[\frac{1}{\Xi_{0}}-\Xi_{3}]\frac{||\bm{\tau}_{dev,n+1}^{trial}||}{\bar{G}}$
Spatial elasticity tensor for hyperelasticity:
$\bm{\mathbbm{c}}_{n+1}^{e,trial}=\frac{1}{J}[\bm{\mathbbm{c}}_{vol,n+1}^{e,trial}+\bm{\mathbbm{c}}_{dev,n+1}^{e,trial}]$
$\bm{\mathbbm{c}}_{dev,n+1}^{e,trial}=2\bar{G}(\bm{I}-\frac{1}{3}\bm{I}\otimes\bm{I})-\frac{2}{3}[(\bm{\tau}_{dev,n+1}^{trial}\otimes\bm{I})+(\bm{I}\otimes\bm{\tau}_{dev,n+1}^{trial})]$
$\bm{\mathbbm{c}}_{vol,n+1}^{e,trial}=K[J^{2}\bm{I}\otimes\bm{I}-(J^{2}-1)\bm{II}]$
Consistent tangent moduli:
$\bm{\mathbbm{c}}_{n+1}^{ep}=\frac{1}{J}[\bm{\mathbbm{c}}_{vol,n+1}^{e,trial}+\bm{\mathbbm{c}}_{dev,n+1}^{e,trial}(1-\Xi_{2})-2\bar{G}\Xi_{3}\bm{n}\otimes\bm{n}-2\bar{G}\Xi_{4}(\bm{n}\otimes\text{dev}[\bm{n}^{2}])^{sym}$

## References

  * [1] Dujc. J and Brank. B.  Stress resultant plasticity for shells revisited.  Computer methods in Applied mechanics and engineering, 247-248:146–165, 2012.
  * [2] Caseiro J. F., Valente R. A. F., Reali A., Kiendl J., Auricchio F, and Alves de Sousa RJ.  Assumed natural strain nurbs-based solid-shell element for the analysis of large deformation elasto-plastic thin-shell structures.  Computer Methods in Applied Mechanics and Engineering, 284:861–880, 2015.
  * [3] Kiendl J., Bletzinger K.-U., J.Linhard, and Wüchner R.  Isogeometric shell analysis with kirchhoff-love elements.  Computer Methods in Applied Mechanics and Engineering, 198:3902–3914, 2009.
  * [4] de Souza Neto EA, Perić D., and Owen D. R. J.  Computational methods for plasticity.  John Wiley & Sons Ltd, 2008.
  * [5] Kiendl J., M.-C. Hsu, M.C.H. Wu, and Reali A.  Isogeometric kirchhoff-love shell formulations for general hyperelastic materials.  Computer Methods in Applied Mechanics and Engineering, 291:280–303, 2015.
  * [6] M. Ambati, Kiendl J., and L. D. Lorenzis.  Isogeometric kirchhoff-love shell formulation for elasto-plasticity.  Comput. Methods Appl. Mech. Engrg, 340:320–339, 2018.
  * [7] E.N. Dvorkin and Klaus-Jürgen Bathe.  A continuum mechanics based four-node shell element for general nonlinear analysis.  Engineering Computations, pages 11367–1378, 1984.
  * [8] Klaus-Jügen Bathe Yeongbin Ki, Phill-Seung Lee.  The mitc4+ shell element and its performance.  Computers and Structures, pages 57–68, 2016.
  * [9] Klaus-Jügen Bathe Yeongbin Ki, Phill-Seung Lee.  The mitc4+ shell element in geometric nonlinear analysis.  Computers and Structures, pages 1–14, 2017.
  * [10] Ted Belytschko Pedro. A. Areias, J.H. Song.  A finite strain quadrilateral shell element based on discrete kirchhoff-love constrains.  International Journal for Numerical Methods in Engineering, pages 1166–1206, 2005.
  * [11] Ted Belytschko Pedro. A. Areias, J.H. Song.  Analysis of fracture in thin shells by overlapping paired elements.  Comput. Methods Appl. Mech. Engrg, pages 5343–5360, 2006.
  * [12] Kiendl J., Bazilevs Y., M.-C. Hsu, Wüchner R., and Bletzinger K.-U.  The bending strip method for isogeometric analysis of kirchhoff-love shell structures comprised of multiple patches.  Computer Methods in Applied Mechanics and Engineering, 199:2403–2416, 2010.
  * [13] Michael J Borden, Michael A Scott, John A Evans, and Thomas JR Hughes.  Isogeometric finite element data structures based on bézier extraction of nurbs.  International Journal for Numerical Methods in Engineering, 2011\.
  * [14] Simo J. C.  A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. part i: continuum formulation.  Comput Methods Appl MechEng, 66:199–219, 1988.
  * [15] Simo J. C.  A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. part ii: computational aspects.  Comput Methods Appl MechEng, 68:1–31, 1988.
  * [16] Simo J. C. and Hughes T. J. R.  Computational Inelasticity.  Springer-Verlag, NewYork, 1998.
  * [17] Dodds Jr. RH.  Numerical techniques for plasticity computations in finite element alnalysis.  Computers and Structures, 1987.
  * [18] Hughes T. J. R., Cottrell J. A., and Bazilevs Y.  Isogeometric analysis: Cad, finite elements, nubrs, exact geometry and mesh refinement.  Computer Methods in Applied Mechanics and Engineering, 194:4135–4195, 2005.
  * [19] Kiendl J.  Isogeometric analysis and shape optimal design of shell strutures.  PhD thesis, Technical University Munich, 2011.
  * [20] Crisfield M. A.  A fast incremental/iterative solution procedure that handles snap through.  Computer and Structures, 13:55–62, 1981.
  * [21] S.H. Lo K.Y. Sze, X.H. Liu.  Popular benchmark problems for geometric nonlinear analysis of shells.  Finite Elem. Anal. Des, page 1551–1569, 2004.
  * [22] F.B. Damjanic B. Brank, D. Peric.  On implementation of a nonlinear four node shell !nite element for thin multilayered elastic shells.  Finite Elem. Anal. Des, page 341–359, 1995.
  * [23] Schwarzea M. and Reese S.  A reduced integration solid-shell fnite element based on the eas and the ans concept - large deformation problems.  International Journal for Numerical Methods in Engineering, 85:289–329, 2011.
  * [24] Brank B., Perić D., and Damjanic F. B.  On large deformation of thin elastoplastic shells: implementation of a finite rotation model for quadrilateral shell elements.  International Journal for Numerical Methods in Engineering, 40:689–726, 1997.
  * [25] P. Areias, T. Rabczuk, J.Cesar de Sa, and J.L. Alves.  Semi-implicit finite strain constitutive integration and mixed strain/stress control based on intermediate configurations.  Engineering Structures, 124:344–360, 2016.
  * [26] Ted Belytschko J.H. Song.  Dynamic fracture of shells subjected to impulsive loads.  Journal of applied mechanics-Transactions of the ASME, 2009.
  * [27] Ted Belytschko J.H. Song, Pedro. A. Areias.  A method for dynamic crack and shear band propagation with phantom nodes.  International Journal for Numerical Methods in Engineering, pages 868–893, 2006.
  * [28] M.A. Msekh P. Areias, T. Rabczuk.  Phase-field analysis of finite-strain plates and shells including element subdivision.  Comput. Methods Appl. Mech. Engrg, pages 322–350, 2016.
  * [29] C. Linder J. Reinoso, M Paggi.  Phase field modeling of brittle fracture for enhanced assumed strain shells at large deformations: formulation and finite element implementation.  Comput Mech, pages 981–1001, 2017.


# Graph-based Neural Modules to Inspect Attention-based Architectures:
A Position Paper

Breno W. Carvalho,1 Artur S. d’Avila Garcez,2 Luís C. Lamb3

###### Abstract

Encoder-decoder architectures are prominent building blocks of state-of-the-
art solutions for tasks across multiple fields where deep learning (DL) or
foundation models play a key role. Although there is a growing community
working on the provision of interpretation for DL models as well as
considerable work in the neuro-symbolic community seeking to integrate
symbolic representations and DL, many open questions remain around the need
for better tools for visualization of the inner workings of DL architectures.
In particular, encoder-decoder models offer an exciting opportunity for
visualization and editing by humans of the knowledge implicitly represented in
model weights. In this work, we explore ways to create an abstraction for
segments of the network as a two-way graph-based representation. Changes to
this graph structure should be reflected directly in the underlying tensor
representations. Such two-way graph representation enables new neuro-symbolic
systems by leveraging the pattern recognition capabilities of the encoder-
decoder along with symbolic reasoning carried out on the graphs. The approach
is expected to produce new ways of interacting with DL models but also to
improve performance as a result of the combination of learning and reasoning
capabilities.

## Introduction

Figure 1: Standard transformer architecture example, adapted from (Vaswani et
al. 2017) with a few possible entry points for introspection modules. Each one
of those prospective entry points carries information about the model state
and context. In 1) the decoder part of the transformer uses the embedded
context information from the encoder to generate the output sequence in the
transduction task. In 2) after the masked multi-head attention module of the
_encoder_ component to inspect which concepts are more strongly related than
others in the construction of contextual information. In 3) information from
the masked multi-head attention modules of the _decoder_ component is used to
inspect which concepts are more relevant to processing the context in light of
the current output state.

We live in hustling times for Artificial Intelligence (AI) research with many
optimistic perspectives from researchers from Computer Science and other
areas. We are entering an era where our models not only can generalize from
examples of a given task, but given the appropriate context and conditions,
they can generalize across different tasks, indicating an emergent phenomenon
that is almost impossible to predict. Such models, notably deep learning (DL)
ones are powerful and influential, yet the current drawbacks in design and
reliability are now obvious. Those models, like fun-house mirrors, can distort
aspects of their training data and generate false affirmations even though
they were fed with the correct facts necessary for the answer. They can also
lean toward unethical biases and make unsupported claims. It is therefore
important to be able to edit a model’s answer and context to first understand
the supporting facts or assumptions behind a given answer and second, correct
or provide new contextual information so as to increase trust in the model’s
outputs.

One approach in this direction is to create introspection entry points into
the models so as to peak into the operations of an intermediate layer (as
depicted in Figure 1), interacting directly with it via graph editing. Those
introspection entry points are neural layers not aimed to contribute to
performance in any single task but to add an extra layer to the developer’s
interaction with the model by providing an abstraction to inspect specific
parts of the network. There are already a few relevant methods in the
literature, some of which are listed in the survey (Lamb et al. 2020) and use
graph representations as part of the network structure itself. At the same
time, the field of neuro-symbolic (NeSy) AI research investigates ways of
merging explicit symbolic knowledge with many such sub-symbolic neural
approaches (d’Avila Garcez and Lamb 2020).

At the Robert S. Engelmore Memorial Lecture, during the AAAI Conference on
Artificial Intelligence, New York, February 10th, 2020, in a talk entitled
_The Third AI Summer_ , Henry Kautz presented a taxonomy for NeSy models, with
six main categories (Kautz 2022). Adding those inspection nodes, which ideally
should not harm model performance at the target task, falls into the _Neuro;
Symbolic_ category in Kautz’s taxonomy. This category includes systems where a
neural component and a symbolic component mutually share information, each
with a different goal. In this way, the inspection module can help edit the
implicit knowledge stored in the DL model and serve even as a debugging tool.
One possible approach to building such inspection modules is based on Graph
Convolutional Neural Networks (GNN) (Liu and Zhou 2020; Wu et al. 2020). By
combining inspection with the model one can perform a secondary task of
internal state-representation of the model (Lamb et al. 2020). As DL models
grow in size and complexity, we believe that those introspection tools in the
network may allow for a more in-depth understanding of such models. This
raises the question of the form of introspection and extensions thereof that
might be expected to produce a better understanding of DL models.

## Short-term Objectives

Our immediate objective is to design and develop an initial tool set for
interpreting encoded knowledge in deep learning (large) language models and
editing it. We aim to add components to infer a graph structure from specific
tensors, so-called entry points, in the model’s architecture. In this way, a
natural step is to evaluate the feasibility of implementing different
approaches for those graph entry points (see Figure 2). We call entry points
compact and non-sparse tensors along the network pipeline that can also be
connected to a GNN to produce a graph abstraction of this stage of the overall
network. For now, we consider only DL architectures based on the overall
encoder-decoder pattern. Those are usually models used in transduction tasks,
which make up a vast class of DL models. Our short-term goals are:

#### To compare different alternatives for implementing graph constraints
without significantly reducing performance.

Creating graph constraints to specific entry points means ensuring that the
underlying tensor always has a meaningful compact graph representation if the
overall model was fed with valid input vectors. To this end, we contemplate
adopting a few alternatives, such as having an auxiliary loss function that
will signal whether the encoder was able to create a valid graph. It will also
signal if the decoder can generate the expected output. Another approach is to
train a GNN decoder and apply it to the output of the encoder part. Aiming at
language tasks supports both kinds of strategies since there are extensive
Language Resources that can be used as an initial graph structure.

#### To find semantic representations for the edges and nodes.

A representation of one embedding without any semantics associated with it
might be as daunting as looking at the model itself. We seek to adapt existing
work in the literature to be able to create a representation with more
meaningful names for entity nodes and connections.

#### To be able to edit the model to increase performance, remove biases or
add new knowledge.

Finally, any editing of the graph representation should be reflected in the
encoder’s weights and influence the final output of the model, ideally in a
way that suggests a causal relationship. This enables us to interact actively
and edit context within the model to ensure that the main component of a given
output is considered to be relevant and correct by an expert in the domain of
interest.

#### To Design a Method for generating interpretations in a local scope.

In the short term, we are concerned with _local interpretations_ where we
interpret the model components’ states and context while processing a specific
input. This approach is already insightful for a series of use cases and might
be the stepping-stone for developing a global interpretation.

Achieving these goals will require the development of a simple but versatile
tool-set that might be useful to the community for inspecting and editing
models in diverse research or industrial areas (similarly to Grad-Cam
(Selvaraju et al. 2017) and Bertology (Rogers, Kovaleva, and Rumshisky 2020)
and many other helpful visualization tools found to be relevant by DL
practitioners).

## Long-term Objectives

Our main goal is to be able to inspect large deep learning language models
through different entry points in such a way as to increase trust in their
output. Transformer or transformer-inspired architectures dominate the state-
of-the-art of such language models. Considering that most transformer
architectures follow an encoder-decoder model, fulfilling the above short-term
objectives should provide an entry point into inspecting the inner state of
the model and its context. We conjecture that having a few editable graph
layers attached to a transformer of considerable size will enable an expert to
study the relationships of the many abstraction levels one might find by
inspecting those layers. Our long-term goals are:

Figure 2: Encoder-decoder example with an entry point for an interpretation
module. One might explore multiple representation languages depending on the
resources available. In this case, the model outputs an incorrect answer for a
question from the TruthfulQA dataset, (Lin, Hilton, and Evans 2021), due to
incomplete context information, that we hope to make explicit with the graph
representation. Figure 3: Illustrative representation of neuro-symbolic
reasoning to complement the pattern learning capabilities of attention-based
models. A GNN entry point is interpreted as an AI system 2 component via the
application of a symbolic reasoner.

#### To have a method to enable experts to visualize, edit and interact
directly with the language models.

One might consider exploring multiple entry points. Manually editing context
or knowledge in the network through a graph interface linked to different
layers might provide much richer insight. Possibilities for graph-entry points
are depicted in Figure 1. The goal is to select tensors that are most
descriptive of the model’s current state. One can also use this method with
visualization methods such as (Rogers, Kovaleva, and Rumshisky 2020).

#### To design a new family of neurosymbolic models based on systematic
interactions with fragments of the underlying architecture.

One might also leverage the power of symbolic reasoners to perform expansion
and inference on the graph context representation to fix or enhance the
model’s output. This approach is much aligned with a fast-slow AI perspective,
where the learning of the language model represents the fast, pattern-
matching, system 1, and the slower, graph-based reasoning represents system 2
(see Figure 3). In this figure, we illustrate that a symbolic expansion of the
graph representation of specific parts of the network might be used to provide
the model with enough information to go from a wrong answer, as depicted in
Figure 2, to the correct answer (see Figure 3).

#### To explore alternative representations of the entry points that are most
useful in a specific context.

As we explore larger language models we might find that different entry points
might benefit from different representation formats. It might even be the case
that reusing smaller transformers will allow us to have natural language
representation for some model entry points. An analogy would be when humans
perform complex tasks and mentally talk to themselves while performing the
task. Another class of representation alternatives to be explored consists of
hyper-graphs that might have a more suitable role in explaining the network’s
inference mechanism since they are able to convey recursive levels of
abstraction in the information being represented.

#### To develop a method where the training data gives the representation
semantics and not necessarily by external language resources.

In an ideal scenario, the only information available during training should be
sufficient to build meaningful representations. Even with available linguistic
resources, expanding the entry point representation with tokens and concepts
from the training data could allow the use of a more robust and extensible
representational language such as a graph, hypergraph, or natural language.

#### To develop a pre-trained zoo of entry point representations for specific
architectures that can be further specialized.

It should be possible to use a library of visualization modules pre-trained on
commonplace text corpora as a baseline, and specialize it with specific domain
corpora. This process might accelerate the inspection process and also
diminish the learning curve of experts from different domains on how to use
those tools. At this point, we are aiming at modularity. We intend not only to
have a methodology but also specific modules that experts from diverse
backgrounds can reuse.

#### To define a method comprising both local and global interpretation scope.

In opposition to local scope interpretations, _global interpretations_ aim at
providing insight into the overall behavior of the model considering an entire
task. For instance, understanding what parts of a large language model are
responsible for summing two numbers or performing sentiment analysis.

To achieve those objectives, one should build a robust methodology to edit and
interact with large language models at multiple levels of abstraction
processed by the model. This methodology might enable further insight into the
inner workings of such models.

## Introspection Modules

Adding introspection modules to current deep learning architectures has the
potential to enable better interpretability, but also unlock new ways to edit
and condition those models. As discussed in detail by this joint initiative
(Bommasani et al. 2021), large-scale language models (also called foundation
models for their re-usability across multiple tasks) pose serious research and
societal challenges. We discuss a few of them in this paper’s _Challenges_
section.

### Fairness and accountability goals

Although somewhat hard to quantify, fairness and accountability are important
to acknowledge in what concerns trust in deep learning models of far-reaching
decisions. Unfairness is characterized in (Du et al. 2020) into two classes:
prediction outcome discrimination, and prediction quality disparity. An
example of the first class is when a model predicts an unfavorable treatment
of a demographic group, such as bias against women. The second class refers to
a model performing poorly for a determined group of individuals. Both cases
raise the need to flag such models and correct them. Accountability might be
understood as a capability for _post-hoc_ inspection of the model in order to
make it available for auditing the causes of the model’s outcome (Chakraborty
et al. 2017).

Following (Bommasani et al. 2021), Section 4.11, which discusses the concept
of interpretability of foundation models, we contemplate three categories of
interpretability: _what_ (the model’s limits and capabilities), _why_ (finding
what in the data set is responsible for the model’s output), and _how_
(gaining an understanding of how certain parts and mechanisms of the model
impact the output). In this work, we focus on the _why_ and _how_
interpretability questions.

As well as being able to interpret those models, one should also seek to
obtain better tools to assist us in accelerating the construction or
adaptation of new models.

### Model editing and design goals

By understanding individual model mechanisms, we seek to build a compositional
understanding of the complex behavior of a foundation model. A possible way is
to interconnect or build larger models from pre-trained components and reduce
societal risks by building more reliable and interpretable models through
“debugging” tools.

## Challenges

Creating those representation/introspection modules raises several challenges.
Some of these, such as the visual interpretation of self-attention heads, have
been addressed to some extent in the literature. Others remain as yet to be
explored.

The _why_ and _how_ interpretability questions from (Bommasani et al. 2021),
as mentioned in the previous section, are essential for the understanding of
large language models, but face many open challenges:

### Challenges on describing _why_ (model behavior) and _how_ (model
mechanisms)

Studying the influence of the model’s input might be performed by carefully
studying how the behavior of the model changes with changes in the input. The
ability to peek inside the model also might be used to provide insight into
which parts of the input are more relevant to the model. The task of
generating interpretations for the inner workings of a language model poses
further challenges:

#### 1) It is not clear how the components of those language models are
interconnected.

Although the overall design of a model’s architecture may be clear, the actual
weights of such units have emergent behavior that might be fairly hard to
predict. A central question to a fair representation of the mechanisms inside
a language model is that of to what extent those mechanisms behave as one
coherent model or as many models. As the model learns to generalize across
multiple tasks, it might be the case that its components find weight sub-
spaces that virtually correspond to many sub-models, possibly one for each
task. It isn’t clear yet to which point in this spectrum current foundation
models belong. This lack of understanding and clarity leads to our next
challenge.

#### 2) Local and global interpretation methods might have pitfalls related to
the complexity of such language models.

Given the highly complex behavior of those models and their emergent
properties, it is unwise to jump to conclusions based on a single global or
local interpretation method, even though there are both local(Fong and Vedaldi
2017; Koh and Liang 2017) and global(Lapuschkin et al. 2019; Thomas et al.
2019) methods for providing insights into task-specific models. For instance,
considering that those large models might find weight sub-spaces for distinct
tasks it might be the case that studying the model for a few of those sub-
tasks doesn’t provide a fair picture of its behavior in other tasks or even
for other similar inputs.

#### 3) The inter-relationship of different components in the model might be
counter-intuitive or hard to determine.

It is not clear whether semantically related mechanisms, e.g. the parts of a
language model responsible for summing up two numbers, are the same as the
ones used to perform a related task, such as e.g. summing two numbers
expressed as numeral nouns.

### Challenges about the introspection modules themselves

Besides the inherent challenges of deriving interpretations from huge language
models, there are also intrinsic points to be evaluated in the design of the
introspection model. There are questions and decisions to be made about the
desired representational language and format, and there are also
considerations about the architecture of those modules and their refinement
procedure.

#### 1) Creating meaningful interpretations assumes an underlying vocabulary
and semantics.

This can be thought of as the symbol grounding problem, i.e. the definition of
a vocabulary for the representation generated by the introspection modules.
This might be achieved through the use of an underlying ontology that will
constrain the interpretation space into a set of predefined concepts and
vocabulary. It remains unclear how the semantics of the ontology will be
transferred to the interpretation. It might be achieved as a secondary task
happening during the training of the language model that we want to inspect.

#### 2) Different abstraction levels must be considered in order to avoid
overly complicated and misguiding interpretations.

This happens because the size of a given interpretation may be prohibitive,
both in terms of meaningful visualization and interpretability. To circumvent
this one needs to employ different levels of abstraction, although those model
layers closer to the output layers seem to pertain to more abstract features
of the data set. Walking back and forth between those different levels of
abstraction reinforces the need to guarantee consistency in the
interpretations’ vocabulary from different entry points.

#### 3) Interpretations should be consistent and reliable across a specific
domain or data set.

This is not an easy requirement to enforce. It might be the case that given
the complexity of those models, consistency can be enforced only partially for
most of the time. We also need to consider consistency across related input or
context. For instance, it should not be the case that two very similar text
prompts result in very different interpretations when looking at the same
entry point in the model. A possible approach is to anchor those graph
representations on concepts and constructs from public linguistic resources,
such as FrameNet (Baker, Fillmore, and Lowe 1998), VerbNet (Palmer, Kipper et
al. 2004), WordNet (Miller 1995) and Propbank (Palmer, Gildea, and Kingsbury
2005), all connected by the SemLink project (Palmer 2009). The process of
generating a consistent interpretation must also be subject to auditing to
ensure the faithfulness of the entire system.

#### 4) Considering that the introspection module can be viewed as a black-box
component, one must devise metrics of trustworthiness specific to it.

The GNN or other DL model used to generate representations from parts of the
language model is also a black-box module that might pose itself some fairness
questions (Jacovi and Goldberg 2020). We hope that it will prove feasible to
keep those models relatively small. It is important to make sure that relying
on those modules will be much less complex than the model being studied. In
such case, existing visualization methods such as (Rogers, Kovaleva, and
Rumshisky 2020; Shi et al. 2021) may suffice to audit them and improve trust.

In summary, the goal of adding a human-editable and interpretable
representation of context and the state of specific parts of a language model
has many open questions that to be solved will require drawing contributions
from different areas of expertise. One way to face some of the challenges
outlined here is to combine approaches that use both linguistic resources to
provide the needed representation building blocks (e.g. concepts and
relations) and to use tokens and concepts from the training data set as
building blocks.

## Outlook

As mentioned in previous sections, this kind of neural-symbolic work draws
from a diverse mixture of research areas in AI and Computer Science more
generally.

### Peeking inside DL models:

While characterizing and describing the behavior of a deep learning model, one
might consider the data set and the predicted outputs, or also the state
variations of the models. It is also possible to build interpretations about
the model behavior on specific data samples or trying to understand how it
behaves in a broader sense when performing a particular task.

Local scope interpretation methods target insight into the model’s behavior on
specific inputs or contrasting a given set of data points. Examples include
(Koh and Liang 2017), an approach that uses influence functions from robust
statistics to understand how perturbations in the training set influence the
model behavior; (Fong and Vedaldi 2017) uses the concept of interpretations as
meta-predictors, i.e. one can use an explanation to predict the output of a
model, and they instantiate this approach using image classification masks. By
contrast, global methods are used to describe general aspects and behaviors of
the model. Interesting examples include (Lapuschkin et al. 2019) and (Thomas
et al. 2019). Both local and global methods treat the language models as
black-box objects to be understood through their inputs and outputs.

Other than looking at the training set and the models’ outputs, another class
of interpretation methods that is more aligned with the work proposed here,
consists of models that extract information from the model weights. This
includes GradCam (Selvaraju et al. 2017), which does backward propagation of
gradients to estimate which pixels of the input where more relevant for a
Convolutional Network output; (Rogers, Kovaleva, and Rumshisky 2020) provides
a tool-set for visualizing the weights of the self-attention heads of BERT-
like models. Although those methods can be very insightful when analyzing
task-specific models, considering that current language models might have
distinct behavior across different tasks, it is unclear how insightful those
methods could be in this scenario. Most likely, the designers and analysts of
such models would need to employ multiple methods of interpretation and
visualization to have a fair understanding of the models’ behavior.

As an alternative to the line of work described in this paper, where we
propose the use of external (and hopefully auditable) modules to peak into the
language models, one can train the model itself to generate an explanation for
its output (Elton 2020). This approach faces some skepticism since language
models can generate plausible, although untrue, outputs and nothing prevents
this behavior from extending to the model’s self-explanations.

### Graph deconvolution and Neurosymbolic AI

We regard GNNs and neuro-graph algorithms in general as a promising general
approach for both local and global interpretation.

Graph convolution is defined by analogy to convolutional layers over Euclidean
data (Scarselli et al. 2008; Sanchez-Lengeling et al. 2021) and it is an
important tool to infuse deep learning models with relational knowledge (Lamb
et al. 2020). Conversely, we conjecture that using specific decoder units to
extract graph representations from the entry point tensors (even if having to
condition the training of the underlying language model) might be a worthwhile
approach to obtain meaningful interpretations from language models.

We might consider each introspection module as a simple decoder unit that maps
the encoded embedding into the interpretation space. As mentioned in the
Challenges section, one might initially use linguist resources conditioned on
the sentence processed by the language model to assemble the vocabulary used
to form the interpretations. A possible approach is to use heuristics to get
concepts from the language model prompt and build a simple graph, even
combining it with a syntax or dependency tree to have an initial graph to use
and evaluate the decoder unit used as an introspection module. In this sense,
there are a few approaches in the literature that one might draw inspiration
from. For instance, there is work on graph encoders (Hamilton, Ying, and
Leskovec 2017), and also deconvolutional graph networks (Li et al. 2021) to
derive a graph from the entry point embedding. At this stage, it might be
unclear how to define the appropriate abstraction level of the interpretation.

Instead of building the modules as decoder units trained on a heuristic graph,
another approach is to adapt variational graph autoencoders (Kipf and Welling
2016). A variation of the encoding part could output a Graph Convolutional
Network-like embedding, and the decoder would try to retrieve the original
input of the language model (assuming that the entry point embedding might
encode a smooth topological representation of concepts). This assumption is
not guaranteed to hold, but it is an interesting issue to investigate. One
might use deconvolution networks (Li et al. 2021) to infer a graph without
prior conditioning. The intuition behind this rationale is that this
deconvolution approach would create representations without any supervision,
which might be too complicated or obscure to be used as interpretation.

As stated in our short and long-term objectives, inferring an arbitrary graph
from those entry points is not enough for our interpretation and editing
purposes. Those graphs also need to have a clear semantics that correlates
with the model’s output and with commonsense and domain knowledge as well.

One might wonder why bother with graph representations and not try to go
directly to building natural language interpretations. We conjecture that
natural language will lack explicit restrictions that are needed for clarity
of a given interpretation and also that the resulting natural language
introspection module would need to be almost as complex as the language model
itself, thus defeating its original purpose.

There are several neurosymbolic methods that aim to infuse the deep learning
model architectures with symbolic or logical structures (d’Avila Garcez and
Lamb 2020), including logic tensor networks (LTN) (Serafini and Garcez 2016)
and logical neural networks (LNN) (Riegel et al. 2020). They can also be used
as parts of the introspection modules described here in order to make DL
models more auditable and reliable.

## Conclusion

Large language models are becoming state-of-the-art in many tasks and fields,
and as they grow in popularity, they also grow in complexity. Those models
usually follow a transduction encoder-decoder pattern based on self-attention
mechanisms that are repeated multiple times to a point where they are so
complex that surprising behavior emerges from them. In this paper, we state
the need for inspection methods to both characterize and describe those
models’ behaviors. Despite the consolidated literature on deep learning
interpretability, the community still needs to build different tools to
inspect different aspects of those models. We believe that graph-based
algorithms, in particular, deconvolution graph networks and variational graph
auto-encoders might be key to generating formal yet flexible interpretations
of parts of such large language models.

Interpretability of language models is a blooming field as those models keep
growing and reaching multiple uses in society. There is a vast literature to
support this field, but there is a flagrant need for new methods and
approaches to deal with such models’ high generalizability. Being able to
audit and edit aspects of the model during development and deployment should
allow the community to correct flawed inferences performed by those models and
adjust to counter certain unethical biases. Thus inspection tools have the
potential for tremendous impact if they can help us develop models we can
trust and rely on to shape the further development of AI.

We envision a roadmap for building tools with which multidisciplinary
communities can start inspecting small transduction models based on self-
attention, and progressively scale up to large-scale language models. By doing
so, those communities are empowered to shape the widely spread use cases of
such models in society. To this end, we hope to start a conversation to unify
our efforts, lower the threshold for other machine learning researchers to
join us, and bring these communities closer together with a common language.

## Acknowledgments

We thank Viviane Torres, Sandro Rama Fiorini, Emílio Ashton Brasil, and Renato
Cerqueira for insightful conversations on the topic. Luis Lamb was supported
in part by CAPES and CNPq, Brazil.

## References

  * Baker, Fillmore, and Lowe (1998) Baker, C. F.; Fillmore, C. J.; and Lowe, J. B. 1998.  The berkeley framenet project.  In _COLING 1998 Volume 1: The 17th International Conference on Computational Linguistics_.
  * Bommasani et al. (2021) Bommasani, R.; Hudson, D. A.; Adeli, E.; Altman, R.; Arora, S.; von Arx, S.; Bernstein, M. S.; Bohg, J.; Bosselut, A.; Brunskill, E.; Brynjolfsson, E.; Buch, S.; Card, D.; Castellon, R.; Chatterji, N. S.; Chen, A. S.; Creel, K.; Davis, J. Q.; Demszky, D.; Donahue, C.; Doumbouya, M.; Durmus, E.; Ermon, S.; Etchemendy, J.; Ethayarajh, K.; Fei-Fei, L.; Finn, C.; Gale, T.; Gillespie, L.; Goel, K.; Goodman, N. D.; Grossman, S.; Guha, N.; Hashimoto, T.; Henderson, P.; Hewitt, J.; Ho, D. E.; Hong, J.; Hsu, K.; Huang, J.; Icard, T.; Jain, S.; Jurafsky, D.; Kalluri, P.; Karamcheti, S.; Keeling, G.; Khani, F.; Khattab, O.; Koh, P. W.; Krass, M. S.; Krishna, R.; Kuditipudi, R.; and et al. 2021.  On the Opportunities and Risks of Foundation Models.  _CoRR_ , abs/2108.07258.
  * Chakraborty et al. (2017) Chakraborty, S.; Tomsett, R.; Raghavendra, R.; Harborne, D.; Alzantot, M.; Cerutti, F.; Srivastava, M.; Preece, A.; Julier, S.; Rao, R. M.; et al. 2017.  Interpretability of deep learning models: A survey of results.  In _2017 IEEE smartworld, ubiquitous intelligence & computing, advanced & trusted computed, scalable computing & communications, cloud & big data computing, Internet of people and smart city innovation (smartworld/SCALCOM/UIC/ATC/CBDcom/IOP/SCI)_, 1–6. IEEE.
  * d’Avila Garcez and Lamb (2020) d’Avila Garcez, A.; and Lamb, L. C. 2020.  Neurosymbolic AI: The 3rd Wave.  _CoRR_ , abs/2012.05876.
  * Du et al. (2020) Du, M.; Yang, F.; Zou, N.; and Hu, X. 2020.  Fairness in deep learning: A computational perspective.  _IEEE Intelligent Systems_ , 36(4): 25–34.
  * Elton (2020) Elton, D. C. 2020.  Self-explaining AI as an alternative to interpretable AI.  In _International conference on artificial general intelligence_ , 95–106. Springer.
  * Fong and Vedaldi (2017) Fong, R. C.; and Vedaldi, A. 2017.  Interpretable explanations of black boxes by meaningful perturbation.  In _Proceedings of the IEEE international conference on computer vision_ , 3429–3437.
  * Hamilton, Ying, and Leskovec (2017) Hamilton, W. L.; Ying, R.; and Leskovec, J. 2017.  Representation learning on graphs: Methods and applications.  _arXiv preprint arXiv:1709.05584_.
  * Jacovi and Goldberg (2020) Jacovi, A.; and Goldberg, Y. 2020.  Towards faithfully interpretable NLP systems: How should we define and evaluate faithfulness?  _arXiv preprint arXiv:2004.03685_.
  * Kautz (2022) Kautz, H. A. 2022.  The Third AI Summer: AAAI Robert S. Engelmore Memorial Lecture.  _AI Mag._ , 43(1): 93–104.
  * Kipf and Welling (2016) Kipf, T. N.; and Welling, M. 2016.  Variational graph auto-encoders.  _arXiv preprint arXiv:1611.07308_.
  * Koh and Liang (2017) Koh, P. W.; and Liang, P. 2017.  Understanding black-box predictions via influence functions.  In _International conference on machine learning_ , 1885–1894. PMLR.
  * Lamb et al. (2020) Lamb, L. C.; d’Avila Garcez, A. S.; Gori, M.; Prates, M. O. R.; Avelar, P. H. C.; and Vardi, M. Y. 2020.  Graph Neural Networks Meet Neural-Symbolic Computing: A Survey and Perspective.  In _IJCAI 2020_ , 4877–4884. ijcai.org.
  * Lapuschkin et al. (2019) Lapuschkin, S.; Wäldchen, S.; Binder, A.; Montavon, G.; Samek, W.; and Müller, K.-R. 2019.  Unmasking Clever Hans predictors and assessing what machines really learn.  _Nature communications_ , 10(1): 1–8.
  * Li et al. (2021) Li, J.; Li, J.; Liu, Y.; Yu, J.; Li, Y.; and Cheng, H. 2021.  Deconvolutional networks on graph data.  _Advances in Neural Information Processing Systems_ , 34: 21019–21030.
  * Lin, Hilton, and Evans (2021) Lin, S.; Hilton, J.; and Evans, O. 2021.  TruthfulQA: Measuring How Models Mimic Human Falsehoods.  _CoRR_ , abs/2109.07958.
  * Liu and Zhou (2020) Liu, Z.; and Zhou, J. 2020.  Introduction to graph neural networks.  _Synthesis Lectures on Artificial Intelligence and Machine Learning_ , 14(2): 1–127.
  * Miller (1995) Miller, G. A. 1995.  WordNet: A Lexical Database for English.  _Commun. ACM_ , 38(11): 39–41.
  * Palmer (2009) Palmer, M. 2009.  Semlink: Linking propbank, verbnet and framenet.  In _Proceedings of the generative lexicon conference_ , 9–15. GenLex-09, Pisa, Italy.
  * Palmer, Gildea, and Kingsbury (2005) Palmer, M.; Gildea, D.; and Kingsbury, P. 2005.  The proposition bank: An annotated corpus of semantic roles.  _Computational linguistics_ , 31(1): 71–106.
  * Palmer, Kipper et al. (2004) Palmer, M.; Kipper, K. L.; et al. 2004.  Verbnet.  _The Oxford Handbook of Cognitive Science_.
  * Riegel et al. (2020) Riegel, R.; Gray, A.; Luus, F.; Khan, N.; Makondo, N.; Akhalwaya, I. Y.; Qian, H.; Fagin, R.; Barahona, F.; Sharma, U.; et al. 2020.  Logical neural networks.  _arXiv preprint arXiv:2006.13155_.
  * Rogers, Kovaleva, and Rumshisky (2020) Rogers, A.; Kovaleva, O.; and Rumshisky, A. 2020.  A primer in bertology: What we know about how bert works.  _Transactions of the Association for Computational Linguistics_ , 8: 842–866.
  * Sanchez-Lengeling et al. (2021) Sanchez-Lengeling, B.; Reif, E.; Pearce, A.; and Wiltschko, A. B. 2021.  A gentle introduction to graph neural networks.  _Distill_ , 6(9): e33.
  * Scarselli et al. (2008) Scarselli, F.; Gori, M.; Tsoi, A. C.; Hagenbuchner, M.; and Monfardini, G. 2008\.  The graph neural network model.  _IEEE transactions on neural networks_ , 20(1): 61–80.
  * Selvaraju et al. (2017) Selvaraju, R. R.; Cogswell, M.; Das, A.; Vedantam, R.; Parikh, D.; and Batra, D. 2017.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In _Proceedings of the IEEE international conference on computer vision_ , 618–626.
  * Serafini and Garcez (2016) Serafini, L.; and Garcez, A. d. 2016.  Logic tensor networks: Deep learning and logical reasoning from data and knowledge.  _arXiv preprint arXiv:1606.04422_.
  * Shi et al. (2021) Shi, X.; Lv, F.; Seng, D.; Zhang, J.; Chen, J.; and Xing, B. 2021.  Visualizing and understanding graph convolutional network.  _Multimedia Tools and Applications_ , 80(6): 8355–8375.
  * Thomas et al. (2019) Thomas, A. W.; Heekeren, H. R.; Müller, K.-R.; and Samek, W. 2019.  Analyzing neuroimaging data through recurrent deep learning models.  _Frontiers in neuroscience_ , 13: 1321.
  * Vaswani et al. (2017) Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, L. u.; and Polosukhin, I. 2017.  Attention is All you Need.  In Guyon, I.; Luxburg, U. V.; Bengio, S.; Wallach, H.; Fergus, R.; Vishwanathan, S.; and Garnett, R., eds., _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc.
  * Wu et al. (2020) Wu, Z.; Pan, S.; Chen, F.; Long, G.; Zhang, C.; and Philip, S. Y. 2020.  A comprehensive survey on graph neural networks.  _IEEE transactions on neural networks and learning systems_ , 32(1): 4–24.