A Unified Finite Strain Theory for Membranes and Ropes
Thomas-Peter Fries Daniel Schöllhammer







A Unified Finite Strain Theory for 
Membranes and Ropes
    T.P. Fries, D. Schöllhammer
    December 30, 2023
=======================================================





Institute of Structural Analysis

Graz University of Technology

Lessingstr. 25/II, 8010 Graz, Austria







The finite strain theory is reformulated in the frame of the Tangential
Differential Calculus (TDC) resulting in a unification in a threefold
sense. Firstly, ropes, membranes and three-dimensional continua are
treated with one set of governing equations. Secondly, the
reformulated boundary value problem applies to parametrized and
implicit geometries. Therefore, the formulation is more general than
classical ones as it does not rely on parametrizations implying curvilinear
coordinate systems and the concept of co- and contravariant base vectors.
This leads to the third unification: TDC-based models are applicable
to two fundamentally different numerical approaches. On the one hand,
one may use the classical Surface FEM where the geometry is discretized
by curved one-dimensional elements for ropes and two-dimensional surface
elements for membranes. On the other hand, it also applies to recent
Trace FEM approaches where the geometry is immersed in a higher-dimensional
background mesh. Then, the shape functions of the background mesh
are evaluated on the trace of the immersed geometry and used for the
approximation. As such, the Trace FEM is a fictitious domain method
for partial differential equations on manifolds. The numerical results
show that the proposed finite strain theory yields higher-order convergence
rates independent of the numerical methodology, the dimension of the
manifold, and the geometric representation type.

Keywords: finite strain theory, ropes, membranes, Trace FEM, fictitious
domain method, embedded domain method, PDEs on manifolds






§ INTRODUCTION


In structural analysis, there are many examples where a dimensional
reduction in the modeling is a key step to achieve simpler models
than considering the full (three-dimensional) continuum. The situation
is simple for flat domains such as straight beams and flat plates.
For curved structures with reduced dimensionality such as curved
beams and shells, modeling becomes considerably more involved and
the governing equations are partial differential equations (PDEs)
on manifolds. The situation is even more complicated for structures
undergoing large displacements such as ropes (cables) and membranes.
These structures do not feature any bending resistance and often deform
largely such that one has to carefully distinguish between the undeformed
and deformed configurations. When equilibrium is to be fulfilled in
the deformed configuration, this is generally referred to as finite
strain theory, the geometrically non-linear situation, or large displacement
analysis.

Herein, we propose a new formulation of the mechanical models for
ropes and membranes in finite strain theory which employs the Tangential
Differential Calculus (TDC) for the definition of geometric and differential
quantities. The geometries of ropes and membranes may be seen as (curved)
manifolds embedded in a higher-dimensional space. Let the dimension
of the surrounding region be d and the dimension of the manifold
be q with q=1 for ropes and q=2 for membranes. The classical
modeling approach to ropes and membranes (including shells) may be
called parameter-based <cit.>, because it relies on
a parametrization which maps from the q-dimensional to the d-dimensional
space. In this case, a curvilinear coordinate system is naturally
implied and co- and contra-variant base vectors are easily defined.
An alternative is to define the geometry implicitly, for example,
based on the level-set method <cit.>.
Then, the zero-level set of a scalar function implies the manifold.
In this case, no parametrization, hence, no curvilinear coordinates
exist on the manifold. This situation cannot be covered by the classical
models but the proposed TDC-based approach can.

When comparing the new TDC-based formulation with the classical parameter-based
formulation, we notice the following important aspects:


  * The first argument is geometrical and is concerned with the
definition of the boundary value problem (BVP) of ropes and membranes.
It is desirable to have a mechanical model which applies to general
geometries no matter whether they are defined in parametric or implicit
form. The classical parameter-based approach does not apply to the
case where the geometry is, for example, defined by a zero-level set.
In that case, a parametrization only becomes available upon meshing
the zero-level set (resulting in an atlas of mappings implied by the
elements). However, then the BVP is rather defined with respect to
a discrete setting than a continuous one. Consequently, one advantage
of the TDC-based formulation is that it allows for a proper definition
of a continuous BVP which applies for both, parametrically and implicitly
defined geometries.

  * The second argument is mechanical. The TDC-based formulation
treats ropes and membranes in a unified sense where all mechanical
quantities such as stress and strain tensors are based on differential
operators formulated with the TDC. These quantities are defined based
on the d-dimensional global coordinate system into which the manifold
is immersed. That is, these tensors are always d× d no matter
whether ropes or membranes are considered. The TDC-based formulation
may also be seen as a special case of a d-dimensional continuum
(non-manifold case) where the manifold-operators become classical
ones. In contrast, in the parameter-based formulation, where the q-dimensional
manifold results from mapping q parameters to d dimensions,
the related mechanical tensors are q× q and thus inherently
different for ropes, membranes and continua <cit.>.
For the TDC-based situation, these tensors result from the same formulas,
however, based on different differential operators depending on the
mechanical situation. 

  * The third argument is numerical and related to discretization
methods. The new formulation allows for two fundamentally different
numerical approaches. One may be seen as the classical approach which
relies on meshing the ropes and membranes by (curved) line or surface
meshes; this is called the Surface FEM herein. It is inherently linked
to the parameter-based formulation of the mechanical models. Due to
the fact that for implicit geometries, one may also provide surface
meshes, and then continue with the parameter-based formulation, there
was not really a strong need for a more general formulation of the
models until recently. However, there is an alternative approach to
solve ropes and membranes which uses d-dimensional, non-matching
background meshes to approximate the displacements instead of conforming,
curved q-dimensional surface meshes. This approach has been labeled
Trace FEM because for the integration of the weak form, the d-dimensional
shape functions are only evaluated on the trace of the manifold <cit.>.
The Trace FEM is inherently linked to implicitly defined manifolds
and does not need any parametrization. The solution of BVPs based
on non-matching background methods is generally referred to as fictitious
(or embedded) domain methods (FDMs) with a large number of variants
existing today. Recently, the Cut FEM has emerged as a popular FDM
allowing for higher-order accuracy <cit.>.
When using the Cut FEM for the solution of PDEs on manifolds as done
herein, the method becomes analogous to the Trace FEM.

Consequently, the TDC-based formulation allows for a unification in
a geometrical (parametric and implicit geometries), mechanical (ropes,
membranes, and d-dimensional continua) and numerical (Surface and
Trace FEM) sense. For similar reasons, the authors have already used
the TDC-based approach to reformulate the mechanical models of linear
Kirchhoff-Love <cit.>
and Reissner-Mindlin shells <cit.>. Using the
TDC for the definition of BVPs on manifolds as discussed herein is
already well-accepted in the definition of transport phenomena on
manifolds where it has replaced the classical parameter-based definition
in many cases <cit.>. It is
thus coming timely and naturally to reformulate mechanical
BVPs on manifolds based on the TDC in the same sense.

The classical theory of large displacement membranes based
on curvilinear coordinates is described, e.g., in <cit.>.
The general equivalence of models using parametrized or implicit manifolds
is outlined in <cit.>. Geometrically exact shell models
based on explicitly defined surfaces and locally using Cartesian coordinates
are defined in the linear setting in <cit.> and for the
non-linear case in <cit.>. Local Cartesian coordinates
are also used in <cit.>, where the initial configuration
is modeled with a stress-free deformation of a flat surface. Another
approach for explicitly defined geometries is the degenerated solid
approach <cit.>. Concerning a TDC-based modeling of large
deformation membranes with the Surface FEM, we emphasize the work
in <cit.>. A TDC-related approach in the field of composite
structures is presented in <cit.> for embedded membranes
using complex material models and with focus on analytical solutions.
In contrast, herein, the novelty is in the unified treatment of ropes
and membranes and the numerical treatment with Surface and
Trace FEM. In this work, we also discuss in detail how ropes and membranes
(in two and three dimensions) are generally defined using the parametric
and the implicit approach, including all relevant geometric and differential
quantities. For example, implicitly defining a rope in three dimensions
requires two level-set functions whereas one level-set function
is sufficient for membranes. Furthermore, the mechanical discussion
presented herein includes stress and strain tensors, so that the manifold
versions of the first and second Piola-Kirchhoff stress tensors, the
Cauchy stress tensor, the Euler-Almansi and Green-Lagrange strain
tensors are explicitly given.

