A nonlinear rotation-free shell formulation with prestressing for vascular biomechanics

We implement a nonlinear rotation-free shell formulation capable of handling large deformations for applications in vascular biomechanics. The formulation employs a previously reported shell element that calculates both the membrane and bending behavior via displacement degrees of freedom for a triangular element. The thickness stretch is statically condensed to enforce vessel wall incompressibility via a plane stress condition. Consequently, the formulation allows incorporation of appropriate 3D constitutive material models. We also incorporate external tissue support conditions to model the effect of surrounding tissue. We present theoretical and variational details of the formulation and verify our implementation against axisymmetric results and literature data. We also adapt a previously reported prestress methodology to identify the unloaded configuration corresponding to the medically imaged in vivo vessel geometry. We verify the prestress methodology in an idealized bifurcation model and demonstrate the significance of including prestress. Lastly, we demonstrate the robustness of our formulation via its application to mouse-specific models of arterial mechanics using an experimentally informed four-fiber constitutive model.


A Axisymmetric theory
In order to validate the proposed framework, it is useful to recall the well-known thick wall cylindrical tube solution. This solution provides, in the most general case, the static mechanical response of a thick axisymmetric tube under combined inflation, extension, and torsion [1,2]. Here, we focus on inflation and extension. The vessel is assumed to be an incompressible thickwalled cylindrical tube that undergoes a finite deformation from an undeformed configuration Ω 0 to the current configuration Ω. For the comparison purposes of the example in section 5.1, no residual stress will be considered. Considering all of the above, configuration Ω 0 is given in terms of cylindrical coordinates (R, Θ, Z) as with R i , R o and L the initial inner and outer radii and initial length of the tube. Similarly, the cylindrical coordinates (r, θ , z) define the current configuration Ω as with r i , r o and l the current inner and outer radii and current length of the tube. In the absence of torsion and bending, the current cylindrical coordinates can be written in terms of the initial coordinates as with λ z the axial stretch. Additionally, the radial and circumferential stretches, λ r , λ θ are defined as The three stretches are related through the incompressibility condition as and allow defining the Green Lagrange strain tensor components as The equilibrium equation, in the absence of body forces, can be written in cylindrical coordinates as Taking into account that the tube is loaded by an internal pressure p i = −σ rr (r = r i ) and that no load is applied on the outer wall, i.e. σ rr (r = r o ) = 0, the above equation can be integrated across the thickness to yield with the stress components σ θ θ , σ rr obtained via some constitutive relationship, i.e.
(S10) Equation (S8), combined with the constitutive relationships in equations (S9) and (S10), yields a nonlinear equation that allows computation of the current internal radius r i as a function of the loading pressure p i and axial stretch λ z . The constitutive relationship in cylindrical coordinates for an incompressible material can be derived using a Lagrange multiplier, as described in Section 2.4 or by enforcing incompressibility directly into the strain energy function via equation (S5). Following the approach described in [1] , a general 3D SEF ψ el can be written in terms of only two principal strains, i.e.
whereψ el is the equivalent of the strain energy function that only depends on E ΘΘ , E RR . Usingψ el and assuming that the principal components of stress and strain coincide, the following is obtained (see [1] for details) which can be used in the nonlinear equation (S8). The partial derivative ofψ el respect to E ΘΘ can be found using the chain rule and equations (S5) and (S6) for an appropriate strain energy function.

B Variation of membrane strain and curvature matrices B.1 Variation of membrane strain
Using the assumed strain approach proposed by Oñate and Flores [3], the membrane strain field is considered as a linear interpolation of the values at the three edge midpoints of the master element. Consequently, for a reduced integration strategy with a single in-plane quadrature point per element, the variation of the Green-Lagrange membrane strain matrix can be expressed as where δ ϕ p denotes the 18 × 1 vector gathering the variation of the positions of the 6 nodes belonging to the element patch and B p m denotes the 3 × 18 matrix given as with (S17)

B.2 Variation of curvature
The variation of curvature is given as Therefore, the first term on RHS of equation (S18) can be written as ∑ a=1 L I ,α N I a,β + L I ,β N I a,α (a 3 · δ ϕ a ). (S19)

3/7
Accordingly, the second term on RHS of equation (S18) whereȃ α are contravariant basis vectors in the master element defined using the linear interpolation and obtained as with a M αβ = ϕ M ,α · ϕ M ,β . Combining these expressions, the variation of curvature is given as (S23)

C Time integration and linearization
We employ a Generalized-α time integration scheme [4] to solve the linearized equation of linear momentum balance. To this end, we consider variations of δW kin , δW int , and δW ext described in section 3 with respect to the discretized displacement vector to obtain the vectors of kinetic, internal, and external nodal forces as where the partial derivatives of E m and χ are obtained, respectively, from equations (S14) and (S23). Using the results above, and upon assembly, the discrete nonlinear equation for linear momentum balance can be written as where R is the assembled residual vector, F ext , F kin and F int are, respectively, the assembled external, kinetic and internal force vectors andÜ,U, U are, respectively, the global vector of nodal accelerations, velocities and displacements. The nonlinear, transient system of equation (S27) is solved using the generalized alpha method, where the residual is evaluated as with α M , α F the time integrator parameters (see reference [4] for details). The above system of equations is solved using Newton-Raphson technique, in the form

4/7
with ∆t the time step size, γ, β the Newmark time integration parameters, M the mass matrix, C the damping matrix and K the stiffness matrix. We will not provide details in the derivation of the mass matrix (standard in Finite Element) nor in the derivation of the damping matrix (the contribution of which comes from the external tissue support). Regarding the stiffness matrix, the elemental contribution (before assembly) for a couple of nodes a, b is given as as a 3 × 3 matrix: and introduce the following thickness integrations of this matrix Using Voigt notation to express the components of various tensors as the above expressions allow us to write the internal component of tangential stiffness matrix as Here, the terms in the first row represent the material component of the stiffness matrix while the terms in the second row represent the geometric component of the stiffness matrix.