A NURBS-based Inverse Analysis for Reconstruction of Nonlinear Deformations of Thin Shell Structures

TL;DR

This paper introduces a NURBS-based inverse analysis method that accurately reconstructs nonlinear deformations and loads of thin shell structures, capturing buckling and snapping phenomena for improved manufacturing processes.

Contribution

It combines NURBS-based inverse analysis with nonlinearities to recover deformations and loads, enabling precise shape reconstruction of shell structures.

Findings

Accurately recovers target shell shapes

Captures buckling and snapping instabilities

Demonstrates applicability in manufacturing

Abstract

This article presents original work combining a NURBS-based inverse analysis with both kinematic and constitutive nonlinearities to recover the applied loads and deformations of thin shell structures. The inverse formulation is tackled by gradient-based optimization algorithms based on computed and measured displacements at a number of discrete locations. The proposed method allows accurately recovering the target shape of shell structures such that instabilities due to snapping and buckling are captured. The results obtained show good performance and applicability of the proposed algorithms to computer-aided manufacturing of shell structures.

Equations

$$\mathbf{x}=\mathbf{x}({\xi }^{\alpha }),\alpha =1,2,$$

$$A_{\alpha \beta }$$

$$a_{\alpha \beta }$$

$$B_{\alpha \beta }$$

$$b_{\alpha \beta }$$

$$H=\frac{1}{2}{a}^{\alpha \beta }{b}_{\alpha \beta },\kappa =\frac{det[b_{\alpha \beta }]}{det[a_{\alpha \beta }]},{\kappa }_{1/2}=H\pm \sqrt{H^2-\kappa }$$

$$\mathbf{C}={\mathbf{F}}^T\mathbf{F}=a_{\alpha \beta }{\mathbf{A}}^{\alpha }\otimes {\mathbf{A}}^{\beta },\mathbf{B}=\mathbf{F}{\mathbf{F}}^T=A^{\alpha \beta }{\mathbf{a}}_{\alpha }\otimes {\mathbf{a}}_{\beta }$$

$$\mathbf{E}$$

$$\mathbf{K}$$

$${\mathbf{T}}_{;\alpha }^{\alpha }+\mathbf{f}=0.$$

$${\mathbf{T}}^{\alpha }=N^{\alpha \beta }{\mathbf{a}}_{\beta }+S^{\alpha }\mathbf{n}.$$

$$N^{\alpha \beta }$$

$$S^{\alpha }$$

$$\mathbf{u}$$

$$\mathbf{n}$$

$$\mathbf{t}$$

$$m_{\tau }$$

$$\delta \Pi =\delta {\Pi }_{\mathrm{int}}-\delta {\Pi }_{\mathrm{ext}}=0\forall \delta \mathbf{x}\in \mathcal{V},$$

$$\delta {\Pi }_{\mathrm{int}}$$

$$\delta {\Pi }_{\mathrm{ext}}$$

$$\Delta \delta {\Pi }_{\mathrm{ext}}={\int }_{\partial \mathcal{S}}{m}_{\tau }\delta {\mathbf{a}}_{\alpha }\cdot \left({\nu }^{\beta }\mathbf{n}\otimes {\mathbf{a}}^{\alpha }+{\nu }^{\alpha }{\mathbf{a}}^{\beta }\otimes \mathbf{n}\right)\Delta {\mathbf{a}}_{\beta }ds$$

$$\Delta \delta {\Pi }_{\mathrm{int}}={\int }_{{\mathcal{S}}_0}\left(c^{\alpha \beta \gamma \delta }\frac{1}{2}\delta a_{\alpha \beta }\frac{1}{2}\Delta a_{\gamma \delta }+d^{\alpha \beta \gamma \delta }\frac{1}{2}\delta a_{\alpha \beta }\Delta b_{\gamma \delta }+{\tau }^{\alpha \beta }\frac{1}{2}\Delta \delta a_{\alpha \beta }+e^{\alpha \beta \gamma \delta }\frac{1}{2}\delta b_{\alpha \beta }\Delta a_{\gamma \delta }+f^{\alpha \beta \gamma \delta }\delta b_{\alpha \beta }\Delta b_{\gamma \delta }+M_0^{\alpha \beta }\Delta \delta b_{\alpha \beta }\right)dA$$

$$c^{\alpha \beta \gamma \delta }:=2\frac{\partial {\tau }^{\alpha \beta }}{\partial a_{\gamma \delta }},d^{\alpha \beta \gamma \delta }:=\frac{\partial {\tau }^{\alpha \beta }}{\partial b_{\gamma \delta }},$$

$$e^{\alpha \beta \gamma \delta }:=2\frac{\partial M_0^{\alpha \beta }}{\partial a_{\gamma \delta }},f^{\alpha \beta \gamma \delta }:=\frac{\partial M_0^{\alpha \beta }}{\partial b_{\gamma \delta }},$$

$${\tau }^{\alpha \beta }=2\frac{\partial W}{\partial a_{\alpha \beta }},M_0^{\alpha \beta }=\frac{\partial W}{\partial b_{\alpha \beta }}.$$

$$W(\mathbf{E},\mathbf{K})=\frac{1}{2}\mathbf{E}:\mathbb{C}:\mathbf{E}+\frac{1}{2}\mathbf{K}:\mathbb{F}:\mathbf{K},$$

$${\tau }^{\alpha \beta }$$

$$M_0^{\alpha \beta }$$

$$\Lambda :=T\frac{2\tilde{\Lambda }\tilde{\mu }}{\tilde{\Lambda }+2\tilde{\mu }},\mu :=T\tilde{\mu },$$

$$\tilde{\mathbf{X}}({\xi }^{\alpha },\xi )$$

$$\tilde{\mathbf{x}}({\xi }^{\alpha },\xi )$$

$${\tau }^{\alpha \beta }$$

$$M_0^{\alpha \beta }$$

$$\tilde{W}({\tilde{I}}_1,\tilde{J})=\frac{\tilde{\Lambda }}{4}\left({\tilde{J}}^2-1-2ln\tilde{J}\right)+\frac{\tilde{\mu }}{2}\left({\tilde{I}}_1-3-2ln\tilde{J}\right),$$