Concerning the numerical approximation of ropes and membranes based
on the Surface and Trace FEM as discussed herein, it is mentioned
that the Trace FEM just as any other FDM simplifies the meshing of
geometries significantly. However, additional effort is required for
(i) the integration of the weak form, wherefore suitable integration
points on the manifolds have to be provided <cit.>,
(ii) the application of boundary conditions which is more involved
because the boundary is within the background elements and instead
of prescribing nodal values one may have to use Nitsche's method <cit.>
or other approaches for enforcing constraints, and (iii) stabilization
terms which are necessary to address the existence of shape functions
with small supports on the manifolds and to find unique solutions
with the background mesh although the BVP is only defined on the manifold
<cit.>. In
spite of the increased effort for implementing the Trace FEM, it may
have significant advantages over the Surface FEM. For example, when
ropes and/or membranes are reinforcing sub-structures embedded in
some three-dimensional continuum, there is no need to consider these
sub-structures in the meshing of the volume. This is the first work
where Trace FEM results for ropes and membranes are shown, enabling
the implicit analysis of these structures.

The paper is organized as follows: In Section <ref>, the
parametric and implicit geometry definitions of manifolds are discussed
in detail for the situation of ropes and membranes. As usual in finite
strain theory, the undeformed and deformed situations are distinguished,
related by the sought displacements. In each configuration, normal
vectors, surface stretches, and differential surface operators are
defined, wherefore it is important to distinguish parametric from
implicit definitions. The mechanical modeling is outlined in Section
<ref> following the classical steps defining stress and
strain tensors and imposing equilibrium in the deformed configuration.
The strong and weak forms of the governing equations are given. It
is noteworthy that ropes and membranes are treated in the same manner
and only the geometry-dependent surface operators defined in Section
<ref> differ. The numerical solution of the governing
equations is considered in Section <ref> where the
discrete weak forms of the Surface and Trace FEM are given. The numerical
results in Section <ref> confirm that both numerical
approaches achieve higher-order convergence rates. It is also confirmed
how ropes and membranes may easily be coupled with the presented formulation.
The paper ends in Section <ref> with a summary and
conclusions.




§ TANGENTIAL DIFFERENTIAL CALCULUS IN FINITE STRAIN THEORY


The focus of this work is on ropes (also called cables) and membranes
undergoing large displacements which is covered by finite strain theory.
Only solids according to the Saint Venant–Kirchhoff material model
are considered herein which may be seen as the simplest extension
of a linear elastic material model to large displacements. Cables
may be modeled as one-dimensional lines in the two- or three-dimensional
space. Membranes are two-dimensional surfaces in the three-dimensional
space. Hence, we may state that cables and membranes are curved
manifolds with a lower dimension q than the surrounding space with
dimension d. In this section, we address the issue of how to define
such manifolds geometrically and formulate (differential) operators
needed later in the mechanical model, see Section <ref>.
For the geometry definition, we separately discuss the situation for
parametrized and implicit manifolds. Implicit manifolds are implied
by level-set(s) of scalar function(s) following the concept of the
level-set method <cit.>. For
the parametric situation, the outline is related to the classical
setup using tensor notation rather than index notation and avoiding
any explicit reference to classical terms in the context of curvilinear
coordinate systems such as co- and contra-variant base vectors and
Christoffel symbols. For the implicit situation, the presented
outline, systematically including all geometric and differential quantitites
for ropes and membranes in finite strain theory, is original.




 §.§ Deformed and undeformed configurations


As usual in finite strain theory, we consider an undeformed material
configuration and a deformed spatial configuration. These are represented
by the q-dimensional manifolds Γ_ X and Γ_ x,
respectively, which are immersed in a d-dimensional space ℝ^d,
herein, d={ 2,3}. The difference d-q is also called
codimension of the manifold. We follow the usual notation to relate
uppercase letters in variable and operator names with the undeformed
configuration and lowercase letters with the deformed one. The displacement
field u( X) relates the two configurations
via 

    x= X+ u( X)  with  X∈Γ_ X⊂ℝ^d and  x∈Γ_ x⊂ℝ^d.






 §.§ Parametrized manifolds


For parametrized manifolds, there exists a map X( r):ℝ^q→ℝ^d
from some lower-dimensional reference domain Ω_ r⊂ℝ^q
to the undeformed configuration Γ_ X⊂ℝ^d.
An important consequence is that local curvilinear coordinate systems
result naturally on the manifolds. It is useful to describe the situation
separately for cables (one-dimensional manifolds) and membranes (two-dimensional
manifolds).







  §.§.§ Two-dimensional manifolds


We start with two-dimensional manifolds (q=2) in the three-dimensional
space (d=3) which is relevant for membranes. Let there be a reference
domain Ω_ r⊂ℝ^2 and a map X( r):Ω_ r→Γ_ X⊂ℝ^3,
see Fig. <ref>. We label the components r=[r,s]^T
and X=[X,Y,Z]^T. The Jacobi matrix 

    J( r)=∂ X( r)/∂ r=∇_ r X( r)=[[ ∂_rX ∂_sX; ∂_rY ∂_sY; ∂_rZ ∂_sZ ]]

has dimensions (3×2). One may easily obtain two
vectors T_1^⋆=∂ X/∂ r and
T_2^⋆=∂ X/∂ s from the columns
of J being tangential to Γ_ X at a mapped point
X. The Jacobi matrix is also used to compute the first fundamental
form G and the operator Q, later needed for defining
surface gradients,

    G= J^T· J     (q× q)-matrix,        Q= J· G^-1 (d× q)-matrix.



Next, we consider a displacement field u( r)
assuming that a point r is given which may also be seen as
a function u( r( X)) when a
point X∈Γ_ X is given (and back-projected to
the reference domain inverting the map X( r)).
We emphasize that in both cases, the displacement field only lives
on the manifold Γ_ X. Hence, no classical partial derivatives
with respect to X may be computed (unless u is smoothly
extended to some neighborhood of X which is not unique and
not considered here) so that the only useful gradient of u
in this context is the surface gradient. For some scalar function
f( r), e.g., each displacement component, the
surface gradient is

    ∇_ X^Γf( r)= Q·∇_ rf( r)  ⇔  [[ ∂_X^Γf; ∂_Y^Γf; ∂_Z^Γf ]]= Q·[[ ∂_rf; ∂_sf ]],

For a vector function u( r)=[u,v,w]^T∈ℝ^3,
we have the directional surface gradient

    ∇_ X^Γ,dir u( r)=[[ (∇_ X^Γu)^T; (∇_ X^Γv)^T; (∇_ X^Γw)^T ]]=[[ ∂_X^Γu ∂_Y^Γu ∂_Z^Γu; ∂_X^Γv ∂_Y^Γv ∂_Z^Γv; ∂_X^Γw ∂_Y^Γw ∂_Z^Γw ]]=∇_ r u( r)· Q^T,

which is to be distinguished from the covariant surface gradient
of a vector field defined in Section <ref>.

