Full Text:   <2784>

Summary:  <1650>

CLC number: TU383

On-line Access: 2014-04-03

Revision Accepted: 2014-01-22

Crosschecked: 2014-03-17

Cited: 5

Clicked: 5809

Citations:  Bibtex RefMan EndNote GB/T7714

 Journal of Zhejiang University SCIENCE A 2014 Vol.15 No.4 P.255-271 http://doi.org/10.1631/jzus.A1300245

An efficient numerical shape analysis for light weight membrane structures*

 Author(s):  Chao Yang, Yan-bin Shen, Yao-zhi Luo Affiliation(s):  . Space Structures Research Center, Zhejiang University, Hangzhou 310058, China Corresponding email(s):   benjamin127@163.com Key Words:  Tension membranes, Finite particle method (FPM), Shape analysis, Explicit time integration, Initial equilibrium shape Share this article to： More <<< Previous Article|Next Article >>>

Chao Yang, Yan-bin Shen, Yao-zhi Luo. An efficient numerical shape analysis for light weight membrane structures[J]. Journal of Zhejiang University Science A, 2014, 15(4): 255-271.

@article{title="An efficient numerical shape analysis for light weight membrane structures",
author="Chao Yang, Yan-bin Shen, Yao-zhi Luo",
journal="Journal of Zhejiang University Science A",
volume="15",
number="4",
pages="255-271",
year="2014",
publisher="Zhejiang University Press & Springer",
doi="10.1631/jzus.A1300245"
}

%0 Journal Article
%T An efficient numerical shape analysis for light weight membrane structures
%A Chao Yang
%A Yan-bin Shen
%A Yao-zhi Luo
%J Journal of Zhejiang University SCIENCE A
%V 15
%N 4
%P 255-271
%@ 1673-565X
%D 2014
%I Zhejiang University Press & Springer
%DOI 10.1631/jzus.A1300245

TY - JOUR
T1 - An efficient numerical shape analysis for light weight membrane structures
A1 - Chao Yang
A1 - Yan-bin Shen
A1 - Yao-zhi Luo
J0 - Journal of Zhejiang University Science A
VL - 15
IS - 4
SP - 255
EP - 271
%@ 1673-565X
Y1 - 2014
PB - Zhejiang University Press & Springer
ER -
DOI - 10.1631/jzus.A1300245

Abstract:
The determination of initial equilibrium shapes is a common problem in research work and engineering applications related to membrane structures. Using a general structural analysis framework of the finite particle method (FPM), this paper presents the first application of the FPM and a recently-developed membrane model to the shape analysis of light weight membranes. The FPM is rooted in vector mechanics and physical viewpoints. It discretizes the analyzed domain into a group of particles linked by elements, and the motion of the free particles is directly described by Newton’s second law while the constrained ones follow the prescribed paths. An efficient physical modeling procedure of handling geometric nonlinearity has been developed to evaluate the particle interaction forces. To achieve the equilibrium shape as fast as possible, an integral-form, explicit time integration scheme has been proposed for solving the equation of motion. The equilibrium shape can be obtained naturally without nonlinear iterative correction and global stiffness matrix integration. Two classical curved surfaces of tension membranes produced under the uniform-stress condition are presented to verify the accuracy and efficiency of the proposed method.

## 1.  Introduction

Membrane structures have attracted considerable attention due to their light weight, high flexibility, and ease of stowing and deployment. The applications range from innovative architectural structures and space structures to medical appliances. A distinguishing feature of this type of structure is that it requires prestress to develop structural stiffness and to maintain a stable shape before any external load is applied. Therefore, for the structural function and strength satisfaction, a key step in the design of a membrane structure is to determine an initial configuration that respects the law of equilibrium for a feasible prestress distribution under given geometrical boundary conditions, especially for the gossamer space structures with high precision requirement of shapes (Jenkins, ).

The topic of shape analysis has become an attractive research field since the pioneering work of Otto Frei in the 1960s and 1970s (Otto and Rasch, ). Early studies were conducted by direct physical modeling for visualization. A representative case was the soap-film analogy (Otto and Rasch, ; Lewis, ), which had the natural property of adjusting to a stable minimal surface with uniform stress in every direction. These physical methods can be tedious and inaccurate for a measurement of the geometry, and involve error magnification when measurements of the model have to be scaled up to a full-size structure. Thus, with the advent of the computer age, many researchers have paid attention to the numerical methods for handling this subject. The analysis methods can be normally categorized into three main families: the stiffness matrix methods, the geometric stiffness methods, and the dynamic equilibrium methods (Veenendaal and Block, ). As a typical example of the stiffness matrix method, the traditional finite element method (FEM) was described in (Haug and Powell, ; Argyris et al., ). It may have some limitations due to the high nonlinearity of membrane systems. Apart from its demand for extensive matrix manipulations, possible singularity of stiffness matrix and divergence may be encountered when the trial geometry is far from the equilibrium position. However, with some special treatments and modifications (Wüchner and Bletzinger, ), it still can be used to solve this problem effectively. As an alternative for finding feasible equilibrium shapes, the geometric stiffness methods can simplify the nonlinear analysis. It is represented by the force density method derived from special discretization and linearization techniques (Schek, ). Although this method enables fast determination of equilibrium surfaces, it should be pointed out that additional iterative procedures are inevitable when it is used to model minimal surface forms (Maurin and Motro, ; Pauletti and Pimenta, ), and thus some of its attractiveness is lost. Unlike the preceding two methods, in the dynamic equilibrium methods (e.g., dynamic relaxation method (DRM)), a steady equilibrium shape for an initially unbalanced structure is achieved by simulating a damped vibration process over time (Lewis and Lewis, ). Particularly, combined with the kinetic damping technique (i.e., resetting the nodal velocity to zero at each kinetic energy peak), this method has been presented as a powerful tool capable of finding forms of tension structures modeled by simple truss networks, linear triangular meshes (Barnes, ; Lewis, ) or high-order elements (Gosling and Lewis, ). Nevertheless, due to the introduction of many hypotheses for simplification (e.g., line elements and triple force elements) and the ignorance of shear deformation, both some errors and the tendency of meshes to distort may induce some inconvenience in the traditional dynamic relaxation-based methods (Wood, ).

In recent years, some improved approaches based on the conventional techniques have also undergone significant development. For example, Brew and Lewis () proposed an alternative methodology to the experimentally physical modeling on the basis of the free hanging concept. Moreover, using particle-spring systems, Kilian and Ochsendorf () presented a novel approach for the exploration of funicular forms. This approach is also a kind of dynamic equilibrium method in nature. It is more suitable for treating the structural components with only axial deformation (e.g., cable or bar members), because for the solids, it is quite difficult to determine the spring stiffness that can represent the true material properties of the body exactly. Even if the body is discretized into equivalent networks, errors are still inevitable.

In addition to the above techniques and approaches which have been extensively applied to shape analysis, the finite particle method (FPM) (Yu, ) that is derived from vector mechanics and superior in the complex structural behavior analysis (Ting et al., ), can be another choice for this subject. In FPM, the analyzed domain is regarded as a body composed of a finite number of particles, instead of a mathematically continuous body adopted in the traditional methods based on analytical mechanics. The motion path of each particle can be modeled via a set of discrete path units. Within each path unit, Newton’s second law, rather than weak variational formulations for the governing partial differential equations (PDEs), is adopted to directly formulate the motion of the particles. Interactions of a particle with its neighbors are determined by the elements connected with them. In the calculation procedure, an incremental theory based on the concept of the convected material frame (Shih et al., ) and an explicit time integration scheme for solving the equations of motion are also adopted. For further specific considerations, the FPM based on Ting et al. ()’s mechanics concept was developed, and was successively applied to the analysis of kinematically indeterminate bar assemblies (Yu and Luo, ), deployable structures (Yu and Luo, ), and progressive failure of truss structures (Yu et al., ).