$${\tilde{\tau }}^{\alpha \beta }=\tilde{\mu }{G}^{\alpha \beta }-\tilde{\mu }\frac{\tilde{\Lambda }+2\tilde{\mu }}{\tilde{\Lambda }{J^{\ast }}^2+2\tilde{\mu }}{g}^{\alpha \beta }.$$

$${\mathrm{N}}_A(\xi ,\eta )=\frac{{\omega }_A{\hat{\mathrm{N}}}_A^e(\xi ,\eta )}{{\sum }_{A=1}^n{\omega }_A{\hat{\mathrm{N}}}_A^e(\xi ,\eta )},$$

Full text

1 Introduction

There exist many examples in the natural world of shape changes of slender structures driven by external sources. For example, biological structures can change shape due to stimuli such as nonuniform heating, local swelling and growth of thin sheets. These geometric modifications formed by such morphological changes can have the benefit of enhancing biological functions[1, 2]. For instance, osmotic swelling leads to snap closing of the leaves of the Venus flytrap, with the resultant configuration being crucial for its nutrition [3]. Since thin structures tend to bend to reduce their stretching energy when subjected to dramatic growth-induced deformations as occurs in growing leaves [4], wrinkling skin [5] and the writhing of tendril-bearing climbers [6], they often morph into nontrivial three dimensional (3D) shapes [7]. Furthermore, large deformations occurring in a growing body result in an interesting shape-altering mechanism, i.e. if specific areas within a thin material are deformed, the entire structure will morph into a new shape. Hence, understanding and controlling the shape changes that result in thin objects subject to external stimuli such as external loads, swelling or temperature-induced thermal strain is crucial for biomimetic engineering. Various studies have been performed to understand the shape change mechanism from swelling or growing of thin, soft bodies [8, 9].

In order to obtain a target shape, the growth-like morphing of a 2D flat sheet into curved, 3D shapes under external stimuli needs to be understood. In particular, the stimuli that is required to produce a target shape needs to be determined and then the resultant configuration of the sheet subject to these stimuli can be reconstructed. Inverse analysis is a widely used approach to determining the external stimuli and the reconstruction of the displacement vector at every material point of the structure [10, 11, 12, 13]. For instance, an inverse method is presented in [14] to determine the active strain needed to deform a bio-hybrid system to a desired curvature that is obtained by appropriate experiments. Inverse techniques are also widely used in aerospace applications. For instance, control-surface morphing usually used in commercial aircraft, employs piezo-electric actuators integrated with shape memory alloys [15]. The method enable identifying a detailed state of structural deformations from a set of experimentally measured displacements. Consequently, other fundamental response quantities such as stress, material and structural failure can be estimated.

Various inverse problems and their applications have been discussed in literature, see [16, 17, 18, 19]. However, very few studies address the reconstruction of three-dimensional deformations of bending structures [20, 21, 22]. Even fewer studies on inverse analysis deal with high-fidelity plate or shell structures as indicated by Tessler et al. [23]. As studied by Jones et al. [24], an inverse approach employing least-squares approach was used to reconstruct deformations of a cantilever plate from a set of strain measurements. Later, Shkarayev [25] developed an inverse approach employing finite element-based parametrization that uses measured surface strains to determine the applied loads, stresses, which then leads to reconstructing structural deformations in an aerospace vehicle. It should be noted that the existing inverse methods are limited to solving small strain elasticity problems, and have not been utilized to account for nonlinear kinematics and large deformations.

For the reconstruction of structural deformations, numerical methods require the solution of inverse problems. Finite element (FE) analysis is commonly used to discretize the spatial domain in combination with the methods of mathematical optimization. Inverse problems may be ill-posed, particularly when considering uniqueness of the solution. In order to obtain a solution, iterative methods of nonlinear optimization are applied. The optimization problem can be solved by using gradient-free or gradient-based methods [26, 27, 28]. In gradient-free methods, only the relationship between an objective and any associated constraint functions and the structural problem are required whereas the gradient-based optimizers request derivatives of the objective and constraint functions. Those can be computed via analytical or semi-analytical sensitivities. Since gradient-free optimization strategies, like genetic algorithms [29] or multilevel coordinate search [30] are often computationally expensive in finding the optimal solution, particularly when the number of parameters to be identified increases [31], gradient-based optimization algorithms are employed in this study.

Kirchhoff-Love shell theory is suitable for thin shells. However, since $C^{1}$-continuity within the computational domain are required, it is rarely used in FE analysis. Meshfree [32, 33, 34] and NURBS-based formulations [35, 36, 37, 38, 39] on the other hand fulfill the $C^{1}$-continuity requirement of Kirchhoff-Love shell theory and allow to avoid rotational degrees of freedom. Due to the lower computational cost and the better opportunity to model complex geometries, we will pursue the NURBS-based approach. NURBS have also been recently used in shape optimization [40, 41] and topology optimization [42, 43]. In this article, we will present a NURBS-based inverse analysis using analytical and semi-analytical sensitivities for shell structures.

To the best of our knowledge, this is the first time an inverse approach is presented for thin shell structures that accounts for both geometric and material nonlinearities. The problem of interest is the determination of the external applied loads and the reconstruction of structural deformations from given data where shape changes due to instabilities like snap-through, snap-back and buckling are allowed.

The article is divided into three major sections. The next section presents a brief description of the theory of thin shells. The governing equation, weak form, and constitutive equations are discussed. The FE formulation is also addressed in this section. The third section describes the inverse analysis based on the gradient-based optimization using analytical and semi-analytical sensitivity approaches. The proposed method is verified and illustrated by various numerical examples in section 4. In particular, elastic instabilities like snap-through, snap-back and buckling will be taken into account.

2 A brief description of rotational-free thin shell theory

2.1 Thin shell kinematics

Consider a general surface $\mathcal{S}$ of the shell in Figure 1 on which a point $\bm{x}$ can be described by the parametric mapping

[TABLE]