Let us next consider the map from the undeformed to the deformed configuration
x( X)= X+ u( X) which
is Γ_ X→ℝ^d. The related Jacobi
matrix is also called surface deformation gradient, 

    F_Γ=∇_ X^Γ,dir x( X)= I+∇_ X^Γ,dir u( X),

where I is a (d× d) identity matrix.

One may obtain all equivalent quantities in the deformed configuration:
The Jacobi-matrix from the reference to the deformed configuration
j= F_Γ· J and the tangent vectors t_1^⋆,
t_2^⋆ to the deformed configuration Γ_ x
at a mapped point x based on the Jacobi matrix j.
Furthermore, the first fundamental form g= j^T· j
and the operator q= j· g^-1 relating the classical
gradient of the reference configuration with the surface gradient
of the deformed configuration as ∇_ x^Γf= q·∇_ rf.

Based on the pairs of tangent vectors in the undeformed and deformed
configuration, one may compute unique normal vectors in each configuration,
N^⋆= T_1^⋆× T_2^⋆ and n^⋆= t_1^⋆× t_2^⋆.
Then, the projectors P( X) and p( x)
are computed as

    P     =     I- N⊗ N  with  N= N^⋆/‖ N^⋆‖,
     p     =     I- n⊗ n  with  n= n^⋆/‖ n^⋆‖.

The same result is obtained when computing a tangent vector T_3^⋆
in the tangent plane spanned by T_1^⋆ and T_2^⋆
being orthogonal to T_1^⋆ using Gram Schmidt orthogonalization,
T_3^⋆= T_2^⋆- T_1^⋆· T_2^⋆/ T_1^⋆· T_1^⋆· T_1^⋆,
then 

    P= T_1⊗ T_1+ T_3⊗ T_3  with  T_1= T_1^⋆/‖ T_1^⋆‖, T_3= T_3^⋆/‖ T_3^⋆‖,

analogously for p. The projector P at some point X
maps an arbitrary vector in ℝ^d to the tangent space
at Γ_ X, hence, P· N=0. P
is symmetric, P= P^T, and idempotent, P· P= P,
which holds analogously for p.

Next, we are interested in the stretch of a differential element of
the membrane when undergoing the deformation. This is interpreted
as an area stretch and defined as

    Λ=√( g)/√( G)=‖ n^⋆‖/‖ N^⋆‖=‖ t_1^⋆× t_2^⋆‖/‖ T_1^⋆× T_2^⋆‖.



Finally, an operator W is introduced which relates surface
gradients of the undeformed and the deformed configuration as

    ∇_ x^Γf= W·∇_ X^Γf  with  W= q( Q^T· Q)^-1 Q^T.

This result is obtained using ∇_ X^Γf= Q·∇_ rf
and ∇_ x^Γf= q·∇_ rf. Note
that Q and q are (d× q)-matrices,
hence, not quadratic and the concept of generalized inverses (or pseudo
inverses) is needed to obtain Eq. (<ref>).







  §.§.§ One-dimensional manifolds


Consider some one-dimensional reference domain Ω_r⊂ℝ
and a map X(r):Ω_r→Γ_ X⊂ℝ^d,
d={ 2,3}, see Fig. <ref>.
This situation applies to cables in two and three dimensions. The
Jacobi matrix J(r)=∇_r X(r)
consists of one column only which implies a tangent vector T∈ℝ^d
being tangential to the undeformed configuration at X. For
the tangent vector in the deformed configuration follows t^⋆= F_Γ· T^⋆.
Most parts of the discussion from Section <ref>
apply accordingly. However, the definition of the projectors changes
to 

    P     =     T⊗ T  with  T= T^⋆/‖ T^⋆‖,
     p     =     t⊗ t  with  t= t^⋆/‖ t^⋆‖.

For the stretch of a differential element of the cable, which may
be seen as a line stretch, there follows

    Λ=√( g)/√( G)=‖ t^⋆‖/‖ T^⋆‖.






 §.§ Implicit manifolds


Implicit manifolds are implied by one or more level-set functions.
Generally speaking, the codimension determines the number of level-set
functions required to define a unique geometry of a manifold. For
the cases relevant in this work, this means that one level-set
function ϕ( X) is required for cables in ℝ^2
and membranes in ℝ^3 which have codimension 1. For
cables in ℝ^3, on the other hand, two level-set
functions ϕ_1( X) and ϕ_2( X)
are needed. We split the discussion depending on the codimension.







  §.§.§ Manifolds with codimension 1


Oriented manifolds with codimension 1 may be defined by one
level-set function ϕ( X). Usually, the zero-level
set of ϕ implies the manifold of interest and there are infinitely
many possible ϕ implying the same geometry. The signed distance
function is a particularly useful concrete example for ϕ( X)
and often used in practice. It is noteworthy that many level-sets
are unbounded in ℝ^d which is not desirable for
the definition of mechanical applications. Fortunately, it is easily
possible to define bounded manifolds by additional (slave) level-set
functions ψ_i( X), see <cit.>.
Then, the undeformed configuration is, e.g., given by 

    Γ_ X={ X∈ℝ^d:ϕ( X)=0, ψ_i( X)≥0, i=1,… n} ,

as shown in Fig. <ref>(a). Even simpler is
to associate a domain of definition Ω_ X⊂ℝ^d
with ϕ, then the bounded manifold results as

    Γ_ X={ X∈Ω_ X:ϕ( X)=0} .

In this case, the boundaries are the intersections of the zero-level
set with the boundary of the domain of definition, see Fig. <ref>(b).
The implicit definition according to Eq. (<ref>)
is used in the following unless noted otherwise.




Next, we focus on the situation in large displacement theory as shown
in Fig. <ref> for this implicit setup. The
normal vector of the undeformed configuration is obtained by the gradient
of the level-set function,

    N^⋆( X)=∇_ Xϕ( X)  for  X∈Γ_ X.



Let there be a displacement field u( X) which
lives in the full d-dimensional space (instead of only the manifold
itself as for parametric manifolds) so that the classical gradient
∇_ X u( X) is available. The resulting
deformation gradient is

    F_Ω=∇_ X x( X)= I+∇_ X u( X)

which is different from the surface deformation gradient F_Γ
in Eq. (<ref>). Based on this, one
may compute the normal vector of the deformed configuration at x= X+ u( X)
as

    n^⋆( x)=∇_ xϕ( X( x))= F_Ω^-T· N^⋆  for  x∈Γ_ x,

which follows by the chain rule. Eqs. (<ref>)
and (<ref>) are used to compute the projectors
P( X) and p( x), respectively.
The surface gradient (with respect to the undeformed configuration)
of a scalar function f( X) with X∈ℝ^d
results as

    ∇_ X^Γf= P·∇_ Xf.

As before, ∇_ Xf is the classical gradient in the d-dimensional
space. The situation is analogous for each component u_i of a
vector function u( X), so that one obtains
for the directional surface gradient

    ∇_ X^Γ,dir u     =    ∇_ X u· P,
    for  u∈ℝ^3:[[ ∂_X^Γu ∂_Y^Γu ∂_Z^Γu; ∂_X^Γv ∂_Y^Γv ∂_Z^Γv; ∂_X^Γw ∂_Y^Γw ∂_Z^Γw ]]     =    [[ ∂_Xu ∂_Yu ∂_Zu; ∂_Xv ∂_Yv ∂_Zv; ∂_Xw ∂_Yw ∂_Zw ]]·[[ P_11 P_12 P_13; P_12 P_22 P_23; P_13 P_23 P_33 ]].

The surface deformation gradient F_Γ follows
using Eq. (<ref>). The stretch of
a differential element upon the deformation is 

    Λ=‖∇_ xϕ‖/‖∇_ Xϕ‖· F_Ω=‖ n^⋆‖/‖ N^⋆‖· F_Ω.