According to the analysis method classification aforementioned, the FPM belongs to the third category, i.e., the dynamic equilibrium methods. In FPM, neither the global material stiffness matrix nor the geometric stiffness matrix needs to be formed. In addition, there are some differences between the FPM and other dynamic equilibrium methods (e.g., DRM, particle-spring method), especially in the fundamental concepts and the procedure of interaction force evaluation. It should be noted that the FPM essentially employs a generalized mechanics frame which is capable of analyzing various complex behaviors. It is not a purely mathematical technique with introducing a fictitious mass and damping to the static equation to increase the numerical stability. Furthermore, in FPM, the convected material frames (i.e., the reference configuration) combined with other concepts, such as the path unit and the physical procedure of fictitious motion, are used to simplify the evaluation of interaction forces, while in some DR methods involving FE techniques, numerous matrix operations are required for calculating the finite rotation and deformation. Besides, the formulations of interaction forces in this method are rigidly derived from the theory of generalized large-deflection analysis of solids, instead of being a result of some simplifications (e.g., expressions in the form of element side/line forces or spring forces) as in some other DR methods. Therefore, a fine description of the internal force state can be obtained, even in some complex cases.

In view of the remarkable performance of the FPM in tackling various nonlinear behaviors of structures, as well as its particular convenience for evaluating the prestress effect on tension structures, the primary purpose of this work is to extend the application of the FPM using a recently-developed 3-node triangular (T3) membrane element to the shape analysis of tension membranes. The remainder of the paper is laid out as follows. Section 2 provides a brief introduction of the fundamentals. Then, in Section 3, techniques for handling the geometrically nonlinear behavior of the T3 facet membrane elements are illustrated, and the formulations of interaction forces are derived. In Section 4, a solution procedure of shape analysis based on the FPM is developed, together with an improved time integration technique to reduce the computational effort. Section 5 presents two classical examples to verify the effectiveness of the proposed method in the shape analysis of stable minimal surfaces of tension membranes. Finally, some concluding remarks of this work are provided in Section 6.

## 2.  Fundamental theory of the FPM

FPM is an intrinsic-modeling-based method designed to describe large rigid body motion and geometrical changes of a system of multiple continuous bodies. For simplicity, in this section, only a brief summary of the fundamental theory is provided by taking a 3D membrane as an example. More details can be found in (Ting et al., ; Yu, ; Yu et al., ).

### 2.1.  Point description

Considering a continuous membrane as shown in Fig. 1, the motion of the membrane can be characterized by the trajectories of a finite number of discrete particles. Assume that the shape of a membrane undergoes a successive transformation from the initial configuration V 0 to the reference configuration Va , and then to the current configuration V. At the initial time t 0, both the position x 0 and velocity v 0 of an arbitrary particle should be predefined. The subject to be addressed is converted to tracing the motion of each particle under the initial condition.

Fig.1
Point description model of a membrane

An important concept in FPM is that particles are adopted as main variables to describe the parameters of a body, such as mass, geometry, and deformation, i.e., so-called point descriptions. At any moment every particle is in an eternal dynamic equilibrium state under resultant forces, whereas elements connecting particles to their neighbors do not carry any mass, and maintain static equilibrium during the motion. Hence, it ensures a strong form of energy balance for the entire system. Besides, particles are mutually independent, and thus all force calculations can be implemented in each element separately.

### 2.2.  Path unit description

Another characteristic of the FPM is the employment of a path unit. Referring to Fig. 2, it is assumed that an arbitrary particle α undergoes a motion during the analysis period, say from time instant t 0 to tn , with its position vector correspondingly changing from x 0 to x n . Such a trajectory can be split into n-discrete segments identified by a series of time points t 0, t 1, t 2, ta , tb , …, tn . Thus, the motion path of each particle within a segment, e.g., ta ttb , is referred to as a path unit.

Fig.2
Discrete path units of a particle

Since the time interval is very small within each path unit, all variables (e.g., the constitutive law of materials, mass and topological relationships) can be regarded as constants with changes only at the boundary of a path unit, namely ta or tb . Consequently, the configuration of a deformable body at time ta can be considered as the reference material frame for the evaluation of interaction forces.

### 2.3.  Governing equations of motion

As the foregoing assumptions in FPM, every particle is considered to be in a dynamic equilibrium state at any instant. Considering an arbitrary particle α, as illustrated in Fig. 3, its motion within each path unit is determined by the equation that follows Newton’s second law, that is $${{\mathbf{M}}_\alpha }{{\mathbf{\ddot d}}_\alpha } = {\mathbf{F}}_\alpha ^{ext} - {\mathbf{F}}_\alpha ^{\{int}}$$, where M α is the particle mass matrix, $${\mathbf{F}}_\alpha ^{ext}$$ represents the acceleration vector, and $${\mathbf{F}}_\alpha ^{ext}$$ and $${\mathbf{F}}_\alpha ^{\operatorname{int} }$$ denote the external force and internal force applied on the particle α, respectively.

Fig.3
Illustration of the forces acting on the particle α

The determination of particle mass is a modeling procedure depending on an individual’s considerations for practical problems, rather than the results of some mathematic techniques. In this study, the mass term M α is prescribed by the physical mass of the membrane, which can be defined as $${\mathbf{M}}_\alpha ^{} = {{\mathbf{m}}_\alpha } + \sum\limits_{i = 1}^{{n_c}} {{\mathbf{m}}_i^e}$$, where m α is the mass attached to the particle α, $${\mathbf{m}}_i^e$$ is the evenly distributed mass of the ith element jointing at the particle α, and n c is the number of elements connected with the particle α.

$${\mathbf{F}}_\alpha ^{\operatorname{int} }$$ represents the internal forces of the particle α with its neighbors and it can be computed through an assembling procedure, that is $${\mathbf{F}}_\alpha ^{int} = \sum\limits_{i = 1}^{{n_c}} {{\mathbf{f}}_i^{int}}$$, where $${\mathbf{f}}_i^{int}$$ is the equivalent internal force contributed from the ith element connected with the particle α. The internal force evaluation is the most important work in FPM, which will be presented at length in Section 3.

$${\mathbf{F}}_\alpha ^{ext}$$ is the sum of the prescribed external forces acting on the particle α, that is $${\mathbf{F}}_\alpha ^{ext} = {{\mathbf{P}}_\alpha } + \sum\limits_{i = 1}^{{n_c}} {{\mathbf{f}}_i^{ext}}$$, where P α is the concentrated force directly applied at the particle α, and $${\mathbf{f}}_i^{ext}$$ is the equivalent force transformed from the distributed loads of the ith element according to the principle of virtual work.

### 2.4.  Description of kinematics

According to continuum mechanics (Bonet and Wood, ), the deformation at an arbitrary point of a membrane can be described by the geometrical change of a differential element with an infinitesimal size. As illustrated in Fig. 4a, at time ta the relative position vector of particle α to its neighborhood is denoted by d x a . In FPM, the configuration at time ta is taken as a convected material reference frame. After some time increments, at time t it changes to d x , that is $${\text{d}}{\mathbf{x}} = {\mathbf{F}}{\text{d}}{{\mathbf{x}}_a}$$, where F is the deformation gradient due to the resultant forces. In terms of the polar decomposition theorem, the deformation gradient can be decomposed into the product of the rigid body rotation matrix R and the pure deformation gradient U . Thus, Eq. (5) can be rewritten as $$d{\mathbf{x}} = {\mathbf{RU}}d{{\mathbf{x}}_a}$$.

Fig.4
Description of kinematics
(a) Different positions and configurations of an element during the motion; (b) Decomposition mechanism