where $\xi^{\alpha}$ denotes coordinates in the parameter domain $\mathcal{P}$. The deformation of $\mathcal{S}$ is described by introducing the reference configuration $\mathcal{S}_{0}$. It is characterized by the parametric description $\bm{X}=\bm{X}(\xi^{\alpha})$. From these mappings, the covariant tangent vectors $\bm{A}_{\alpha}=\frac{\partial\bm{X}}{\partial\xi^{\alpha}}$ and $\bm{a}_{\alpha}=\frac{\partial\bm{x}}{\partial\xi^{\alpha}}$ can be determined, respectively. The respective surface normal vectors of the reference and current configurations are then defined by $\bm{N}=\frac{\bm{A}_{1}\times\bm{A}_{2}}{\parallel\bm{A}_{1}\times\bm{A}_{2}\parallel}$ and $\bm{n}=\frac{\bm{a}_{1}\times\bm{a}_{2}}{\parallel\bm{a}_{1}\times\bm{a}_{2}\parallel}$.

Given the covariant tangent vectors, the components of the covariant metric tensor of the midsurface

[TABLE]

and the components of the contravariant metric tensors $[A^{\alpha\beta}]=[A_{\alpha\beta}]^{-1}$ and $[a^{\alpha\beta}]=[a_{\alpha\beta}]^{-1}$ can be evaluated. From the contravariant metric tensors, the contravariant tangent vectors are defined as $\bm{A}^{\alpha}=A^{\alpha\beta}\bm{A}_{\beta}$ and $\bm{a}^{\alpha}=a^{\alpha\beta}\bm{a}_{\beta}$. The covariant components of the curvature tensor are defined as

[TABLE]

The mean, Gaussian and principal curvatures of surface $\mathcal{S}$ are then given by

[TABLE]

The surface deformation gradient tensor $\bm{F}$ describing the mapping from $\mathcal{S}_{0}$ to $\mathcal{S}$ is defined as $\bm{F}=\bm{a}_{\alpha}\otimes\bm{A}^{\alpha}$. Given $\bm{F}$, the right and left Cauchy-Green surface tensors are introduced by:

[TABLE]

The surface Green-Lagrange strain tensor and the curvature tensor can then be expressed as follows [45]

[TABLE]

2.2 Balance laws

Given a surface $\mathcal{S}$ subjected to a prescribed body force $\bm{f}$, the quasi-static equilibrium of the thin shell is described by

[TABLE]

Here, $\bm{T}^{\alpha}_{;\beta}$ denotes the covariant derivative of $\bm{T}^{\alpha}$ defined by $\bm{T}^{\alpha}_{;\beta}=\bm{T}^{\alpha}_{,\beta}+\Gamma^{\alpha}_{\beta\gamma}\bm{T}^{\gamma}$, where $\Gamma^{\alpha}_{\beta\gamma}$ are the Christoffel symbols defined by $\Gamma^{\gamma}_{\alpha\beta}=\bm{a}^{\gamma}\cdot\bm{a}_{\alpha,\beta}$. $\bm{T}^{\alpha}$ denotes the traction defined on the surface normal to $\bm{a}^{\alpha}$, given by

[TABLE]

Here the respective in plane and shear stress components of the distributed sectional force are defined by

[TABLE]

where $\sigma^{\alpha\beta}$ and $M^{\alpha\beta}$ are the respective membrane stress and bending moment components [46]. The boundary conditions at the boundary $\partial\mathcal{S}$ are expressed as

[TABLE]

with $\bar{\bm{u}}$, $\bar{\bm{n}}$, $\bar{\bm{t}}=\bar{t}^{\alpha}\bm{a}_{\alpha}$, and $\bar{m_{\tau}}$ being the respective prescribed displacement, rotation, boundary traction and bending moment. $m_{\tau}$ and $m_{\nu}$ are the bending moment components parallel and perpendicular to boundary $\partial\mathcal{S}$, respectively.

2.3 Weak form

It can be shown that the weak form of Equation (7) is given by the statement [47]: Find $\bm{x}\in\mathcal{V}$ $\forall\delta\bm{x}\in\mathcal{V}_{0}$ such that

[TABLE]

with the internal and external virtual work given by

[TABLE]

where the last term in Equation (13) refers to the virtual work of the point loads $m_{\nu}\bm{n}$ applied at corners on the boundary $\partial_{m}\mathcal{S}$ (in case $m_{\nu}\neq 0$). It can be shown that after linearization the external and the internal virtual work are characterized by the increments

[TABLE]

[TABLE]

in which $c^{\alpha\beta\gamma\delta}$, $d^{\alpha\beta\gamma\delta}$, $e^{\alpha\beta\gamma\delta}$, and $f^{\alpha\beta\gamma\delta}$ are defined by

[TABLE]

and where $\tau^{\alpha\beta}=J\sigma^{\alpha\beta}$, $M^{\alpha\beta}_{0}=JM^{\alpha\beta}$, and $da=J~{}dA$. The Newton-Raphson and the arc length methods [48] are employed to solve the linearized system of equations.

2.4 Constitutive equations

A Neo-Hookean material model is used to describe nonlinear stress-strain relationship of the materials under large deformations. This model is non-dissipative. For hyperelastic material models, constitutive equations for stretching and bending are derived from a strain energy function $W$, i.e.

[TABLE]

2.5 Koiter shell material model

According to the Koiter material model [49, 50], the surface strain energy for initially curved shells is defined as

[TABLE]

in which $\mathbb{C}=\Lambda\bm{I}\odot\bm{I}+2\mu\left(\bm{I}\otimes\bm{I}\right)$ and $\mathbb{F}=\frac{T^{2}}{12}\mathbb{C}$. Consequently, the Kirchhoff stress and moment in Equation (17) are

[TABLE]

where $\Lambda$ and $\mu$ can be determined in different ways. As presented in [44], we can integrate analytically over the thickness of the 3D Saint Vernant-Kirchhoff material model to obtain

[TABLE]

with $\tilde{\Lambda}$ and $\tilde{\mu}$ being the classical 3D Lamé constants in linear elasticity.

2.6 Shell constitutive model derived from 3D material model

In this section, a Kirchhoff-Love shell constitutive model is derived from a three-dimensional continuum mechanics. The 3D model is projected onto the surface $\mathcal{S}$ to obtain the 2D surface shell model as shown in Figure 1. The actual metric at an arbitrary material point $P$ in the shell has to be mapped onto the midsurface in the respective reference and deformed configuration [51] as

[TABLE]

where $\xi\in[-T/2~{}T/2]$ is the coordinate along the shell thickness and $\bm{d}:=\lambda_{3}\bm{n}$ with $\lambda_{3}$ being the stretch in the normal direction for the Kirchhoff-Love shell model. The components of the projected Kirchhoff stress and moment in Equation (17) can be expressed as