Finally, for the operator W relating surface gradients of
the undeformed and deformed configuration, one obtains

    ∇_ x^Γf= W·∇_ X^Γf  with  W= p· F_Ω^-T.

Note that W· P= W, hence, when using the classical
derivatives one obtains ∇_ x^Γf= W·∇_ Xf.




  §.§.§ Manifolds with codimension 2


The focus is on one-dimensional manifolds in ℝ^3 such
as cables in the three-dimensional space. As mentioned before, two
level-set functions ϕ_1 and ϕ_2 are needed for the
geometry definition,

    Γ_ X={ X∈Ω_ X:ϕ_1( X)=0 and ϕ_2( X)=0} ,

see Fig. <ref>. The d-dimensional displacement
field u is given as before with classical derivatives ∇_ X u( X)
and the related deformation gradient F_Ω as in Eq. (<ref>).




The two normal vectors associated to the deformed and undeformed configuration
each, are given for i={ 1,2} as

    N_i^⋆    =    ∇_ Xϕ_i( X)  for  X∈Γ_ X,
     n_i^⋆    =    ∇_ xϕ_i( X( x))= F_Ω^-T· N_i^⋆  for  x∈Γ_ x.

One may then compute the unique tangent vectors 

    T^⋆    =     N_2^⋆× N_1^⋆,
     t^⋆    =     n_2^⋆× n_1^⋆.

The projectors P and p follow using Eqs. (<ref>)
and (<ref>), respectively. It is noted that the
same projector P is obtained when using N_1^⋆
and the orthogonalized normal vector N_3^⋆= N_2^⋆- N_1^⋆· N_2^⋆/ N_1^⋆· N_1^⋆· N_1^⋆
and

    P= I- N_1⊗ N_1- N_3⊗ N_3  with  N_1= N_1^⋆/‖ N_1^⋆‖, N_3= N_3^⋆/‖ N_3^⋆‖,

analogously for p. As before, these projectors are used to
determine surface gradients of scalar and vector functions as in Eqs. (<ref>)
and (<ref>). The line stretch is given
as

    Λ=‖ t^⋆‖/‖ T^⋆‖· F_Ω=‖ n_2^⋆× n_1^⋆‖/‖ N_2^⋆× N_1^⋆‖· F_Ω.






 §.§ Similarities and differences in the parametric and implicit descriptions


In this section, purely based on geometric considerations related
to the undeformed and deformed situation, a number of useful quantities
and operators are given following the concept of the TDC. The situation
is summarized in Table <ref> for parametric
manifolds and in Table <ref> for implicit ones.
For the mathematical equivalence of these two descriptions and more
details, see, e.g., <cit.>. The order of the rows in
the tables is determined by the information available for parametric
and implicit manifolds. Thereby, it is made sure that the (geometric
and differential) quantities of interest may be computed in this order.







It was found above that for parametrized manifolds, tangent vectors
result naturally as primary quantities through the existence of Jacobi
matrices. For problems with codimension 1, one may then obtain
normal vectors as secondary quantities (e.g., by a cross product of
tangent vectors), however, for higher codimensions, no unique normal
vector exists. For implicit manifolds, the situation is rather the
opposite: Here, normal vectors result naturally through the gradients
of level-set functions. A unique tangent vector is only applicable
for one-dimensional manifolds and may, for example, be computed by
a cross product of the normal vectors. This situation is an example
of the duality in the parametric and implicit description of manifolds.




 §.§ Further definitions


Based on the previous definitions, some additional differential operators
are introduced which apply to parametric and implicit manifolds equivalently. 




  §.§.§ Covariant surface gradient and divergence


The covariant surface gradient of a vector function u( X):Γ_ X→ℝ^d
is based on the projection of the directional one (see Eqs. (<ref>)
and (<ref>) for the parametric and implicit
situation, respectively) onto the tangent space. With respect to the
undeformed configuration, it is defined as 

    ∇_ X^Γ,cov u( X)= P·∇_ X^Γ,dir u( X).

Concerning the surface divergence of vector functions u( X)
and tensor functions A( X):Γ_ X→ℝ^d× d,
there holds 

    Div_Γ  u( X)     =    tr(∇_ X^Γ,dir u)=tr(∇_ X^Γ,cov u)∇_ X^Γ· u,
    Div_Γ  A( X)     =    [[ Div_Γ(A_11,A_12,A_13); Div_Γ(A_21,A_22,A_23); Div_Γ(A_31,A_32,A_33) ]]∇_ X^Γ· A.

The divergence operator with respect to the deformed configuration
follows accordingly as div_Γ  u( X)=∇_ x^Γ· u
and div_Γ  A( X)=∇_ x^Γ· A,
respectively.




  §.§.§ Conormal vector on the boundary


Unit normal and tangent vectors on manifolds have already been used
before and exist with respect to the deformed and undeformed configuration,
respectively. For example, the unit normal vector on an undeformed
membrane is N( X) for all X∈Γ_ X
and, after the deformation, n( x) for all
x∈Γ_ x. It is important to note that, for physical
reasons, the manifolds used in this work are bounded. The boundary
of the undeformed configuration is labeled ∂Γ_ X
and in the deformed situation ∂Γ_ x.




There exists a conormal unit vector N_∂Γ
along the boundary ∂Γ_ X which is in the tangent
plane of the manifold yet normal to ∂Γ_ X. This
vector points out of the manifold in the direction which naturally
extends the manifold, see Fig. <ref>. The
computation of the conormal vectors is straightforward (often using
cross products) and depends on q and d. In the deformed configuration,
the situation is similar for computing n_∂Γ( x)
for all x∈∂Γ_ x. In the context of the
definition of boundary value problems on manifolds, the conormal vectors
play a crucial role for the consideration of boundary conditions as
shall be seen in Section <ref>.




  §.§.§ Divergence theorem on manifolds


To derive the weak form of the governing equations later on, the following
divergence theorem on manifolds is needed <cit.>,

    ∫_Γ_ X u·Div_Γ A dΓ=-∫_Γ_ X∇_ X^Γ,dir u: A dΓ+∫_Γ_ Xϰ· u· A· N dΓ+∫_∂Γ_ X u· A· N_∂Γ d∂Γ,

where ∇_ X^Γ,dir u: A=tr(∇_ X^Γ,dir u· A^T)
is a matrix scalar product. The mean curvature is ϰ=tr( H)
with H=∇_ X^Γ,dir N=∇_ X^Γ,cov N
being the second fundamental form. For in-plane tensor functions
with A= P· A· P, the term involving the
curvature ϰ vanishes and one finds ∇_ X^Γ,dir u: A=∇_ X^Γ,cov u: A.




§ MECHANICAL MODEL AND GOVERNING EQUATIONS


In Section <ref>, a number of geometric quantities (such
as normal vectors, projectors, area/line stretches etc.) and differential
operators related to (surface) gradients are introduced. It was shown
how these quantities are obtained for parametrized and implicitly
defined manifolds. The focus is now turned to the mechanics and the
procedure follows the classical outline, however, it is based on the
TDC here. It is emphasized that all tensors considered in the
following refer to the parametric as well as implicit situation. These
tensors have dimensions d× d (with d being the dimension
into which the cable or membrane is immersed). A tensor A
is called “in-plane” or “tangential” to the undeformed configuration
Γ_ X if A= P· A· P and to
the deformed configuration Γ_ x if A= p· A· p.
An in-plane (d× d)-tensor has q non-zero eigenvalues
representing the principal mechanical quantity (with q being the
dimension of the manifold: q=1 for cables, q=2 for membranes).