It is worth noting that R is only theoretically existent, because the rigid motion and deformation are fully coupled and it is rather difficult to get the exact deformation. In fact, only an approximation to the real one with the same order of magnitude as that of the residual rotation will be sufficient to obtain a convergent solution. Besides, there is an inherited predictor-corrector mechanism in FPM that prevents divergence even if the evaluation of the interaction forces is not absolutely accurate. Thus, in this contribution, a simple scheme for the measure of the average rigid body rotation is developed, which will be explained in detail in Section 3. Referring to Fig. 4b, let particles along with the element undergo a fictitious reverse rotation, and then the motion and deformation of the differential element can be expressed as $${\text{d}}{{\mathbf{x}}_r} = {\mathbf{R}}_r^T{\text{d}}{\mathbf{x}} = ({\mathbf{R}}_r^T{\mathbf{R}}){\mathbf{U}}{\text{d}}{{\mathbf{x}}_a} = {{\mathbf{F}}_r}{\text{d}}{{\mathbf{x}}_a}$$, where R r represents an approximate average rotation, and F r denotes the deformation gradient between Va and Vr .

As illustrated in Fig. 4b, Vd is a fictitious configuration after the reverse rotation with the accurate rotation matrix R T. The difference between Vd and Vr is merely a residual rotation, viz. $$({\mathbf{R}}_r^T{\mathbf{R}})$$, the order of which should be equal to or less than that of a pure deformation. Thus, it is not necessary to take into account the effect of the residual rotation on the calculation of strain. Moreover, for the convenience of defining independent degrees of freedom (DOFs), a set of deformation coordinates according to F r are assumed. The corresponding strain, stress, and virtual work in the global coordinate system are $${\mathbf{E}} = {{\mathbf{E}}_r} = ({{\mathbf{U}}^2} - {\mathbf{I}})/2 = [{({{\mathbf{F}}_r})^T}{{\mathbf{F}}_r} - {\mathbf{I}}]/2$$, $${\mathbf{S}}({\mathbf{E}}) = {{\mathbf{S}}_r}({{\mathbf{E}}_r})$$, $$\delta U = \int_{{V_a}} {tr({\mathbf{S}} \otimes \delta {\mathbf{E}})} {\text{d}}V = \int_{{{\hat V}_a}} {tr({{\mathbf{S}}_r} \otimes \delta {{\mathbf{E}}_r})} {\text{d}}V$$.

With the transformation matrix Q that relates the global coordinates to the deformation coordinates at time t, Eqs. (8)–(10) can be rewritten in the deformation coordinate system as $${\mathbf{\hat E}} = {\mathbf{QE}}{{\mathbf{Q}}^T} = {\mathbf{Q}}{{\mathbf{E}}_r}{{\mathbf{Q}}^T} = {{\mathbf{\hat E}}_r}$$, $${\mathbf{\hat S}} = {\mathbf{QS}}{{\mathbf{Q}}^T} = {\mathbf{Q}}{{\mathbf{S}}_r}{{\mathbf{Q}}^T} = {{\mathbf{\hat S}}_r}$$, $$\delta U = \int_{{{\hat V}_a}} {tr({\mathbf{\hat S}} \otimes \delta {\mathbf{\hat E}})} {\text{d}}V$$.

In addition, the relationship between the Cauchy stress and the second Piola-kirchhoff stress can be expressed as $${{\mathbf{\hat \sigma }}_r} = \frac{1}{{\det |{{{\mathbf{\hat F}}}_r}|}}{{\mathbf{\hat F}}_r}{{\mathbf{\hat S}}_r}{\mathbf{\hat F}}_r^T = {\mathbf{Q}}\frac{1}{{\det |{{\mathbf{F}}_r}|}}{{\mathbf{F}}_r}{\mathbf{SF}}_r^T{{\mathbf{Q}}^T} = {\mathbf{Q}}\frac{1}{{\det |{\mathbf{U}}|}}({\mathbf{R}}_r^T{\mathbf{R}}){\mathbf{US}}{{\mathbf{U}}^T}{({\mathbf{R}}_r^T{\mathbf{R}})^T}{{\mathbf{Q}}^T}$$. According to the principle of material frame indifference, the stress-strain expression should be formulated based on the pure deformation state. This requirement can be approximately satisfied when the strain, the stress and their relationships are defined at the fictitious state Vr which is very close to Vd , the exact state after pure deformation. Thus, if the geometric change is negligible as compared with the initial configuration at time ta , the incremental strain $$\Delta {\mathbf{\hat E}}$$ can be simplified using the following infinitesimal strain $$\Delta {\mathbf{\hat E}} = \Delta {{\mathbf{\hat E}}_r} \cong \Delta {\mathbf{\hat \varepsilon }} = \frac{1}{2}\left\{ {\frac{{\partial (\Delta {\mathbf{\hat u}})}}{{\partial {\mathbf{\hat x}}}} + {{\left[ {\frac{{\partial (\Delta {\mathbf{\hat u}})}}{{\partial {\mathbf{\hat x}}}}} \right]}^T}} \right\}$$, where $$\Delta {\mathbf{\hat u}}$$ is the incremental deformation displacement vector. Given that $$\Delta U \cong I$$ and $${\mathbf{R}}_r^T{\mathbf{R}} \cong I$$, a simple relationship for the incremental stress can be obtained from Eq. (14) $$\Delta {{\mathbf{\hat \sigma }}_r} \cong {\mathbf{Q}}\Delta {\mathbf{S}}{{\mathbf{Q}}^T} = \Delta {\mathbf{\hat S}}$$.

From Eq. (16), it can be concluded that complicated stress updates for satisfying the principle of objectivity are not required, and the incremental stress can be defined as the engineering stress. Then, the total second Piola-Kirchhoff stress at time t, defined in compliance with the material reference frame, can be expressed in the deformation coordinate as $${\mathbf{\hat S}} = {{\mathbf{\hat S}}_a} + \Delta {\mathbf{\hat S}} \cong {{\mathbf{\hat \sigma }}_a} + \Delta {{\mathbf{\hat \sigma }}_r}$$.

The actual Cauchy stress at the current state can be obtained from the total second Piola-Kirchhoff stress through a pure rotation, that is $${\mathbf{\sigma }} = {{\mathbf{R}}_r}{{\mathbf{\sigma }}_r}{\mathbf{R}}_r^T \cong {{\mathbf{R}}_r}{{\mathbf{Q}}^T}{\mathbf{\hat SQR}}_r^T = {{\mathbf{R}}_r}{\mathbf{SR}}_r^T$$.

Substituting Eqs. (15) and (17) into Eq. (10), the virtual strain energy at time t is $$\delta U = \int_{{{\hat V}_a}} {tr\left[ {{{{\mathbf{\hat \sigma }}}_a} \otimes \delta {\text{(}}\Delta {\mathbf{\hat \varepsilon }}{\text{)}}} \right]} {\text{d}}V + \int_{{{\hat V}_a}} {tr\left[ {\Delta {{{\mathbf{\hat \sigma }}}_r} \otimes \delta {\text{(}}\Delta {\mathbf{\hat \varepsilon }}{\text{)}}} \right]} {\text{d}}V$$.

In Eqs. (17) and (19), the Cauchy stress $${{\mathbf{\hat \sigma }}_a}$$ is a known quantity calculated by Eq. (18) in the previous path unit.

Note that in FPM, the evaluation of the virtual work is only used for deriving the reaction forces acting on particles from each element, instead of solving the equivalent weak form of the PDEs of the entire system in terms of variational principles as in the FEM (Shih et al., ).

## 3.  Interaction forces in membrane structures

In this section, a T3 facet membrane element is used to evaluate the interaction forces between particles within the membrane structures. Instead of adopting large strain tensors or simple coordinate transformations as in the common approaches, a physical modeling procedure is developed to handle the geometric nonlinearity.

For instance, consider an arbitrary T3 membrane element (1-2-3) with thickness h as shown in Fig. 5a, in which each node has three translational DOFs in the global coordinate system. Because the nodes of the element are rigidly connected with the particles, the position vectors of node i at times ta and t can be defined as $$x_i^a$$ and x i (i=1, 2, 3), respectively. Since the time interval between ta and t is very small, both the material properties and the configuration are assumed to be unchanged within this time segment. If the configuration of the element at time ta is chosen as the reference material frame, the displacement increment of node i between ta and t is $$\Delta {u_i} = {x_i} - x_i^a,i = {\text{1}},{\text{ 2}},{\text{ 3}}$$.