[TABLE]

where $\tilde{\tau}^{\alpha\beta}$ denote the 3D Kirchhoff stress components. The detailed projection procedure can be found in [44].

2.7 Compressible Neo-Hookean materials

The 3D compressible Neo-Hookean strain energy per reference configuration is provided as follows:

[TABLE]

with $\tilde{I}_{1}$ and $\tilde{J}$ being invariants of the 3D Cauchy-Green tensor $\tilde{\bm{C}}$. The 3D Kirchhoff stress is then given by [44] as

[TABLE]

where the respective covariant and contravariant metric tensors are defined as $g_{\alpha\beta}:=\bm{g}_{\alpha}\cdot\bm{g}_{\beta}$ and $g^{\alpha\beta}:=\bm{g}^{\alpha}\cdot\bm{g}^{\beta}$, where $\bm{g}_{\bullet}$ and $\bm{g}^{\bullet}$ denote the covariant and contravariant tangent vectors at P, respectively. Further $J^{*}=\sqrt{\frac{g}{G}}$ with $G:=\mathrm{det}[G_{\alpha\beta}]$ and $g:=\mathrm{det}[g_{\alpha\beta}]$.

By substituting the stress in Equation (24), the stress and moment shown in Equation (22) can be evaluated. In this article, we will consider two constitutive material models: (1) the Koiter material model based on Equation (19) and (2) the projected material model based on Equations (22) and (24).

2.8 FE approximation

Because the Kirchoff-Love shell formulations involve second derivatives, $C^{1}$-continuity of the surface is required. Hence, NURBS-based shape functions proposed by [52] are used to discretize the surface $\mathcal{S}$ and to solve Equation (11). The NURBS basis functions are defined by

[TABLE]

in which $\{\hat{\mathrm{N}}^{e}_{A}\}^{n}_{A=1}$ is the B-spline basis function written in terms of Bernstein polynomials as

[TABLE]

with $\hat{\mathrm{N}}^{e}_{A}$ being entries of matrix $\hat{\mathbf{N}}^{e}$. The geometry of the respective reference and deformed surfaces are approximated by interpolation of the control point coordinates $\bm{X}_{e}$ and $\bm{x}_{e}$ using the FE approach as follows:

[TABLE]

where $\mathbf{N}(\xi,\eta):=[N_{1}\bm{1},~{}N_{2}\bm{1},~{}...,~{}N_{n}\bm{1}]$ contains the nodal shape functions of element $\Omega^{e}$ shown in Equation (25). The covariant tangent vectors of the surface are then determined by

[TABLE]

Variation of $\bm{x}$ and $\bm{a}_{\alpha}$ are expressed as

[TABLE]

Consequently, all kinematical quantities presented in section 2.1 and their variation can be determined, see [47, 44] for details.

2.9 Discretized weak form

The weak form in Equation (11) is discretized using the above interpolations. This leads to the approximation

[TABLE]

with $n_{el}$ being the number of elements. The internal and external virtual work contributions of element $e$ are given by

[TABLE]

in which the internal FE force vectors due to the membrane stress $\tau^{\alpha\beta}$ and the bending moment $M^{\alpha\beta}_{0}$ are defined by [44]

[TABLE]

and the external force FE vectors due to a constant body force $\bm{f}_{0}$, external pressure perpendicular to the surface $\mathcal{S}$, boundary traction $\bm{t}$ and boundary moment $m_{\tau}$ are defined by [44]

[TABLE]

The linearization of $\delta\Pi^{e}_{\mathrm{int}}$ and $\delta\Pi^{e}_{\mathrm{ext}}$ is shown in A.

2.10 Lagrange multiplier method for rotational constraints

This section presents an approach used to enforce $G^{1}$-continuity between patches. The rotational constrains is introduced by adding the constraint potential, see [44]

[TABLE]

to the shell formulation, with $p$ being the Lagrange multiplier and

[TABLE]

where $\cos\theta_{0}:=\bm{N}\cdot\bar{\bm{N}}$ and $\cos\theta:=\bm{n}\cdot\bar{\bm{n}}$. Here $\bm{N}$, $\bar{\bm{N}}$ are the two surface normals of neighboring patches in the reference configuration, and $\bm{n}$, $\bar{\bm{n}}$ are the corresponding normals in the deformed configuration. Details of the variation, linearization and FE discretization of $\Pi_{n}$ can be found in [44]. The particular patch boundary denoted as $\mathcal{L}_{0}$ in Equation (34) refers to the reference configuration.

3 Inverse analysis

This section presents an inverse analysis procedure for the isogeometric rotation-free thin shell considering hyperelastic material behavior. Nonlinear kinematics and material response are considered. The boundary conditions, i.e. prescribed boundary displacements, tractions and moments, as well as surface loads like out-of-plane pressure, that are applied to produce a desired shape of thin shell structures, can be identified. The considered technique is based on gradient-based strategies, where the required derivatives are determined either analytically, by the finite difference method or by the adjoint approach. Inverse problems are often solved iteratively in form of the optimization problem [53]

[TABLE]

in which $s$ are the design variables; $u$ are the state variables describing the structural response, e.g. the displacements; $\mathrm{J}$ is the objective function; $g_{j}$ are inequality constraints imposed on these variables to restrict quantities describing the shell response; $s_{l}$ and $s_{u}$ are the respective side constraints which restrict the design variables. In this study, the Dirichlet boundary conditions - the parameters $s$ - can be determined based on the measured displacements of the surface. The forward operator which maps those parameters to measurements is shown by [31] as follows:

[TABLE]

where $X$ refers to the parameter space, i.e. the range of Neumann or Dirichlet boundary conditions (external loads, prescribed displacements etc.) and $Y$ refers to the measurement space. Given the measurements with noise $u^{meas}$, the inverse problem corresponds to determining $s$ from the following equation

[TABLE]

which is rewritten in terms of a least-squares functional form

[TABLE]

The inverse analysis will be employed to identify unknown applied loads and reconstruct nonlinear deformations of the shell structure where the shape changes due to instabilities (i.e. snapping or buckling) are allowed. The above objective function in Equation (41) is expressed in terms of the FE discretization as