Starting point is the surface deformation gradient F_Γ( X)
at X∈Γ_ X, specified previously in Eq. (<ref>).
It may also be seen as a geometrical quantity mapping tangent vectors
from the undeformed to the deformed configuration. It is also noted
that the situation also applies to the “volumetric” case (where
q=d=2 are flat shells in two dimensions and q=d=3 are volumetric
continua in three dimensions). In this case, F_Γ= F_Ω
as in Eq. (<ref>), and for the projectors,
P= p= I. In that sense, the presented mechanical outline
below applies to cables, membranes and continua in a unified sense.




 §.§ Governing equations in strong form





  §.§.§ Strain tensors


Based on the surface deformation gradient, the directional and tangential
Cauchy-Green strain tensors are defined as 

    E_dir    =    1/2·( F_Γ^T· F_Γ- I),
     E_tang    =     P· E_dir· P,

respectively. The Euler-Almansi strain tensors are

    e_dir    =    1/2·( I-( F_Γ· F_Γ^T)^-1),
     e_tang    =     p· e_dir· p,

where e_tang is tangential to the deformed configuration
Γ_ x. As usual, there holds e_dir= F_Γ^-T· E_dir· F_Γ^-1
(which is not true for the tangential versions of these strain tensors).




  §.§.§ Stress tensors


Conjugated stress tensors are introduced next and only the tangential
versions are considered. Generally speaking, we assume some hyper-elastic
material with an elastic energy function Ψ( E_tang)
and obtain the second Piola-Kirchhoff stress tensor as S=∂Ψ/∂ E_tang.
For simplicity, only Saint Venant–Kirchhoff solids are considered
herein and there follows

    S     =    λ·trace( E_tang)· P+2μ E_tang,
        =     P·(λ·trace( E_dir)· I+2μ E_dir)· P,

with S being tangential to Γ_ X. λ
and μ are the Lamé constants; for cables λ becomes
0. On the other hand, the Cauchy stress tensor reads

    σ=1/Λ· F_Γ· S· F_Γ^T,

where Λ is a line stretch for cables and an area stretch
for membranes when undergoing the displacement, see Section <ref>.
For the volumetric case (q=d), Λ= F_Ω is
the volumetric stretch. The Cauchy stress is tangential to the deformed
configuration Γ_ x since F_Γ· P= p· F_Γ· P
and P· F_Γ^T= P· F_Γ^T· p,
hence σ= p·σ· p. The first Piola-Kirchhoff
stress tensor is given by

    K= F_Γ· S

and there holds K= K· P= p· K.




  §.§.§ Relation of stress and strain tensors


For every point X∈Γ_ X and its mapped counterpart
x( X)∈Γ_ x, we have the following
equality,

    S( X): E_tang( X)=(σ( x): e_tang( x))·Λ( X),

where : represents the matrix scalar product operator. In this
sense the two stress tensors S and σ are conjugated
to their related strain tensors E_tang and e_tang,
respectively. It is noted that

    S: E_tang= S: E_dir  and  σ: e_tang=σ: e_dir

which will be important later. Furthermore, the result of these matrix
scalar products may also be derived by the non-zero eigenvalues S_i,
E_tang,i, σ_i, e_tang,i, i=1,…,q,
of the tangential tensors S, E_tang, σ,
e_tang, respectively. Hence, we obtain

    S: E_tang=∑_i=1^qS_i· E_tang,i  and  σ: e_tang=∑_i=1^qσ_i· e_tang,i.






  §.§.§ Equilibrium


A crucial aspect of finite strain theory is that equilibrium is to
be fulfilled in the deformed configuration which is expressed
in strong form as

    div_Γ σ( x)=- f( x)  ∀ x∈Γ_ x,

where f are body forces. Recall from (<ref>)
that div_Γ σ=∇_ x^Γ,dir·σ=∇_ x^Γ,cov·σ
is the divergence of the Cauchy stress tensor with respect to Γ_ x.
Furthermore, we have the identity

    Div_Γ  K( X)=div_Γ σ( x( X))·Λ( X)

with Div_Γ  K=∇_ X^Γ,dir· K=∇_ X^Γ,cov· K
being the divergence of the first Piola-Kirchhoff stress tensor from
Eq. (<ref>) with respect to Γ_ X.
In order to transform the derivatives in the divergence operators
from the undeformed to the deformed situation, use Eqs. (<ref>)
and (<ref>) for parametric and implicit
manifolds, respectively. Due to F( X)= f( x( X))·Λ( X),
the equilibrium in Γ_ x can be stated equivalently to
Eq. (<ref>) based on quantities in the undeformed
configuration as

    Div_Γ  K( X)=- F( X)  ∀ X∈Γ_ X.






  §.§.§ Boundary conditions


The domain of interest is a bounded manifold where the boundary ∂Γ
falls into the two non-overlapping parts ∂Γ_D
and ∂Γ_N, which holds in the deformed and
undeformed configuration Γ_ X and Γ_ x,
respectively. Hence, the boundary conditions in the deformed configuration
are 

    u( x)     =ĝ( x) on ∂Γ_ x,D,
    σ( x)· n_∂Γ( x)     =ĥ( x) on ∂Γ_ x,N,

where ĝ are prescribed displacements and ĥ
are tractions (force per area for q=3, force per length for q=2
or a single force for q=1). Note that for ropes and membranes,
ĥ must be in the tangent space of the deformed manifold
in order to satisfy the equilibrium due to the absence of bending
stresses or transverse shear stresses. The equivalent boundary conditions
formulated in the undeformed configuration are 

    u( X)     =Ĝ( X) on ∂Γ_ X,D,
     K( X)· N_∂Γ( X)     =Ĥ( X) on ∂Γ_ X,N ,

where Ĝ and Ĥ have similar interpretations
as before. The relation between ĥ and Ĥ
is Ĥ(X)=Λ̅(X)·ĥ(x),
where 

    Λ̅=
    1    for q=1, d={ 2,3} (cables),
    line stretch along the boundary   for q=2, d=23 (membranes),
    area stretch of the face at the boundary   for q=3, d=3 (continuum).

Further information about boundary conditions for ropes and membranes
are given in <cit.>.

With the boundary conditions above, the complete second-order boundary
value problem (BVP) is defined in the deformed and undeformed configuration.
The obtained BVP in the frame of the TDC is valid for explicitly and
implicitly defined manifolds and does not rely on curvilinear coordinates
implied by a parametrization, which are typically used in classical
approaches, see, e.g., <cit.>. Therefore,
the proposed formulation of ropes and membranes based on the TDC is
more general compared to the classical theory.




 §.§ Governing equations in weak form


For stating the governing equations in weak form, the following test
and trial function spaces are introduced 

    𝒮_ u    ={ v∈[ℋ^1(Γ_ X)]^d:  v=Ĝ on ∂Γ_ X,D} ,
    𝒱_ u    ={ v∈[ℋ^1(Γ_ X)]^d:  v=0 on ∂Γ_ X,D} ,

where ℋ^1 is the Sobolev space of functions with square
integrable first derivatives. The task is to find u∈𝒮_ u
such that for all w∈𝒱_ u, there holds 

    η·∫_Γ_ X∇_ X^Γ,dir w: K( u) dΓ=η·∫_Γ_ X w· F dΓ+∫_∂Γ_ X,N w·Ĥ d∂Γ.

where 

    η=
    A    for q=1, d={ 2,3} is the cross section of the cable,
    
    t    for q=2, d={ 2,3} is the thickness of the shell/membrane,
    
    1    for q=3, d=3, i.e., a continuum.