Fig.5
Motion of a 3D triangular membrane element
(a) Positions and displacement increments of the particles at three corners of an element; (b) Rigid body rotation and relative displacement increments of the particles

Note that the total nodal displacement increment Δ u i includes pure deformation, rigid body translation, and rotation. For computational expediency, an efficient kinematical strategy is suggested here to extract deformation. As illustrated in Fig. 5b, it is assumed that the total rigid body rotation can be divided into two components, namely out-plane rotation θ 1 and in-plane rotation θ 2.

Firstly, the angle of out-plane rotation can be conveniently measured by $${\theta _{\text{1}}} = {\text{arcsin}}(\left| {{n_a} \times n} \right|)$$, where n a and n denote the normal vectors of a membrane element at time ta and t, respectively. The corresponding rotation axial vector is $${{\mathbf{n}}_1} = \frac{{{{\mathbf{n}}_a} \times {\mathbf{n}}}}{{\left| {{{\mathbf{n}}_a} \times {\mathbf{n}}} \right|}}$$.

Secondly, to account for the in-plane rotation, the element (1-2-3) is moved to make the centroids Ca and C coincide (Ca and C denote the centroids of the elements at time ta and t, respectively), and then every node (1, 2, 3) of the element at time t is projected onto the plane (1 p , 2 p , 3 p ) with normal vector n a (Fig. 6a). Consequently, the projective nodes (1 p , 2 p , 3 p ) and the element nodes (1 a , 2 a , 3 a ) should be in a same plane. Then, the rotation angle with regard to node i can be evaluated according to the change of orientation for a line connecting node i with the centroid of the element, as shown in Fig. 6b, $${\theta _{2i}} = arc\sin \left( {\frac{{\Delta {\mathbf{x}}_{ci}^a \times \Delta {\mathbf{x}}_{ci}^p}}{{\left| {\Delta {\mathbf{x}}_{ci}^a \times \Delta {\mathbf{x}}_{ci}^p} \right|}} \cdot {{\mathbf{n}}_a}} \right){\text{, }}\;\;i = 1,{\kern 1pt} {\kern 1pt} {\kern 1pt} 2,{\kern 1pt} {\kern 1pt} {\kern 1pt} 3$$, where $$\Delta x_{ci}^a$$ and $$\Delta x_{ci}^p$$ are the position vectors of node i with relative to the centroids of the element at time ta and t, respectively, i.e., $$\left\{ {\begin{array}{*{20}{c}} {\Delta {\mathbf{x}}_{ci}^a = \Delta {\mathbf{x}}_i^a - \Delta {\mathbf{x}}_c^a,} \\ {\Delta {\mathbf{x}}_{ci}^{} = \Delta {\mathbf{x}}_i^{} - \Delta {\mathbf{x}}_c^{},} \\ {\Delta {\mathbf{x}}_{ci}^p = \Delta {{\mathbf{x}}_{ci}} - (\Delta {{\mathbf{x}}_{ci}} \cdot {{\mathbf{n}}_a}){{\mathbf{n}}_a},} \end{array}} \right.$$ $$\left\{ {\begin{array}{*{20}{c}} {\Delta {\mathbf{x}}_c^a = \frac{1}{3}\sum\limits_{i = 1}^3 {({\mathbf{x}}_i^a - {\mathbf{x}}_1^a)} ,} \\ {\Delta {\mathbf{x}}_c^{} = \frac{1}{3}\sum\limits_{i = 1}^3 {({\mathbf{x}}_i^{} - {\mathbf{x}}_1^{})} .} \end{array}} \right.$$

Fig.6
Measure of in-plane rotation angle
(a) Position vectors $$\Delta x_{ci}^p$$ and Δxci ; (b) Nodal rotation angle θ 2i in pane na

Then, the angle of in-plane rotation is approximately measured by the average of three node rotations, that is $${\theta _2} = \frac{1}{3}\sum\limits_{i = 1}^3 {{\theta _{2i}}}$$.

So far, the total rigid body rotation vector of the element (1-2-3) from ta to t can be expressed as $$\theta = \theta {{\mathbf{n}}_\theta } = {\theta _{\text{1}}}{n_{\text{1}}} + {\theta _{\text{2}}}{n_{\text{2}}}$$, where θ is the angle of total rigid body rotation, and n θ is the corresponding rotation axial vector, n θ =[l m n]T.

Next, assume that the element (1-2-3) translates first with a fictitious reverse displacement (−Δ u 1) to make the nodes 1 and 1′ coincide, and then undergoes a fictitious reverse rotation (− θ ) to arrive at the fictitious position Vr with the nodal labels changed to (1 r -2 r -3 r ), as shown in Fig. 7a. When the fictitious reverse motion is completed, the pure deformation of the analyzed element can be obtained from a comparison of the configurations between the fictitious state Vr (1 r -2 r -3 r ) and the reference state Va (1 a -2 a -3 a ). Therefore, the deformation displacement vector of node i from ta to t is $${\mathbf{\eta }}_1^d = 0$$, $${\mathbf{\eta }}_i^d = {\mathbf{\eta }}_i^{} + {\mathbf{\eta }}_i^r = (\Delta {{\mathbf{u}}_i} - \Delta {{\mathbf{u}}_1}) + {{\mathbf{R}}^*}( - \theta )({{\mathbf{x}}_i} - {{\mathbf{x}}_1}),i = 2,\;\;3$$, where η i and $${\mathbf{\eta }}_i^r$$ are the total displacement increments of node i relative to node 1 and the rigid body displacement increment due to the fictitious reverse rotation of node i, respectively, and $${{\mathbf{R}}^*}( - \theta ) = {\mathbf{R}}( - \theta ) - {\mathbf{I}} = \left[ {1 - \cos ( - \theta )} \right]{\Omega ^2}({n_\theta }) + \sin ( - \theta )\Omega {\text{(}}{n_\theta }{\text{)}}$$, where R is a rotation matrix, I is a 3×3 identity matrix, and Ω ( n θ ) is an angular velocity matrix expressed in terms of the components of the axial vector n θ (Belytschko et al., ), $$\Omega \;({n_\theta }) = \left[ {\begin{array}{*{20}{c}} 0\;{ - n}\;m \\ n\;0\;{ - l} \\ { - m}\;l\;0 \end{array}} \right]$$.

Fig.7
Evaluation of interaction forces
(a) Fictitious reverse translation (−Δu 1) and rotation (−θ); (b) Deformation displacements in the deformation coordinate system

To evaluate the work equivalent internal forces acting on particles, a shape function with the same form as that developed in FEM is introduced for describing the strain distribution within the element. While different from the total displacements employed in FEM, the nodal variables here only account for the deformation. Therefore, the redundant DOFs corresponding to the rigid body modes must be eliminated to reduce the total independent variables to the correct number. To this end, a set of deformation coordinate systems are defined, with the origin locating at node 1 and the $$\hat x$$ axis parallel to the deformation displacement vector of node 2, $${\mathbf{\eta }}_2^d$$ (also can be replaced by $${\mathbf{\eta }}_3^d$$), as shown in Fig. 7b. The basic vectors of $$\hat x$$, $$\hat y$$, and $$\hat z$$ axes are denoted by $${{\mathbf{\hat e}}_x}$$, $${{\mathbf{\hat e}}_y}$$, and $${{\mathbf{\hat e}}_z}$$, respectively, that is $${{\mathbf{\hat e}}_x} = \frac{1}{{\left| {{\mathbf{\eta }}_2^d} \right|}}{\mathbf{\eta }}_2^d,{\text{ }}{{\mathbf{\hat e}}_z} = {{\mathbf{n}}_a},{\text{ }}{{\mathbf{\hat e}}_y} = \frac{{{{{\mathbf{\hat e}}}_3} \times {{{\mathbf{\hat e}}}_1}}}{{\left| {{{{\mathbf{\hat e}}}_3} \times {{{\mathbf{\hat e}}}_1}} \right|}}$$.

Before the evaluation of the interaction forces can be implemented in the deformation coordinates, the deformation vector $${\mathbf{\eta }}_i^d$$ should be transformed to the new coordinate system first by $${{\mathbf{\hat u}}_i} = {\left[ {{{\hat u}_i}\;\;{{\hat v}_i}\;{\text{ }}{{\hat w}_i}} \right]^T} = {\mathbf{\hat Q\eta }}_i^d,\;\;\;\;\;i = 1,\;\;2,\;\;3$$, where $${\mathbf{\hat Q}} = {\left[ {{\mathbf{\hat e}}_x^{}{\text{ }}{\mathbf{\hat e}}_y^{}{\text{ }}{\mathbf{\hat e}}_z^{}} \right]^T}$$ represents the coordinate transformation matrix from the global coordinates (x, y, z) to the local deformation coordinates $$(\hat x,\;\hat y,\;\hat z)$$.

In the new coordinate system, six components (i.e., $${\hat u_1},{\text{ }}{\hat v_1},{\text{ }}{\hat w_1},{\text{ }}{\hat v_2},{\text{ }}{\hat w_2},{\text{ }}{\hat u_3},{\text{ }}{\hat v_3},{\text{ }}{\hat w_3}$$) are zero, which implies that six DOFs are deleted. The remaining three non-zero deformation components are independent, and for convenience, they can be expressed in a condensed form, i.e., $${{\mathbf{\hat d}}^*} = {\left[ {\begin{array}{*{20}{c}} {{{\hat u}_2}}\;{{{\hat u}_3}}\;{{{\hat v}_3}} \end{array}} \right]^T}$$. Thus, the deformation displacements at an arbitrary point within the membrane element can be described by a set of linear interpolation functions similar to those in FEM, $${\mathbf{\hat u}} = {{\mathbf{\hat N}}^*}{{\mathbf{\hat d}}^*} = \sum\limits_{i = 2}^3 {{{\hat N}_i}} {\mathbf{\hat u}}_i^*$$, where $${\mathbf{\hat u}}_i^* = {\left[ {{{\hat u}_i}\;\;{{\hat v}_i}} \right]^T}$$ is a condensed nodal deformation vector, and $${\hat N_i}$$ and $${{\mathbf{\hat N}}^*}$$ are the shape function and the corresponding matrix form, respectively, where $${\hat N_i} = ({\beta _i}\hat x + {\gamma _i}\hat y)/{\alpha _1}$$ (i=2, 3), $${\alpha _1} = {\hat x_2}{\hat y_3} - {\hat x_3}{\hat y_2},{\beta _2} = {\hat y_3},{\beta _3} = - {\hat y_2},{\gamma _2} = - {\hat x_3},{\gamma _3} = {\hat x_2}$$. With the deformation distribution functions, one can formulate the incremental strain and stress in terms of the elastic mechanics theory, and hence $$\Delta {\mathbf{\hat \varepsilon }} = {{\mathbf{\hat B}}^*}{{\mathbf{\hat d}}^*} = \frac{1}{{{\alpha _1}}}\left[ {\begin{array}{*{20}{c}} {{\beta _2}}\;{{\beta _3}}\;0 \\ 0\;0\;{{\gamma _3}} \\ {{\gamma _2}}\;{{\gamma _3}}\;{{\beta _3}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\hat u}_2}} \\ {{{\hat u}_3}} \\ {{{\hat v}_3}} \end{array}} \right]$$, $$\Delta {{\mathbf{\hat \sigma }}_r} = {{\mathbf{\hat D}}_a}\Delta {\mathbf{\hat \varepsilon }} = {{\mathbf{\hat D}}_a}{{\mathbf{\hat B}}^*}{{\mathbf{\hat d}}^*}$$, where $${{\mathbf{\hat D}}_a}$$ denotes the constitutive matrix of the material under the plane stress condition referring to the stress state, and $${{\mathbf{\hat \sigma }}_a}$$ is the stress value at time ta . If an isotropic elastic material is used, the matrix $${{\mathbf{\hat D}}_a}$$ can be written as $${{\mathbf{\hat D}}_a} = \frac{{{E_a}}}{{1 - {\nu ^2}}}\left[ {\begin{array}{*{20}{c}} 1\;\nu \;0 \\ \nu \;1\;0 \\ 0\;0\;{\frac{{1 - \nu }}{2}} \end{array}} \right]$$. Then, substituting Eqs. (32) and (33) into Eq. (19), one can obtain the equivalent internal forces, $${{\mathbf{\hat f}}^*} = {h_a}\int_{{A_a}} {{{({{{\mathbf{\hat B}}}^*})}^T}} {{\mathbf{\hat \sigma }}_a}{\text{d}}\hat A + ({h_a}\int_{{A_a}} {{{({{{\mathbf{\hat B}}}^*})}^T}} {{\mathbf{\hat D}}_a}{{\mathbf{\hat B}}^*}{\text{d}}\hat A){{\mathbf{\hat d}}^*}$$, where Aa and ha denote the area and the thickness of the membrane element at time ta , respectively.