[TABLE]

where $\bm{\mathrm{u}}$ is the vector of nodal displacements on the surface; $\bm{\mathrm{u}}^{meas}$ is the vector of the displacements measured on the target shape. The design variables are chosen as the applied loads. The state variables are the nodal displacements. The gradient based method of moving asymptotes (MMA) [54] is employed together with analytically derived sensitivities. Thus, the differentiation of the objective function $\mathrm{J}$ and constraint $g$ with respect to (w.r.t.) design variables $\bm{s}$ is required. Using the chain rule of differentiation, the total derivative of the objective in Equation (42) w.r.t. the design variable $s_{i}$ is expressed as

[TABLE]

The partial derivative w.r.t. the design variable $\frac{\partial\mathrm{J}}{\partial\bm{s}}$ is zero according to Equation (42). The equilibrium equations for nonlinear FE analysis is given by

[TABLE]

where $\bm{f}_{\mathrm{int}}$ and $\bm{f}_{\mathrm{ext}}$ denote the internal force and the external applied load, respectively. Differentiating Equation (44) w.r.t. the design variables $\bm{s}$ gives

[TABLE]

Rewriting this equation gives the nonlinear state derivative w.r.t. the design variables $\frac{\partial\bm{u}}{\partial\bm{s}}$ as

[TABLE]

where $\bm{K}_{T}$ is the tangent stiffness matrix; $\bm{f}^{*}_{\mathrm{nln}}$ is known as the nonlinear pseudo load vector [53]. In optimization algorithms, the derivatives of the displacement $\frac{\partial\bm{u}}{\partial s_{i}}$ are not computed explicitly. Instead, they are obtained by using the adjoint method that defines an adjoint vector $\bm{\lambda}$ [55] which is the solution of the equation

[TABLE]

where $\bm{z}$, with $z_{i}=\frac{\partial J}{\partial u_{i}}$, is the vector of derivatives of the objective function w.r.t. the displacement. The solution $\bm{\lambda}$ of Equation (47) is a FE solution for displacements and is obtained by solving the governing equation (7) with Neumann boundary condition defined by

[TABLE]

and the Dirichlet boundary conditions are the same as those of the actual problem shown in Equations (10.1). Then, Equation (43) is written as

[TABLE]

As observed, the derivatives of $\bm{f}_{\mathrm{int}}$ and $\bm{f}_{\mathrm{ext}}$ w.r.t. the design variables $s_{i}$ in Equation (49) are required in the adjoint method. While it is possible to use numerical sensitivities, it is more efficient to use analytical sensitivities. The analytical sensitivities will be presented in the following section 3.1 whilst the semi-analytical sensitivities will be described in section 3.2.

3.1 Analytical sensitivities

In this section, the analytical sensitivities of the quantities in Equation (49) w.r.t. the respective prescribed displacement ($\bm{u}$), traction ($\bm{t}$), point load ($\bm{f}_{0}$), moment ($m_{\tau}$) and lateral pressure ($q$) are presented.

3.1.1 Prescribed displacement

Assume a vector of displacements is prescribed at $\bm{X}_{A}$ along $\bm{e}_{i}$, the design variable is defined as $\bm{s}=\bm{u}_{A}$ such that $\bm{u}_{A}=u~{}\delta(\bm{X}_{A})~{}\bm{e}_{i}$. The nonlinear pseudo load vector $\bm{f}^{*}_{\mathrm{nln}}$ is given by

[TABLE]

where $\bm{K}_{T,A}$ refers to all columns of the tangent stiffness matrix $\bm{K}_{T}$ defined in Equation (46) that correspond to the degrees of freedom (dofs) in the design vector $\bm{s}$. Note that here $\delta$ is the Dirac-Delta distribution. Implementation should be done at the global level (i.e. after elemental assembly). By doing so, the sensitivity for Equation (49) is obtained.

3.1.2 Traction

For an unknown edge traction along direction $\bm{e}_{i}$ the design variable $s$ is

[TABLE]

The nonlinear pseudo load vector $\bm{f}^{*}_{\mathrm{nln}}$ is obtained by

[TABLE]

where $\bm{f}^{e}_{\mathrm{ext}}$ is given by Equation (33.3). Since now $\bm{t}=s\bm{e}_{i}$, we find

[TABLE]

More general, if $\bm{s}=\bm{t}$, the pseudo load vector is

[TABLE]

And if $s=\bm{t}_{0}$, where $\bm{t}_{0}=\bm{t}\frac{ds}{dS}$, we find

[TABLE]

3.1.3 Point load

Considering a shell structure under a point load $\bm{F}$, at $\bm{X}_{A}$ along $\bm{e}_{i}$, i.e. $\bm{F}=F~{}\delta(\bm{X}_{A})~{}\bm{e}_{i}$, the design variable $s$ is equal to $F$. The nonlinear pseudo load vector $\bm{f}^{*}_{\mathrm{nln}}$ is obtained by

[TABLE]

where $\bm{f}^{e}_{\mathrm{ext}}$ is given by Equation (33.1). The only non-zero entry within Equation (33.1) is $\bm{f}_{\mathrm{ext}A}=\mathrm{N}_{A}~{}F~{}\bm{e}_{i}$ where $A$ denotes the FE node at $\bm{X}_{A}$. Hence, Equation (56) becomes

[TABLE]

More general, in case $\bm{s}=\bm{f}_{0}$, $\bm{f}^{*}_{\mathrm{nln}}$ is given by

[TABLE]

3.1.4 Uniaxial load and bending moment

Let consider a shell structure subjected to uniaxial traction and bending moment as shown in example 4.3. The two design variables are $s_{1}=\bm{t}\cdot\bm{e}_{1}$ and $s_{2}=m_{\tau}$. For the uniaxial traction, the derivative $\bm{f}^{*}_{\mathrm{nln1}}=\frac{\partial\bm{f}^{e}_{\mathrm{ext}}}{\partial s_{1}}$ is given in Equation (53). Taking the derivative of $\bm{f}^{e}_{\mathrm{ext}}$ w.r.t. $s_{2}$, we obtain

[TABLE]