The integrals in Eq. (<ref>) are one-, two-, or
three-dimensional for cables, membranes and continua, respectively.
The multiplication with η ensures that the units of the integration
are always consistent. Hence, it is possible to naturally consider
situations where cables, membranes and continua are coupled in one
setup by simply adding up the corresponding integrals as in Eq. (<ref>)
for each structure. In order to obtain Eq. (<ref>),
we applied the usual procedure for converting the strong form to a
weak form: Multiply Eq. (<ref>) with test
functions, integrate over the domain Γ_ X, and apply
the divergence theorem from Eq. (<ref>). It is noteworthy
that the curvature term from Eq. (<ref>) vanishes
also for cables and membranes due to K· N=0.

The weak form stated above is related to energy minimization in the
sense that

    ∫_Γ_ X∇_ X^Γ,dir w: K( u) dΓ=∫_Γ_ Xδ E_tang( u): S( u) dΓ,

where δ is the variational operator.




  §.§.§ Energy relation


An immediate consequence of Eq. (<ref>) is that
one may obtain the same stored potential energy of the deformed body
by integrating over the undeformed or deformed configuration as follows

    𝔢( u)     =    1/2η·∫_Γ_ x e_tang( u):σ( u) dΓ,
        =    1/2η·∫_Γ_ X E_tang( u): S( u) dΓ.






§ DISCRETIZATION AND NUMERICAL METHODS


In order to solve the boundary value problem, i.e., to approximate
the sought displacements, one may use two fundamentally different
approaches. Possibly the more intuitive one is to discretize the cable(s)
or membrane(s) by curved (q-dimensional) line or surface elements,
respectively. Then, each element results from an (often isoparametric)
map of some reference element so that this approach is naturally linked
to the parametric description of manifolds as discussed in Section
<ref>. This is the classical approach labeled
Surface FEM herein. An alternative is to use a (d-dimensional)
background mesh for the approximation of the weak form. That
is, higher-dimensional shape functions (than the dimension of the
manifold) are used and evaluated on the trace of the manifold only.
This approach is called Trace FEM or Cut FEM. It is naturally related
to an implicit manifold description as discussed in Section <ref>.

It is important to note that the classical definition of finite strain
theory based on curvilinear coordinates does not cover the latter
approach. This is another reason why the presented formulation based
on the TDC is more general as it supports both, the Surface and Trace
FEM. For the discussion below, we assume the manifold case, hence,
q<d. With q=d, the situation results in the standard FEM which
is not further outlined here.




 §.§ Surface FEM





Starting point is the discretization Γ_ X^h of the
undeformed cable (q=1) or membrane (q=2) by a line or surface
mesh, respectively. Herein, we use higher-order q-dimensional Lagrange
elements with equally spaced nodes in the reference element. The nodal
coordinates in the undeformed configuration are labeled X_i
with i=1,…,n_q and n_q being the number of nodes in
the mesh, see Fig. <ref>. The resulting shape functions
M_i^q( X) span a C^0-continuous finite
element space as 

    𝒬_Γ_ X^h:={ v_h∈ C^0(Γ_ X^h): v_h=∑_i=1^n_qM_i^q( X)·v̂_i with v̂_i∈ℝ}⊂ℋ^1(Γ_ X^h) .

M_i^q( X) are obtained by isoparametric mappings from the
q-dimensional reference element to the physical elements in d
dimensions. Based on Eq. (<ref>), the following
discrete test and trial function spaces are introduced 

    𝒮_Γ_ X^h    ={ v_h∈[𝒬_Γ_ X^h]^d:  v_h=Ĝ on ∂Γ_ X,D^h} ,
    𝒱_Γ_ X^h    ={ v_h∈[𝒬_Γ_ X^h]^d:  v_h=0 on ∂Γ_ X,D^h} .

The discrete weak form of Eq. (<ref>) reads as follows:
Given Lamé constants (λ,μ)∈ℝ^+, body forces
F∈ℝ^d on Γ_ X^h, tractions Ĥ∈ℝ^d
on ∂Γ_ X,N^h, find the displacement
field u_h∈𝒮_Γ_ X^h such that
for all test functions w_h∈𝒱_Γ_ X^h
there holds in Γ_ X^h

    η·∫_Γ_ X^h∇_ X^Γ,dir w_h: K( u_h) dΓ=η·∫_Γ_ X^h w_h· F dΓ+∫_∂Γ_ X,N^h w_h·Ĥ d∂Γ.

The sought discrete displacement field u_h( X) is obtained
solving a non-linear system of equations for the n_DOF=d· n_q
nodal values (degrees of freedom) as usual in the context of finite
strain theory.




 §.§ Trace FEM








Let there be a d-dimensional background mesh into which the manifold
is completely immersed. Only those elements and corresponding nodes
are considered that are cut by the manifold, see Fig. <ref>.
They may be labeled “active” elements and nodes, all others are
neglected. The shape functions of the active nodes are constructed
by (often isoparametric) mappings from a d-dimensional reference
element, but the shape functions are only evaluated on the q-dimensional
manifold.

The Trace FEM is a fictitious domain method (FDM) for PDEs on manifolds
<cit.>.
As in any FDM, there is no boundary-conforming mesh but a background
mesh, herein further complicated by the fact that background mesh
and manifold have different dimensions. In general, the following
issues have to be properly addressed:


  * Integration points have to be defined for the integration of the weak
form of the governing equations—only at these points, shape functions
are evaluated. This requires the identification of the zero-level
set of some level-set function ϕ( X) within each
active background element, see Fig. <ref>.
The situation may be further complicated, when the boundary of the
manifold is within the background element, which may be defined by
additional slave level-set functions ψ_i( X)
as mentioned in Section <ref>. The placement
of integration points is an important and challenging task, in particular
with higher-order accuracy. For an overview, we refer to the references
<cit.>. It is
useful to evaluate given level-set functions at the active nodes and
interpolate them in-between. The identification of the zero-level
sets and placement of integration points may then be achieved in the
d-dimensional reference element which simplifies the evaluation
of shape functions in the background elements. Of course, the interpolated
zero-level set is only a (higher-order) approximation of the exact
manifolds Γ_ X and Γ_ x and labeled Γ_ X^h
and Γ_ x^h , respectively. In the Trace FEM context,
Γ_ X^h and Γ_ x^h may be seen as integration
cells, see Fig. <ref>(b), to define integration
points, they do not imply shape functions.

  * The treatment of boundary conditions is a challenging task in FDMs
as it is impossible to directly prescribe values of the nodes in the
active background elements. The additional constraints may, in principle,
be enforced using penalty methods, Lagrange multiplier methods, or
Nitsche's method. The latter has been developed to be a standard choice
in FDMs because the equations are formulated in a consistent way without
needing additional degrees of freedom <cit.>. In the
standard (symmetric) form of the Nitsche's method, stabilization parameters
are required. In the simplest case, these parameters may be set to
a fixed user-defined number <cit.>, however,
for background elements which are cut by a tiny fraction of the manifold,
resulting in small supports of the shape functions, this may lead
to unsatisfactory results. An alternative is to compute the stabilization
parameters based on global or element-wise generalized eigenvalue
problems <cit.>. However, for the awkward
cut situations, this may result in unbounded values resulting in similar
problems known for penalty methods <cit.>. Therefore,
we prefer the non-symmetric Nitsche's method herein with the
main advantage that an additional stabilization is not required
for imposing boundary conditions <cit.>.

  * Stabilization is still necessary when applying FDMs in the context