The other three components of nodal forces corresponds to the DOFs that have been eliminated in the deformation coordinate system, and thus they should be determined by the static equilibrium conditions of the element in the $$\hat x - \hat y$$ plane, namely $$\sum {{{\hat M}_z} = 0, \; } {\hat f_{2y}} = ({\hat f_{3x}}{\hat y_3} - {\hat f_{2x}}{\hat y_2} - {\hat f_{3y}}{\hat x_3})/{\hat x_2}$$, $$\sum {{{\hat F}_x} = 0,\;\;\;\;} {\hat f_{1x}} = - ({\hat f_{2x}} + {\hat f_{3x}})$$, $$\sum {{{\hat F}_y} = 0,\;\;\;\;} {\hat f_{1y}} = - ({\hat f_{2y}} + {\hat f_{3y}})$$.

As a matter of fact, this equilibrium condition must be imposed to ensure that the internal virtual work remains unchanged even if subjected to a superposed rigid body motion, which is a requirement of the principle of objectivity for a continuum.

Next, to assemble the equivalent interaction forces acting on a conjunct particle (i.e., summed by Eq. (3)), all the force components must be transformed to the global coordinate system. In addition, the element needs to undergo a forward motion, including a translation (+Δ u 1) and a rotation (+ θ ), and to move back to the position at time t. Therefore, the actual internal forces at time t can be determined by $${\mathbf{f}}_i^{int} = {\left[ {\begin{array}{*{20}{c}} {{f_{ix}}}\;{{f_{iy}}}\;{{f_{iz}}} \end{array}} \right]^T} = {\mathbf{R}}( + \theta ){{\mathbf{Q}}^T}{{\mathbf{\hat f}}_i}{\text{,}}\;\;\;\;i = 1,\;\;2,\;\;3$$, where $${{\mathbf{\hat f}}_i} = {\left[ {\begin{array}{*{20}{c}} {{{\hat f}_{ix}}}\;{{{\hat f}_{iy}}}\;0 \end{array}} \right]^T}$$ can be derived from Eqs. (35) and (36), R is a rotation matrix from Eq. (28) except that the rotation angle should be replaced by (+θ), and Q is a coordinate transformation matrix as that in Eq. (30).

In the above formulations, magnitudes and directions of internal forces depend on the estimated rigid body rotation that may be not an exact value. A correction measure in displacements and rotations is thus generated in the subsequent calculation. In other words, if a particle moves too far from the equilibrium position due to the inaccurate evaluation in the previous step, the errors may exert effect on the calculation in the following step, which will increase the internal force in the opposite direction and pull the particle back to the equilibrium position. In short, the estimated rotation values, the equilibrium requirements of the elements and the equations of motion compose a self-predictor-corrector system in a natural way, which prevents the divergence due to the error accumulation.

## 4.  Computational procedure of shape analysis based on the FPM

### 4.1.  Solution strategy for governing equations