where $m_{\tau}$ denotes the distributed moment per current length. Assuming that $m_{\tau_{0}}$ = distributed moment per reference length is the design variable, then $\bm{f}^{e}_{\mathrm{ext}}=\int_{\partial_{m}\Omega^{e}_{0}}\mathbf{N}^{T}_{,\alpha}\nu^{\alpha}m_{\tau_{0}}\bm{n}~{}dS$ and

[TABLE]

3.1.5 Uniaxial load and lateral pressure

Let us consider a shell structure that is subjected to uniaxial traction and lateral pressure as shown in example 4.4. In this case, the two design variables will be $s_{1}=\bm{t}\cdot\bm{e}_{1}$ and $s_{2}=q$. The derivatives $\bm{f}^{*}_{\mathrm{nln1}}=\frac{\partial\bm{f}^{e}_{\mathrm{ext}}}{\partial s_{1}}$ is given in Equation (53) and $\bm{f}^{*}_{\mathrm{nln2}}=\frac{\partial\bm{f}^{e}_{\mathrm{ext}}}{\partial s_{2}}$ is given by

[TABLE]

3.2 Semi-analytical sensitivities

The semi-analytical sensitivities of the quantities in Equation (49) are approximated by the finite difference method as follows:

[TABLE]

where $f_{\bullet}(\bm{u},\bm{s}+\Delta\bm{s}_{i})$ refers to perturbed quantities in direction of variable $s_{i}$. The small step size $\Delta s_{i}=0.001\times s_{i}$ is selected to reduce the truncation error to acceptable value [55]. The details of the selection of the step size can be found in [56, 57].

4 Numerical examples

Since experimental measurements are not available for this study, the proposed approach is validated using a two-step procedure. First, a FE model on the basis of the isogeometric approach [58] is developed for a given shell structure under the applied loads. Both nonlinear kinematics and nonlinear material constitutive laws are accounted for. Surface displacements $\bm{\mathrm{u}}^{iga}$ are then computed at discrete locations. Random noise representing small disturbances is added to the displacement vector $\bm{\mathrm{u}}^{iga}$, giving $\bm{\mathrm{u}}^{meas}:=\bm{\mathrm{u}}^{iga}(1+0.01\gamma)$, with $\gamma$ being a random number in the interval $[-1,1]$. Since this disturbed data could originate from an experiment we call it experiment-like. Secondly, using this experiment-like data, we apply gradient-based optimization methods on Equation (42). The approximation of the applied loads are inversely solved and the structural deformations are reconstructed. Solution of the inverse analysis is achieved when insignificant changes occur in the configuration between iterations. Particularly, the algorithm is completed where the convergence criterion, $\epsilon_{obj.func.}=|\frac{J^{I}-J^{I-1}}{J^{0}}|\leqslant 10^{-3}$, is satisfied with $I$ being the iteration number and $J^{0}$ is the objective value computed at the initial input.

Four numerical examples on cantilevered and hinged cylindrical shells governed by geometrically nonlinear mechanics are presented. Hyperelastic constitutive laws which account for stretching and bending are considered. The first example is a cantilever shell subjected to prescribed end displacement (also end shear traction), the second is the hinged cylindrical shell subjected to a prescribed center displacement. In the two last examples, imperfections, which are essential for instability shape changes are considered. The inverse analysis of nonlinear instabilities of shell structures under compression using both the Koiter material model and the projected shell material model using the compressible Neo-Hookean formulation are examined. Note that hereafter we use the Koiter model to refer to the former and the projected model for the latter.

4.1 Cantilever subjected to prescribed end displacement and end shear traction

In this example, the prescribed displacement $w_{tip}$ or the shear traction $t$ applied to the free end of a cantilever shell is identified using experiment-like data, see Figures 2(a) and 2(b). The experiment-like data is generated by considering the large deflection of a cantilever shell, with dimensions $L\times W\times T=1\times 10\times 0.1$ mm, under the prescribed displacement and the end shear traction as shown in Figure 2. The shell is discretized by $4\times 40$ NURBS elements. The projected model using the compressible Neo-Hookean formulation, according to sections 2.6 and 2.7 is used with $E=1.2\times 10^{6}~{}N/mm^{2}$ and $\nu=0$. The force-displacement responses obtained from the forward problem are compared to that reported by Sze et al. [59]. Convergence of the energy norm error determined by Equation (63) is shown in Figure 3 to assess the FE results. The error $e_{E}$ is calculated as

[TABLE]

where $E^{num}$ and $E^{ref}$ are the numerical and the reference strain energies, respectively. It should be noted that since the analytical solution is not available, the reference strain energy obtained from a fine mesh of $20\times 200$ NURBS-based elements is used. The experiment-like data, i.e. the nodal displacements, is obtained by using NURBS-based finite elements. The target configuration corresponds to point A on the solid blue curve in Figure 2(c) where $w^{meas}_{tip}=4.0$ or $t^{meas}=1.5~{}N/mm$.

Inverse analysis based on the target configuration is then carried out corresponding to Neumann (traction $t$) and Dirichlet (prescribed displacement $w_{tip}$) boundary conditions (BCs). The element derivatives in Equation(49) are obtained by using analytical and semi-analytical sensitivities. Detailed implementation for the analytical sensitivities when the cantilever shell subjected to the displacement $w_{tip}$ or the end shear traction $t$ can be found in sections 3.1.1 and 3.1.2, respectively. Accordingly, the convergence of the objective functions and $L^{2}$ error norms in spaces $Y$ and $X$, see Equation (38), w.r.t. iterations is shown in Figures 4 and 5. The inverse solutions $t^{inverse}=1.4927~{}N/mm$ (for Neumann BCs) and $w^{inverse}_{tip}=4.0077$ (for Dirichlet BCs) are obtained respectively. It can be seen that they are in good agreement with $t^{meas}=1.5~{}N/mm$ $w^{inverse}_{tip}=4$.

In addition, inverse approach using the semi-analytical sensitivity is performed to reconstruct the target configuration of the cantilever shell. Figure 6 shows the excellent accuracy of the proposed inverse method. The convergence of the objective function and $L^{2}$ error norm over iterations is illustrated in Figures 6(a+b), respectively. Very good agreement of the convergence of the objective function w.r.t. the number of iterations is obtained between the analytical and semi-analytical sensitivities, see Figure 6(c). The deformed configuration under the inverse solution $w^{inverse}_{tip}=3.9696$ shows excellent agreement with the target configuration in Figure 2(d). As observed, the maximal error ($err_{\bm{w}}=\frac{4-3.9696}{4}=0.76\%$) between two configurations occurs at the nodes along the free end.