of PDEs on manifolds to ensure the regularity of the resulting system
of equations. This may be traced back to two sources: One is found
in the shape functions with small supports and the other in the fact
that the approximated displacement field on the manifold may not necessarily
be represented by a unique set of nodal values in the background
mesh. That is, the background shape functions restricted to the manifold
build a frame but not necessarily a basis <cit.>.
Fortunately, different stabilization approaches exist to cure both
issues and we refer to the overview given in <cit.>
for the Trace FEM. Herein, we use the “normal derivative volume
stabilization”, introduced for scalar-valued problems in <cit.>
and in <cit.> for vector-valued problems. This stabilization
technique enables higher-order accurate results, does not change the
sparsity pattern of the stiffness matrix, and only first derivatives
are needed.

With these comments made, we are ready to define the discrete weak
form for the Trace FEM. Let Ω_ X^h be the background
mesh into which the undeformed configuration Γ_ X is
immersed, only active elements and nodes are present in Ω_ X^h.
The resulting shape functions M_i^d( X) span
a C_0-continuous finite element space on the manifold as 

    𝒬_Ω_ X^h:={ v_h∈ C^0(Ω_ X^h): v_h=∑_i=1^n_dM_i^d( X)·v̂_i with v̂_i∈ℝ}⊂ℋ^1(Ω_ X^h) .

Based on Eq. (<ref>), the following discrete
trace test and trial function spaces are introduced 

    𝒯_Γ_ X^h    ={ v_h∈[𝒬_Ω_ X^h]^d:  v_h only on Γ_ X^h} ,
    𝒰_Γ_ X^h    =𝒯_Γ_ X^h.

For the Trace FEM, the discrete weak form of Eq. (<ref>)
is: Given Lamé constants (λ,μ)∈ℝ^+, body
forces F∈ℝ^d on Γ_ X^h, tractions
Ĥ∈ℝ^d on ∂Γ_ X,N^h,
and stabilization parameter ρ, find the displacement field u_h∈𝒯_Γ_ X^h
such that for all test functions w_h∈𝒰_Γ_ X^h
there holds in Γ_ X^h

    η·∫_Γ_ X^h∇_ X^Γ,dir w_h: K( u_h) dΓboundary term due to  on ∂Γ_ X,D^h-∫_∂Γ_X,D^h w_h·[ K( u_h)· N_∂Γ_ X] d∂Γ+
    Nitsche term∫_∂Γ_ X,D^h( u_h-Ĝ)·[ K( w_h)· N_∂Γ_ X] d∂Γ+stabilization termρ∫_Ω_ X^h( N^e·∇_ X u_h)·( N^e·∇_ X w_h) dΩ
    
    =η·∫_Γ_ X^h w_h· F dΓ+∫_∂Γ_ X,N^h w_h·Ĥ d∂Γ .

In comparison to Eq. (<ref>), additional terms
occur in the discrete weak form due to the weak enforcement of essential
boundary conditions with Nitsche's method and the stabilization. The
sought discrete displacement field u_h( X) is obtained
solving a non-linear system of equations for the n_DOF=d· n_d
nodal values (degrees of freedom). 

Note that the stabilization term is the only term which is not evaluated
on the manifold Γ_ X^h but in the volumetric background
mesh Ω_ X^h (using standard Gauss integration). Therefore,
one has to extend the normal vector N( X)
from the undeformed manifold Γ_ X^h to the neighborhood,
resulting in N^e( X) for all X∈Ω_ X^h.
This is particularly simple for implicitly defined manifolds, e.g.,
using N^e( X)=∇_ Xϕ^h( X)/‖∇_ Xϕ^h( X)‖
for all X∈Ω_ X^h. Furthermore, in the stabilization
term, the classical gradient operator ∇_ X
is used instead of the surface operators used in the other
terms. In <cit.>, it is recommended that the stabilization
parameter should be chosen in the range O(h)≲ρ≲ O(h^-1),
where h is the element size in the active background mesh.

A remark is added concerning slip supports because the above
mentioned weak form rather expects that all displacement components
are prescribed through Ĝ along the Dirichlet boundary
∂Γ_ X,D^h. In the Nitsche's method,
displacement constraints in selected, arbitrary unit directions v_d
with magnitude Ĝ may be prescribed by replacing the corresponding
terms in Eq. (<ref>) with

    -∫_∂Γ_ X,D^h( w_h· v_d) [ K( u_h)· N_∂Γ_ X]· v_d d∂Γ
    
    +∫_∂Γ_ X,D^h( u_h· v_d-Ĝ) [ K( w_h)· N_∂Γ_ X]· v_d d∂Γ .






§ NUMERICAL RESULTS


A number of test cases for ropes and membranes in two and three dimensions
are considered in this section. The numerical results focus on the
convergence rates for two different types of errors. The “energy
error” ε_𝔢 compares the approximated stored
elastic energy with the analytical one,

    ε_𝔢=|𝔢( u)-𝔢( u_h)|,

with 𝔢 computed based on Eq. (<ref>).
The analytical energy 𝔢( u) may also
be computed by an overkill approximation, i.e., based on an extremely
fine mesh with higher-order elements. Provided that geometry and boundary
conditions allow for sufficiently smooth solutions, the expected convergence
rates in this error norm are p+1 with p being the order of the
elements.

The “residual error” ε_res integrates
the error in the equilibrium as stated in Eq. (<ref>),
that is,

    ε_res=√(∑_e=1^n_el∫_Γ_ X^h,e𝔯( u_h)·𝔯( u_h) dΓ)  with 𝔯( u_h)=div_Γ σ( u_h)+ f( x)

This error obviously vanishes for the analytical solution. It is important
to note that the integrand in (<ref>) involves second-order
derivatives. Therefore, the integral must not be carried out over
the whole (discretized) domain Γ_ X^h but integrated
element by element as indicated by the summation. That is, element
boundaries, where already the first derivatives of the C^0-continuous
shape functions feature jumps, are neglected in computing ε_res.
Due to the presence of second-order derivatives, the expected convergence
rates are p-1 which also indicates that higher-order elements are
crucial for convergence in ε_res.




 §.§ Membrane with given deformation


In the first test case, a membrane in the shape of a half sphere with
radius r=1.0 undergoes a prescribed displacement and the stored
elastic energy is computed from the viewpoint of the Surface FEM and
the Trace FEM. That is, in the Surface FEM, surface meshes
with different resolutions and element orders are generated. See Fig. <ref>(c)
for some example mesh composed by quadratic elements. The displacements

    u( X)=[[          5/2+1/5(X+1);                  -3/2; -1/2[1-(X^2+Y^2)]-3/2 ]]

are evaluated at the nodes and interpolated based on the shape functions
implied by the surface meshes, yielding u_h( X);
see Fig. <ref>(c) for the resulting deformed membrane.
Then, the elastic energy of the deformed configuration 𝔢( u_h)
is computed with Eq. (<ref>).

For the Trace FEM viewpoint, background meshes of different
resolutions and orders are generated in Ω_ X=[-1,1]×[-1,1]×[0,1]
and the geometry is defined based on the level-set function ϕ( X)=‖ X‖ -r.
See Fig. <ref>(c) for a sketch of the situation using
quadratic background elements. Then, the (active) nodes of the background
mesh are deformed by the given displacement field (<ref>),
yielding u_h( X) based on the shape functions
implied by the background meshes. This displacement field living
in the whole background mesh is only evaluated on the membrane surface
in order to compute the stored energy 𝔢( u_h)
according to the Trace FEM.




We set the Lamé constants to λ=3 and μ=2. The resulting
energy is given by the value 𝔢( u)=1.642871443585262.
In Fig. <ref>, the convergence results for the various
meshes are shown for the Surface and the Trace FEM. It is seen that
in both cases optimal convergence results are achieved. The
energy error converges one order higher than expected for even element
orders. In the Trace FEM, the convergence curves are less smooth than
in the Surface FEM because the approximation spaces are not nested
upon refinement; this is well-known for results obtained with FDMs
in general.