Since every item in the governing equations of motion is explicitly expressed (Eqs. (1) to (4)), the structural response under given initial conditions may be explicitly obtained by replacing the nonlinear iterations with the recursions over time. Particularly, if a damping item is added to Eq. (1) to dissipate the kinetic energy of the system, the unbalanced force acting on each particle will be gradually reduced. In other words, every particle can naturally approach the equilibrium position and finally come to rest under the damping effect when the oscillations are progressively reduced. However, choosing a proper damping that is able to attenuate the pseudo-dynamic response in the shortest time is not easy. Moreover, in this study we mainly focus on the accuracy of the steady-state solution, while the transient process has no physical meaning and thus is not concerned. Therefore, to obtain fast convergence to the equilibrium shape, a modified integral-form time integration strategy is suggested here.

Firstly, instead of directly solving the second-order force equation, we convert Eq. (1) into a first-order one (i.e., so-called a momentum equation of motion (Chang et al., )) by integrating it with respect to time once, i.e., $${{\mathbf{M}}_\alpha }{{\mathbf{\dot d}}_\alpha } = {\mathbf{\bar F}}_\alpha ^{{\text{ext}}} - {\mathbf{\bar F}}_\alpha ^{{\text{int}}}$$, where $${\mathbf{\bar F}}_\alpha ^{{\text{ext}}} = \int {{\mathbf{F}}_\alpha ^{{\text{ext}}}} {\text{d}}t$$ and $${\mathbf{\bar F}}_\alpha ^{{\text{int}}} = \int {{\mathbf{F}}_\alpha ^{{\text{int}}}} {\text{d}}t$$ represent the time integral of the external force and the time integral of the internal force, namely the external momentum and the internal momentum, respectively. In the case of shape analysis, no external forces are applied, so $${\mathbf{\bar F}}_\alpha ^{{\text{ext}}} = 0$$.

Then, assume that every particle satisfies a special condition similar to that of an initial value problem as follows: $${{\mathbf{\bar d}}_n} = 0$$, $${{\mathbf{\bar d}}_{n - 1}} = 0$$, where $${{\mathbf{\bar d}}_{n - 1}}$$ and $${{\mathbf{\bar d}}_n}$$ are used to denote the approximations to $${\mathbf{\bar d}}({t_{n - 1}})$$ and $${\mathbf{\bar d}}({t_n})$$, the theoretical integrals of the displacements from t 0 to tn- 1 and from t 0 to tn , respectively. This condition implies that the displacement can be considered to be zero before time tn , namely d n =0 and d n −1=0. Thus, the calculation of the integral function of the displacement at the following step, $${\mathbf{\bar d}}({t_{n + 1}})$$, is independent of the previous trajectory history. In the physical perspective, it suggests that the motion of the particles can be deemed to re-start from rest at each step. In this manner, the pseudo dynamic response of the discretized system can be effectively attenuated.

It should be mentioned that the direct integration methods can be employed to computer the solutions of the new governing equation, and their basic properties in stability, accuracy, and convergence remain unchanged. If an explicit central difference formulation is adopted, an approximation to the velocity of an arbitrary particle α at time tn can be given by $${{\mathbf{\dot d}}_n} = ({{\mathbf{\bar d}}_{n + 1}} - 2{{\mathbf{\bar d}}_n} + {{\mathbf{\bar d}}_{n - 1}})/\Delta {t^2}$$.

Substituting Eqs. (39) and (40) into Eq. (38) and rearranging terms gives the time integral of displacement at time tn explicitly, $${{\mathbf{\bar d}}_{n + 1}} = \Delta {t^2}{\mathbf{M}}_\alpha ^{ - 1}({\mathbf{\bar F}}_n^{{\text{ext}}} - {\mathbf{\bar F}}_n^{{\text{int}}})\; \cdot S = - \Delta {t^2}{\mathbf{M}}_\alpha ^{ - 1}{\mathbf{\bar F}}_n^{{\text{int}}} \cdot S$$.

Then, the displacement of particle α at time tn+ 1 can be approximated in a finite difference form as $${{\mathbf{d}}_{n + 1}} = \frac{{({{{\mathbf{\bar d}}}_{n + 1}} - {{{\mathbf{\bar d}}}_{n - 1}})}}{{2\Delta t}} = - \frac{{\Delta t}}{2}{\mathbf{M}}_\alpha ^{ - 1}{\mathbf{\bar F}}_n^{{\text{int}}} \cdot S$$.

In Eqs. (41) and (42), S is a positive parameter for accelerating the convergence procedure. The physical meaning of factor S is illustrated in Fig. 8. With the aid of this scale parameter, it can be expected that the rate of convergence will be improved. After a large number of trial calculations, it is also found that an empiric value of S on the order of 102 is appropriate when a similar rate of convergence can be obtained.

Fig.8
Illustration of the fast solution strategy

The resulting formulation (Eq. (42)) seems similar to that of Euler’s method since they both employ the explicit difference format, but here we involve some special physical and mechanical considerations, rather than a pure mathematical technique for solving the differential equations of motion.

It also can be found that the proposed time integration strategy only involves simple algebraic manipulations for vectors and at each step the solution may proceed particle by particle in a numerically uncoupled form. As a result, no global stiffness matrices need to be formed, and the requirement for the memory storage is very low.

### 4.2.  Procedure for initial stable shape analysis

During the shape analysis implemented by the FPM, the calculation of equilibrium is always made based on a set of convected material frames, adopting a constant prescribed stress σ 0 as the Cauchy stress (in Eq. (35)) at the beginning of each path unit. Moreover, we set the fictitious elastic modulus of the membrane model to be a very small value to ensure that the stress distribution in the equilibrium configuration is close enough to that of the specified prestress. However, the value of the modulus should not be zero, otherwise the self-error-correction mechanism (whereby the geometric change of an element can be suppressed by the interaction forces in the opposite direction induced by deformation) will become invalid, possibly leading to surface inaccuracy and irregular locations of particles that may be unusable for the subsequent load analysis in the design. After providing an appropriate initial arrangement of particles, we start the solution for the equilibrium shape of a membrane structure under given boundary conditions in which the displacements of the particles attached to the supporting systems are constrained. A recurrent calculation has to be carried out until the system reaches a steady equilibrium configuration according to the convergence criteria for the out-of-balance forces and the stress deviations, given by $$\frac{{{{\left\| {{\mathbf{F}}_n^{{\text{int}}}} \right\|}_2}}}{{{{\left\| {{\mathbf{F}}_0^{{\text{int}}}} \right\|}_2}}} \leqslant {\varepsilon _{\text{F}}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{{\left\| {{\mathbf{\sigma }}_n^{} - {\mathbf{\sigma }}_0^{}} \right\|}_2}}}{{{{\left\| {{\mathbf{\sigma }}_0^{}} \right\|}_2}}} \leqslant {\varepsilon _{\text{S}}}$$, where the tolerances ε F and ε S are set to be very small values, so as to minimize the admissible influence of the dynamic response.

The whole solution procedure outlined above for finding initial equilibrium shapes of tension membranes is summarized in the flowchart in Fig. 9. A self-designed C++ program has been developed to implement the proposed method that has been used for the numerical tests in Section 5.

Fig.9
Procedure of shape analysis performed by the FPM

## 5.  Examples

To illustrate the performance of the proposed methodology, two well-known examples are described in this section: a catenoid and a Scherk-like surface. Both cases preserve the principle of uniform and isotropic tensile stress for driving the membrane towards an ideal equilibrium shape that minimizes the surface area under given boundary conditions, as observed in the soap film models (Lewis and Lewis, ; Wüchner and Bletzinger, ; Pauletti and Pimenta, ). The self-coded program for the numerical tests is run on a desktop PC with a CPU of Inter Core processor i5-2400 3.1 GHz and a memory of 4.00 GB RAM. All the results are obtained with the fixed tolerance values of ε F=ε S=10−3 for the convergence criteria.

### 5.1.  Catenoid