In order to estimate the effect of the noise on accuracy of the algorithm, systematic noise ($\gamma$ is taken as deterministic numbers) is examined. As shown in Figure 7(a), the accuracy of the algorithm decreases with an increase in systematic noise $\gamma$. The convexity of the objective function w.r.t. the end displacement $w_{tip}$ can be seen in Figure 7(b).

Generally, inverse problems may be ill-posed, i.e., that one of the following conditions is violated: Existence, uniqueness and stable dependency on the data. As we do not aim to solve the inverse problem directly (Equation (40)) but make use of a sum of least-squares approach (Equation (41)) the existence is guaranteed, at least in the least-squares sense. Proving uniqueness analytically is challenging and most probable not possible for the complex nonlinear forward problem. So consequently, we followed a numerical approach where values of the objective are plotted over different values for the excitations, see Figure 7(b). As a clear minimal point is indicated, we can assume uniqueness at least for the simplified problems. For problems with more degrees of freedom a generalization is unfortunately impossible.

To solve now the inverse problem, iterative methods are applied. Regarding error norms in $X$ and $Y$ space we see good convergence and no strong tendency for divergence, even for higher error amounts in the data (here the prescribed shape), see Figures 4, 5 and 7(a).

4.2 Hinged cylinderical shell subjected to prescribed center displacement and central point force

The second benchmark example consists of a cylindrical shell of dimensions $R=2540~{}mm$, $L=B=504~{}mm$ and thickness $T=12.7~{}mm$. The longitudinal edges of the shell are pinned, its curved edges are free and the shell is subjected to a prescribed center displacement $w_{cen}$. Due to the symmetry, a quarter of the shell is modeled using $8\times 10$ quadratic NURBS elements as shown in Figure 8(a). The symmetry boundary conditions are enforced by the Lagrange multiplier method presented in section 2.10.

The forward problem is first solved based on the isogeometric analysis approach. The material is described by the Koiter model and the projected model with compressible Neo-Hookean formulation using the Young’s modulus $E=3105~{}N/mm^{2}$ and the Poisson’s ratio $\nu=0.3$. The force-displacement response is compared in Figure 8(c) to the one reported by Sze et al. [59], which considered the Saint Venant-Kirchhoff formulation with the same material properties. Good agreement between the response obtained by using the Koiter model with the reference result demonstrates the accuracy of the model. The deformed configuration corresponding to point A on the solid blue curve in Figure 8(c) is then taken as the target shape.

Next, the determination of the prescribed center displacement $w_{cen}$ using the target shape shown in Figure 8(b) is performed in the framework of compressible Neo-Hookean hyperelasticity to take nonlinear bending and stretching behaviors into account. Analytical and semi-analytical sensitivities are utilized to compute the derivatives in Equation(49). Computation of the analytical sensitivities for prescribed displacement is provided in section 3.1.1. As shown in Figure 9, the objective function and the $L^{2}$ norm converge after $20$ iterations. The semi-analytical sensitivities are also used to provide the derivatives in Equation (62) based on how the inverse problem is carried out. The accuracy of the inverse method is highlighted in Figures 10(a+b) where convergence is gained after $20$ iterations. The solution of the inverse problem over iterations is shown in Figures 10(c+d+e+f). The reconstructed deformation after $20$ iterations shows good agreement with the target shape: $w^{meas}_{cen}=-24.98$ compares well to the maximal displacement $w^{inverse}_{cen}=-24.90$ at the central point. It should be noted that the reconstruction of the target shape is performed in the presence of a snap-through instability.

In order to take the uniqueness of the solution into consideration, the hinged shell problem with smaller thickness ($T=6.35~{}mm$) and snap-back behavior is also simulated. For this model, the prescribed center displacement $w_{cen}$ in Figure 8 is replaced by a central point load $F=700~{}N$. The forward problem is implemented using the arc length method. The force-displacement responses for the Koiter model and the projected model with compressible Neo-Hookean material are shown in Figure 11. The target configuration corresponding to point A on the blue force-displacement curve is obtained and used for the identification of the applied external load.

Inverse analysis using the analytical sensitivity presented in Equation (57) is then performed based on the experiment-like data obtained from the forward problem. The projected model using the compressible Neo-Hookean material is employed. The solution to the inverse problem is obtained after $20$ iterations as depicted in Figure 12. The reconstructed configuration shown in Figure 12(d) under the inverse solution $F^{inverse}=678.3~{}N$ shows good agreement with the target configuration. It should be noted that this configuration is reconstructed under a snap-back instability.

4.3 A flat strip subjected to compressive displacement and bending moment

In reality, a perfect structure with perfect loading does not exist. Inevitable imperfections may be related to factors such as initial curvature, load eccentricity or small disturbing loads. Small imperfections do not matter in linear analysis, however in nonlinear problems, they are quite important. Hence, instabilities of shell structures that are sensitive to geometrical and loading imperfections must be considered. The shape changes occurring due to nonlinear buckling will be examined in the following examples.

We consider a pinned-end flat strip that it is perfectly straight before the displacement is prescribed. The strip is prescribed by compressive displacement and end distributed moment as described in Figure 13. The moment $m$ is applied first and we then apply the compressive displacement $v_{right}$. A FE solution for of $8\times 24$ quadratic NURBS elements is first determined for the above structure with both the Koiter and projected shell models using the compressible Neo-Hookean formulation. The 3D Young’s modulus $E=3\times 10^{6}~{}N/mm^{2}$ and Poisson’s ratio $\nu=0.3$ are adopted. The FE force vector associated with the distributed moment is given by Equation (33) and the corresponding tangent matrix is given by Equation (69.1). The effect of the distributed moment on the force-displacement response for both models is indicated in Figures 13(b+c), respectively. The unstably buckled shape under $v^{meas}_{Koit,right}=-0.5$ and $m^{meas}_{Koit}=12\times 10^{3}~{}N.mm/mm$ is illustrated in Figure 14 for the Koiter model. Furthermore, Figure 16 shows the unstably buckled shape under $v^{meas}_{proj,right}=-0.4$ and $m^{meas}_{proj}=10\times 10^{3}~{}N.mm/mm$ for the projected model. The nodal displacements $\bm{u}^{meas}$ corresponding to point A on the curves in Figures 14(b) and 16(b) are then computed which are used as the experiment-like data.