It is thus seen that the Surface FEM as well as the Trace FEM have
the potential to achieve optimal results. For all other test cases
below, we shall now obtain the discrete displacement fields based
on solving the non-linear systems of equations resulting from the
weak forms given in Section <ref>.




 §.§ Rope with Surface and Trace FEM





The second test case considers a rope in two dimensions as shown in
Fig. <ref>. The cross section of the rope is A=0.01,
Young's modulus is E=10 000, and the Lamé constants are μ=1/2E
and λ=0. The left support is located at (0,1/2)
and the right at (1,0). The geometry may be given in
parametric form for r∈(0,1) as 

    X(r)=[[ X(r); Y(r) ]]=[[                     r; 1/2(1-r)-1/7sin(π· r) ]].

Alternatively, the same geometry may be implied by the level-set function

    ϕ( X)=Y-[1/2(1-X)-1/7sin(π· X)]  ∀ X∈(0,1).

The structure is loaded by its own weight with F( X)=[0,-2 000· A]^T
for all X∈Γ_ X. The deformed rope is illustrated
in Fig. <ref>(a) where the color information indicates
the principal Cauchy stress in the rope. The stored elastic energy
is 𝔢=0.7528302283000 and the length of the rope is increased
by a factor of 1.1053648264108 due to the deformation.




The displacements are computed with the Surface and Trace FEM, respectively.
For the Trace FEM, the stabilization parameter in Eq. (<ref>)
is set to ρ=1000/h. The mesh resolutions and orders
are systematically varied and convergence results are seen in Fig. <ref>.
Figs. <ref>(a) and (b) show the error in the elastic
energy ε_𝔢 and Figs. <ref>(c)
and (d) the residual error ε_res. Of course,
for some given element length h, the use of (curved) line meshes
in the Surface FEM results in considerably less degrees of freedom
than using background meshes in the Trace FEM. Therefore, the convergence
studies for the Surface FEM are realized for up to 1024 elements
whereas the background meshes in the Trace FEM feature up to 320×160
elements (of which only those cut by the rope are active and, hence,
taken into account for the simulation). It is apparent from the convergence
results in Fig. <ref> that optimal higher-order
convergence rates are achieved. It is noteworthy that the Surface
FEM with even element orders converges one order higher than expected
in ε_𝔢, whereas this is not the case for
the Trace FEM. This can already be traced back to the accuracy in
the numerical integration: Simply integrating the length of the rope
using the Surface or Trace FEM viewpoint shows that an extra order
in the accuracy is achieved with the Surface FEM for even element
orders. For the convergence results in ε_res,
results for linear meshes are omitted because second order derivatives
are needed for this error measure, see Eq. (<ref>).
In this error norm, Surface and Trace FEM converge optimally with
p-1.




 §.§ Membranes with Surface FEM





For the next test case, the deformation of membranes is approximated
with the Surface FEM. Two different undeformed configurations are
considered which result from a map of a unit circle for map
A and of a unit square for map B. The undeformed configurations
Γ_ X are given by the parametrizations

    X( r)=[[      3/2· r;           s; c·sin(r· s) ]] with √(r^2+s^2)∈(0,1)    for map A,
    
    r, s∈(-1,1)    for map B,

where c∈ℝ is a scaling parameter in vertical direction.
The resulting configurations for map A and B are illustrated in Figs. <ref>(a)
and (d) for c=0.4, respectively. An important difference is that
the boundary is smooth for map A but involves corners for map B, later
resulting in different convergence behaviors. The thickness of the
membrane is t=0.01, Young's modulus is E=1 000 and Poisson
ratio ν=0.3, which is easily converted into the Lamé parameters.
The loading is gravity acting on the membrane surface with F( X)=[0,0,-200· t]^T
for all X∈Γ_ X. The whole boundary is treated
as a Dirichlet boundary with prescribed zero-displacements. The deformed
configurations are displayed in Figs. <ref>(b)
and (e) with computed von-Mises stresses based on the Cauchy stress
tensor. The vertical displacement field is given in Fig. <ref>(c)
and (f) for the two maps in top view, respectively.




Convergence results are given in Fig. <ref>. For
map A, where the boundary is smooth, optimal convergence rates are
found in the energy error ε_𝔢 and residual
error ε_res. For map B, where corners are
present in the membrane geometry, it is seen that the convergence
rates are bounded. Only linear and quadratic elements converge optimally
in ε_𝔢, higher orders improve the error
level, however, not the convergence rates. It is thus confirmed that
corners in membranes have the potential to reduce the convergence
rates.




 §.§ Coupled ropes and membranes





It was mentioned several times that the proposed framework allows
for a unified treatment of ropes and membranes (and even continua).
This shall be confirmed with this last test case where cables and
membranes are coupled. The situation may be described as an inflated
ball with embedded reinforcement cables, see Fig. <ref>(a).
The original radius of the ball is r=1.0 and the inner pressure
is p=20. This pressure converts to a loading of f( x)=p· n( x)
which depends on the normal vector of the deformed configuration,
hence on the displacements u, thereby adding a new source
of non-linearity which has to be properly considered in the Newton-Raphson
loop. The membrane is defined by the parameters t=0.01, E=1 000
and ν=0.3. The cables which are on the gray planes in Fig. <ref>(a)
feature a Young's modulus of E=1 000 000 and a cross section
of A=0.0001. All other cables have a Young's modulus of E=500 000.
For the modeling, only one eighth of the initial sphere is considered,
see Fig. <ref>(b) for the deformed and undeformed
situations. Symmetry boundary conditions are applied. The resulting
von-Mises stresses in the membrane are shown in Fig. <ref>(c).
The elastic energy stored in the membrane and ropes is given
as 𝔢( u)=2.9802127651.




The convergence results are seen in Fig. <ref> for
the energy error ε_𝔢 and residual error
ε_res. It is again seen that the convergence
rates in ε_𝔢 are optimal for linear and
quadratic elements. Higher orders achieve better results, however,
do not improve the convergence rates. The reason for the bounded convergence
rates is found in the stress concentrations where the embedded cables
meet, see Fig. <ref>(c).




§ CONCLUSIONS


The modeling of ropes and membranes leads to partial differential
equations on manifolds. Herein, a framework is proposed which unifies
the mechanical modeling of ropes and membranes undergoing large displacements
and, furthermore, applies to parametric and implicit definitions of
manifolds. The fundamental ingredient is the use of Tangential Differential
Calculus (TDC) for the definition of geometric and differential quantities
on the manifolds. The proposed TDC-based formulation is more general
than the classical parametric formulation based on curvilinear coordinates
because it allows for two different numerical approaches, the Surface
and the Trace FEM. The classical approach is restricted to the Surface
FEM and cannot handle manifolds which are defined implicitly (unless
discretized by a surface mesh).

Even in the classical setting using the Surface FEM with parametric
formulations, the new TDC-based formulation leads to a significantly
different implementation, however, ultimately achieving the same results
(as expected). The advantage of the TDC-based formulation is that
it also applies immediately to implicit definitions and the Trace
FEM. More technically speaking, we have one element integration routine
which applies to ropes and membranes (and even continua) no matter
whether we are using the Surface FEM or the Trace FEM. Of course,
in the Trace FEM, additional terms have to be considered for the stabilization.

The numerical results consider various test cases with ropes and membranes
in two and three dimensions. The Surface and Trace FEM with higher-order
elements are used successfully, achieving higher-order convergence
rates. It is confirmed that the smoothness of the geometry and the
coupling of ropes and membranes has an important influence on the
convergence behavior.

schanz
tocsection
 