The catenoid is the only minimal surface of rotational symmetry that can be described analytically, and thus it is commonly used as a benchmark shape-analysis case (Lewis and Lewis, ; Wood, ; Ye et al., ). The theoretical equation of a catenoid surface is given by $$z = h - a\left[ {\ln (\sqrt {{x^2} + {y^2}} + \sqrt {{x^2} + {y^2} - {a^2}} ) - \ln a} \right]$$. In this example, a=12 m, h=27.5092 m are used and the sketch for the curvature shape of the catenary along the cross section is shown in Fig. 10.

Fig.10
A sketch for the curvature shape of a catenary

The initial configuration of the analysis model is a conoid shape with fixed constraints imposed on the particles lying along the edges of two parallel rings with 12 and 60 m in radius and a distance of h apart to satisfy the boundary conditions determined by Eq. (44). Its geometrical and material properties are described in Fig. 11. A state of isotropic prestress, i.e., σ 0=[σx σy τxy ]T=[σ 0 σ 0 0]T with σ 0=20 kN/m (force per unit width), is re-applied during the analysis process in order to obtain a close approximation to the minimal surface. Given the symmetric condition of the problem, only a quarter of the entire membrane is solved. To monitor the trend of convergence, three cases with different densities of particles (i.e., 49, 100, and 196 particles) are considered. An artificial scale factor S=100 for increasing the accumulative displacements, a time step Δt=10−4 s and a fictitious elastic modulus $$\tilde E = 1\;\;kN/{m^2}$$ are employed for these numerical implementations. According to the given convergence criteria, the program respectively takes 182, 278, and 394 steps (corresponding to the three different models) to complete the solutions.

Fig.11
Catenoid: initially assumed analysis model

Fig. 12 illustrates the stable-found shape of each model after analysis. It is noticed that all the particles are reasonably distributed across the final surfaces with regular spacing. The surface areas are 12 253.217 m2, 12 186.362 m2, and 12 149.755 m2, respectively, giving errors within 1.5% of the theoretical value of 12 118.302 m2. At the same time, the maximum root mean square (RMS) error (in the model composed of 49 particles) in the z-direction is only 0.431 m, i.e., 1.567 % of the height.

Fig.12
Form-found shapes of tension membranes in the form of catenoid
(a) Model 1 with 49 particles; (b) Model 2 with 100 particles; (c) Model 3 with 196 particles

The geometrical profiles of the sections of final surfaces across the diameter of the upper/lower ring are compared against the analytical solution in Fig. 13. It can be observed that the discrepancy is small even in the case of a coarse arrangement of particles, and the convergence trend is apparent as the surface discretization is refined.

Fig.13
Profiles of the form-found shapes in comparison with the theoretical solution

For comparison, DRM with kinetic damping (Lewis and Lewis, ) is also programmed and applied to this example. Note that the same level of discretization is employed by the DRM. The results of the comparative study on accuracy of FPM against DRM are displayed in Table 1, where the coordinates of the points marked with 1–m observed in Fig. 12 are presented and the errors are evaluated against the exact positions obtained from the theoretical equation. It can be found that both methods can provide accurate solutions. The results predicted by the FPM are a little closer to the theoretical values, though the improvement is very limited. However, the DRM has to monitor the kinetic energy of the entire system as the iterations proceed. Moreover, an extra effort for estimating the precise point of maximum kinetic energy and corresponding corrections for displacements are repeatedly required whenever an energy peak occurs. For the FPM, by contrast, it is not necessary to take into account any additional condition during the whole analysis process. Furthermore, with the aid of the fast integration strategy for static solutions, the FPM can handle the shape analysis of tension membrane structures more effectively, which will be further investigated in the next example.

#### Table 1

Comparison of the membrane surface coordinates of the FPM with the exact solutions
 Point No. Radius, r (m) Z-component (m) FPM error (%) DRM error (%) FPM results DRM results Analytic values 1 12.461 24.210 24.215 24.195 0.063 0.086 2 13.836 20.985 20.990 20.953 0.153 0.179 3 15.982 18.028 18.034 17.985 0.241 0.275 4 18.753 15.356 15.367 15.310 0.302 0.371 5 22.041 12.954 12.964 12.904 0.391 0.468 6 25.733 10.797 10.811 10.750 0.439 0.564 7 29.732 8.873 8.883 8.825 0.543 0.658 8 34.048 7.113 7.121 7.068 0.640 0.753 9 38.669 5.489 5.496 5.450 0.726 0.848 10 43.563 3.985 3.991 3.954 0.781 0.942 11 48.764 2.573 2.579 2.552 0.803 1.038 12 54.305 1.235 1.238 1.224 0.877 1.136

• Note: the numerical results listed above are corresponding to the fairly fine analysis model with a density of 196 particles

• ### 5.2.  Scherk-like surface

Like the catenoid, the Scherk-like minimal surface (or referred to as a Box surface) is also a classic shape-analysis case, which has been previously studied by Maurin and Motro (), Brew and Lewis (), and Yuan et al. (). In this example, the initial configuration of the analysis model is a box shape which consists of three rectangular planes with the straight boundaries following the edges of a cube. The geometrical and material properties of the box model are specified in Fig. 14. To achieve the stable minimal surface, an isotropic and uniform prestress with σ 0=20 kN/m and a scale factor S=100 are also employed. The fictitious elastic modulus is set as $$\tilde E = 1\;\;kN/{m^2}$$ and the size of time step is Δt=5×10−4 s.

Fig.14
Scherk-like surface: initially assumed analysis model

In the case of a dense arrangement of particles (425 particles), the final equilibrium surface spanning the rigid boundary frame is represented in Fig. 15. The final area of this surface is 246.28 m2, which shows good coincidence with the exact geometry, as that obtained by Yuan et al. ().

Fig.15
Form-found shape of tension membrane in the form of Scherk-like minimal surface

Here, since the theoretical solution is not known, the results generated by DRM based on Barnes’s formulations (Barnes, ) are used for the assessment of accuracy. Figs. 16a and 16b show the profiles of the sections of the surfaces across the x-z and y-z central planes, respectively. A visual comparison shows that the discrepancy between the results predicted by FPM and those generated by DRM is not apparent. It also can be found that the center point at the top of the profiles lies very closely to midway between the top and bottom boundaries, with a slight deviation of only 0.536% of the height range, as necessitated by the zero mean curvature required for the minimal surface and the symmetry of the boundaries as well.

Fig.16
Profiles of the sections of the Scherk-like surface
(a) Across the x-z central plane; (b) Across the y-z central plane

For a better investigation of the solution process, the variations of kinetic energy (KE) of the entire system with the iteration history and the shapes of the membrane structure at different steps are plotted in Fig. 17. As predicted, the KE becomes smaller and eventually approaches zero, as more computation steps are taken in the implementations of both methods. In details, it is noticed that the decrease trends of KE during the solution process are almost the same for both methods, but with the FPM, a much smoother convergence process for attenuating the pseudo-dynamic oscillations can be achieved within a shorter iteration history, when compared with the DRM with kinetic damping.

Fig.17
Convergence process for kinetic energy (KE) of the entire system

Table 2 gives the results concerning the relative efficiency of the two methods for increasing DOFs, in terms of the number of iterations for convergence, the amount of time per iteration, and the total runtime required by FPM and DRM, respectively. One can observe that the slightly more computational effort of the FPM per step can be compensated well by the fairly decreased number of iterations for completing the solutions. Obviously, the FPM using the proposed formulations combined with the fast time integration strategy exhibits a higher convergence rate, so that a maximum reduction of solution time up to around 40% has occurred while a result of equal accuracy still can be achieved as compared to the DRM. Based on these results, it can be reasonably deduced that the FPM has a potential for an efficient analysis of large and complex structures with sufficiently fine discretization for ensuring accuracy. This statement is also documented in a comparison of the performance between the two methods as shown in Fig. 18, where the relationship between the influence of the level of discretization input and the time taken for the calculation to terminate is further illustrated. This comparison reveals that the time cost increases approximately log-linearly with the DOFs for both methods, but at a little lower rate for the FPM. To sum up, the proposed method can mimic the solutions generated by the DRM very well while considerably reducing the computational costs, especially when the size of the problem is enormous.

#### Table 2