In the second step, an inverse approach based on the nodal displacement of target configuration (with additional noise) is employed to identify the applied displacement $v_{right}$ and moment $m$ for the two shell models: (1) the Koiter model and (2) the projected model with the compressible Neo-Hookean formulation. For the shell structure modeled by the Koiter model as shown in Figure 15(a), convergence is gained after $40$ iterations and the optimal shape in Figure 15(e) recovered after $40$ iterations, corresponding to the inverse solution $v^{inverse}_{Koit,right}=-0.5055$ and $m^{inverse}_{Koit}=9.06\times 10^{3}~{}N.mm/mm$, which is in excellent agreement with the target shape in Figure 14(a).

Similar observation can be made for the results obtained from the projected model. The convergence of the objective function over iterations and the reconstructed configurations with iterations are depicted after $40$ iterations Figures 17. The recovered shape in Figure 17(e), corresponding to the inverse solution $v^{inverse}_{proj,right}=-0.3989$ and $m^{inverse}_{proj}=11.13\times 10^{3}~{}N.mm/mm$, is almost identical with the measured shape in Figure 16(a). In summary, the proposed inverse method is able to accurately determine the prescribed displacement $v_{right}$ and moment $m$ and the specified final shape can be reconstructed where the unstable deformation due to buckling are captured. Note that nodal displacements in the $z$ direction for both models are restricted to be larger than [math] to ensure the convex shape of the measured shells in Figures 14(a) and 16(a), respectively.

4.4 A flat strip subjected to compressive displacement and disturbing pressure

The pinned-end flat strip, with dimensions $L\times W\times T=5\times 1\times 0.2~{}mm$ is considered again. We assume its initial shape is perfectly straight but subjected to axially prescribed displacement $v_{right}$ and lateral pressure $q$ as depicted in Figure 18. First, NURBS-based finite elements are used to solve the forward problem. This example was analyzed using $8\times 40$ NURBS elements with the material constants: $E=3\times 10^{6}~{}N/mm^{2}$ and $\nu=0.3$. The FE force vector due to lateral pressure is given by Equation (33.2) and the corresponding stiffness matrix $k^{e}_{\mathrm{ext}p}$ is determined by (69.2). The Koiter model and the projected model using the compressible Neo-Hookean formulation are considered. Figures 18(b+c) respectively show the influence of the disturbing pressure on the force-displacement response for both models. The measured shapes and force-displacement response, corresponding to $v^{meas}_{Koit,right}=-0.2$, $q^{meas}_{Koit}=30~{}N/mm^{2}$, are shown in Figure 19 for the Koiter model and in Figure 21 for the projected model with $v^{meas}_{proj,right}=-0.2$, $q^{meas}_{proj}=30~{}N/mm^{2}$, respectively.

Based on the nodal displacements $\bm{u}^{meas}$ corresponding to point A on the curves obtained from the forward problem, we conduct an inverse analysis for the strip modeled by (1) the Koiter model and (2) the projected model with the compressible Neo-Hookean formulation. The convex shape of the reconstructed deformations is ensured by restricting nodal displacements in the $z$ direction to be larger than [math]. For the former model, convergence is reached after $40$ iterations, see Figure 20(a) and the resultant shape under the identified $v^{inverse}_{Koit,right}=-0.1992$, $q^{inverse}_{Koit}=35.18~{}N/mm^{2}$ is accurately reconstructed after $40$ iterations in Figure 20(e) compared to the one in Figure 19(a). Similarly, the objective function versus iterations and the reconstructed configurations are illustrated in Figure 22 for the latter. By comparing Figure 22(e) to Figure 21(a), it can be seen that the specified target shape under the identified $v^{inverse}_{proj}=-0.1999$, $q^{inverse}_{proj}=35.74~{}N/mm^{2}$ can be recovered nearly exactly. Overall, the inverse formulation correctly recovers the prescribed displacement $v_{right}$, the lateral pressure $q$ and the measured deformation due to buckling.

5 Conclusions

In this article, we have developed a new method to combine inverse analysis with material and geometric nonlinearities to identify the unknown applied loads from given displacements and reconstruct the deformations of thin shell structures. NURBS-based finite elements, which offer geometrically exact discretization, are used to ensure $C^{1}$-continuity of the 3D surface. Both a Koiter material model and a projected shell material model with compressible Neo-Hookean formulation are considered. Gradient-based optimization algorithms with analytical and semi-analytical sensitivities were utilized to tackle the inverse problems for given experiment-like displacements at specified locations. Various numerical simulations are performed to capture either the stable or unstable (i.e. snap-through, snap-back and buckling) shape changes. The solution of the inverse analysis demonstrates that the proposed method is able to recover target shapes with high accuracy. We believe this research will open a new path for computer-aided manufacturing of shell structures.

6 Acknowledgments

We gratefully acknowledge the support by ERC COMBAT project (project number 615132). We would like to thank Prof. Krister Svanberg from Royal Institute of Technology for providing the MMA code.

Appendix A FE tangent matrices

According to [44], the linearization of $\delta\Pi^{e}_{\mathrm{int}}$ and $\delta\Pi^{e}_{\mathrm{ext}}$ is expressed as

[TABLE]

where the material stiffness matrices are defined by

[TABLE]

where $\mathbf{N}_{,\alpha}(\bm{\xi}):=[N_{1,\alpha}\bm{1},N_{2,\alpha}\bm{1},...,N_{n,\alpha}\bm{1}]$ and

[TABLE]

with $\mathbf{N}_{,\alpha\beta}(\bm{\xi}):=[N_{1,\alpha\beta}\bm{1},N_{2,\alpha\beta}\bm{1},...,N_{n,\alpha\beta}\bm{1}]$. The geometric stiffness matrices are shown as

[TABLE]

in which

[TABLE]

The external moment and surface pressure components of the tangent matrix are given by [44, 60]

[TABLE]

where $\xi$ refers the convective coordinate of the curve $\partial_{m}\Omega^{e}$.