Comparison of the efficiency between FPM and DRM
 No. of DOFs DRM FPM Improvement, $$\frac{{DRM - FPM}}{{DRM}}$$ (%) RMS deviation in z-direction (m) No. of iterations Average runtime per iteration (s) Total runtime (s) No. of iterations Average runtime per iteration (s) Total runtime (s) 288 592 0.0042 2.483 447 0.0046 2.051 17.4 3.46×10−1 399 648 0.0063 4.075 481 0.0068 3.270 19.8 1.82×10−1 528 697 0.0089 6.204 503 0.0097 4.894 21.1 1.09×10−1 675 745 0.0122 9.079 524 0.0133 6.988 23.0 6.81×10−2 1023 887 0.0203 18.029 554 0.0218 12.098 32.9 3.95×10−2 1275 993 0.0264 26.215 561 0.0281 15.792 39.8 2.73×10−2

Fig.18
Relationship between computational time and degrees of freedom for both numerical methods

## 6.  Conclusions

This study presents a general FPM framework for investigating the shape analysis of tension membrane structures under given boundary conditions. An efficient numerical procedure of generating viable initial equilibrium shapes is given. If this method is applied iteratively with constant and isotropic stress prescribed at any given point on the surface, it can lead to a stable minimal surface. The analysis of two classical tension membranes characterized by minimal surfaces has been conducted using a program following the proposed formulations and solution strategies.

In FPM, the rigid body motion can be simply decomposed from the total displacement increments by fictitious reverse motion without complicated and time-consuming matrix operations, which makes it convenient to deal with nonlinear problems. At the same time, the modified explicit time integration strategy which is expressed in the form of momentums in place of forces, can efficiently reduce the computational effort in search of a steady equilibrium shape. Besides, during the entire analysis, no global stiffness matrices need to assembled or solved, and no iterations to follow nonlinear laws are necessary as well. Moreover, the effect of prestress can be readily exerted on the particles, and the final equilibrium shape can be naturally obtained as the limit of the pseudo oscillations without introducing any special treatment or modification to cope with the common numerical instabilities. Thus, it has inherent advantages in the shape analysis of tension membranes.

From the numerical examples, the ability of the proposed method to obtain the accurate minimal surfaces of tension membrane structures has been verified by the comparisons with the analytical solutions and the numerical solutions of the DRM. It also has been proved that the FPM is much more efficient, in terms of the total CPU time required for convergence, than the DRM. In addition, no significant mesh distortion that may lead to poor representation of a surface occurs in this study. All these facts indicate a great promise for developing FPM into a general and effective method for shape analysis and further shape optimization of large-scale complex tensile structure systems that have potential applications in civil and aerospace engineering.

## Acknowledgements

The authors would like to express their gratitude to Prof. Edward TING of Purdue University, USA for his advice and encouragement during the preparation of this work.

* Project supported by the National Natural Science Foundation of China (Nos. 51025858 and 51178415)

## References

[1] Argyris, J.H., Angelopoulos, T., Bichat, B., 1974. A general method for the shape finding of lightweight tension structures. Computer Methods in Applied Mechanics and Engineering, 3(1):135-149.

[2] Barnes, M., 1999. Form finding and analysis of tension structures by dynamic relaxation. International Journal of Space Structures, 14(2):89-104.

[3] Belytschko, T., Liu, W.K., Moran, B., 2000.  Nonlinear Finite Elements for Continua and Structures. Wiley,New York :544-549.

[4] Bonet, J., Wood, R.D., 1997.  Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press,Cambridge, UK :57-113.

[5] Brew, J.S., Lewis, W.J., 2003. Computational form-finding of tension membrane structures: non-finite element approaches: Part 1. Use of cubic splines in finding minimal surface membranes. International Journal for Numerical Methods in Engineering, 56(5):651-668.

[6] Brew, J.S., Lewis, W.J., 2007. Free hanging membrane model for shell structures. International Journal for Numerical Methods in Engineering, 71(13):1513-1533.

[7] Chang, S., Tsai, K., Chen, K., 1998. Improved time integration for pseudodynamic tests. Earthquake Engineering & Structural Dynamics, 27(7):711-730.

[8] Gosling, P.D., Lewis, W.J., 1996. Optimal structural membranes-I. Form-finding of prestressed membranes using a curved quadrilateral finite element for surface definition. Computers & Structures, 61(5):885-895.

[9] Haug, E., Powell, G.H., 1972. Finite element analysis of nonlinear membrane structures. , IASS Pacific Symposium Part II: on Tension Structures and Space Frames, Tokyo and Kyoto, Japan, 93-102. :93-102.

[10] Jenkins, C.H.M., 2001. Gossamer spacecraft: membrane and inflatable structures technology for space applications.  Progress in Astronautics and Aeronautics. AIAA,Reston, VA 191:

[11] Kilian, A., Ochsendorf, J., 2005. Particle-spring systems for structural form finding. Journal of the International Association for Shell and Spatial Structures, 46(2):77-84.

[12] Lewis, W.J., 2003.  Tension Structures: Form and Behaviour. Thomas Telford,London :

[13] Lewis, W.J., Lewis, T.S., 1996. Application of formian and dynamic relaxation to the form-finding of minimal surfaces. Journal of the International Association for Shell and Spatial Structures, 37(3):165-186.

[14] Maurin, B., Motro, R., 1998. The surface stress density method as a form-finding tool for tensile membranes. Engineering Structures, 20(8):712-719.

[15] Otto, F., Rasch, B., 1995.  Finding Form: Towards an Architecture of the Minimal. Edition Axel Menges,Fellbach, Deutschland :

[16] Pauletti, R.M.O., Pimenta, P.M., 2008. The natural force density method for the shape finding of taut structures. Computer Methods in Applied Mechanics and Engineering, 197(49-50):4419-4428.

[17] Schek, H.J., 1974. The force density method for form finding and computation of general networks. Computer Methods in Applied Mechanics and Engineering, 3(1):115-134.

[18] Shih, C., Wang, Y.K., Ting, E.C., 2004. Fundamentals of a vector form intrinsic finite element: Part III. Convected material frame and example. Journal of Mechanics, 20(2):133-143.

[19] Ting, E.C., Shih, C., Wang, Y.K., 2004. Fundamentals of a vector form intrinsic finite element: Part I. Basic procedure and a plane frame element. Journal of Mechanics, 20(2):113-122.

[20] Veenendaal, D., Block, P., 2012. An overview and comparison of structural form finding methods for general networks. International Journal of Solids and Structures, 49(26):3741-3753.

[21] Wood, R.D., 2002. A simple technique for controlling element distortion in dynamic relaxation form-finding of tension membranes. Computers & Structures, 80(27-30):2115-2120.

[22] Wchner, R., Bletzinger, K.U., 2005. Stress-adapted numerical form finding of pre-stressed surfaces by the updated reference strategy. International Journal for Numerical Methods in Engineering, 64(2):143-166.

[23] Ye, J., Feng, R., Zhou, S., 2011. The modified dynamic relaxation method for the form-finding of membrane structures. Advanced Science Letters, 4(8-10):2845-2853.

[24] Yu, Y., 2010.  Progressive Collapse of Space Steel Structures based on the Finite Particle Method. PhD Thesis, (in Chinese), Zhejiang University,Hangzhou, China :

[25] Yu, Y., Luo, Y.Z., 2009. Finite particle method for kinematically indeterminate bar assemblies. Journal of Zhejiang University-SCIENCE A, 10(5):667-676.

[26] Yu, Y., Luo, Y.Z., 2009. Motion analysis of deployable structures based on the rod hinge element by the finite particle method. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 223:156-171.

[27] Yu, Y., Luo, Y.Z., Paulino, G.H., 2011. Finite particle method for progressive failure simulation of truss structures. ASCE Journal of Structural Engineering, 137(10):1168-1181.

[28] Yuan, S., Liu, X., Ye, K., 2010. FEMOL solution for minimal surface form finding of tensile membrane structures. China Civil Engineering Journal, (in Chinese),43(11):1-7.