膜構造研究論文集 第30号, 9-14, 2016年
Research Paper on Membrane Structures, No.30, 9-14, 2016

Nonlinear finite element procedure for membrane structures using ETFE films
considering friction contact condition

Quang Hieu BUI 1*1
Masaya KAWABATA 2*2

Synopsis
In frame-supported membrane structures type, one of techniques to introduce the pre-stress on the membranes is extending theboundary. In this case, the friction contact between the membrane and the frame-supported, such as steel arch, was observed.
Therefore, the analytical method of membrane structures considering the friction contact is essential.
This paper describes a nonlinear finite element procedure, which includes the non-linear geometry of membrane elements andfriction contact condition. In this analysis, the friction contact between membrane elements and frame-supported elements wasreplaced by Gap elements. The proposed method will be applied to investigate the sliding between ETFE film and steel arch of theparallel steel arch model under forced displacement condition. The accuracy and applicability of proposed procedure wasconfirmed by an experiment.

are the Lagrange multiplier method and Penalty method. The later
method were attracted researchers to develop the Gap/Friction
element. Contact element stiffness matrices are symmetric for
frictionless contact problems. Unfortunately, problems involving
sliding friction result in an unsymmetrical stiffness matrix whichmakes the problems in convergence criteria. In 1995 Ju et al suggested
the symmetrical stiffness matrix for 3D frictional contact element
which can be implemented into an existing nonlinear finite element
program to solve the nonlinear friction contact [4][5].
In this paper, the FEM procedure which includes the non-linear
geometry of membrane elements and nonlinear friction contact
problem using Ju’s method, is proposed. The proposed method will
be applied to investigate the sliding between ETFE film and steel archof the parallel steel arch model under forced displacement condition.
The accuracy and applicability of proposed procedure was confirmedby an experiment.

1. Introduction
Recently, ETFE film is well used in membrane structures because
its advantages, such as transparent, light weight, easy to fabricate alarge-sized panel by heat-sealing and excelling in the durability overultraviolet rays. In term of design of ETFE membranes structure,since the ETFE film possesses no flexural stiffness, its shape, the loadon structure and the internal stresses interact in a nonlinear manner tosatisfy the equilibrium equations. Therefore, various kinds of analysis
procedure have been proposed to solve those geometrical nonlinear
problems. One of them was the finite element technique. The
procedure of this technique was described clearly in many previous
reports [1][2].
In term of suspension and frame-supported ETFE membranestructures, one of techniques to introduce the pre-stress on themembranes is extending the boundary. In this case, the friction contact
between the films and the frame-supported was observed. Therefore,
the analytical method of membrane structures considering the friction
contact is essential. Computationally, the analysis of frictional contact
problems is extremely difficult, even for static cases involving thesimplest constitute relation [3]. The reasons could be explained as theboundary conditions are unknown as a priori and the contacting
surfaces change in size and shape as load is applied. There had two
main approaches to solve the contact problem in FEM analysis, which*1*2

2. Finite element formulation
2.1 Membrane elements
The triangular membrane elements with three nodes and ninedegrees of freedom per element were used in this paper. The
equilibrium equation for a single element in local coordinate system
may be obtained via the principle of virtual work, as follows:

Graduate student, Graduate School of Urban Innovation, Yokohama National University, Japan
Associate Prof., Graduate School of Urban Innovation, Yokohama National University, Japan9∫ δ𝛜 𝛔dV − δ𝐮 𝐩 = 𝟎TT

here, the matrix M can be also referenced in [2].
2.2 Gap elements
A Gap element is considered to connect nodes A and B. The local
coordinate xyz of this element is shown in Fig. 1. The normal
direction of Gap element is assumed in z direction, while the slidingdirection, when Gap element is in sliding status, is assumed in x
direction. The y direction is determined from x and z directions

(2.1)

Ve

where, 𝛔: stress, 𝐩: the external node force vector, 𝛜: strain, u: thenode displacement vector, 0: zeroes vector and Ve : the element
volume. The relation between the strain 𝛜 and node displacement
vector u can be described as follows:𝛜 = 𝐁0 𝐮 + 0.5𝐀𝛉; δ𝛜 = (𝐁0 + 𝐀𝐆)δ𝐮xrxA(2.2)
A

here, the matrices 𝐁0 , A, 𝛉, G can be referenced in [2]. In this paper,
the large displacements and small strains are considered, so theconstitute relations for plane stress analysis may be used as:

yrxB
rzA
B

ryA

rzBzryB

Fig.1 The local coordinate of Gap elements
𝛔 = 𝐃𝛜 + 𝛔0(2.3)
The relation between the internal node force vector andincremental displacement vector for Gap element via the penalty
method can be determined as follows:

where, 𝛔0 denotes the initial stress vector and D is the elastic matrix
defined as:
1𝐸𝜗𝐃=[1 − 𝜗20𝜗
10001 − 𝜗]
2𝐫 = 𝐤 𝐠 ∆𝐮here, ∆𝐮 = [∆𝑢𝑨 ∆𝑣𝑨 ∆𝑤𝐴 ∆𝑢𝑩 ∆𝑣𝑩 ∆𝑤𝑩 ]𝐓 is
incremental node displacement vector, in which ∆u, ∆v and ∆w are
incremental displacements in x, y and z directions, respectively; 𝐫 =[𝑟𝑥𝐴 𝑟𝑦𝐴 𝑟𝑧𝐴 𝑟𝑥𝐵 𝑟𝑦𝐵 𝑟𝑧𝐵 ]𝐓 is the internal node force
vector of Gap element; and kg is Gap element stiffness matrix anddetermined as equation (2.9) [4]:here, E and 𝜗 are the Young’s modulus and Poisson’s ratio.
Substituting (2.2) and (2.3) into (2.1) and eliminating δ𝐮T , wehave:∫ (𝐁0 + 𝐀𝐆)T [𝐃(𝐁0 𝐮 + 0.5𝐀𝛉) + 𝛔0 ]dV − 𝐩 = 𝟎(2.8)(2.4)

Ve𝐤𝐤𝐠 = [−𝐤The above equations must be transformed to global coordinatesand finally assembled to obtain the global equilibrium equations.
Since the global equations will be solved iteratively by NewtonRaphson method, the linearize process of the governing equation, atthe element level, can be described after the ith iteration as follow:
𝐤 𝒊𝒎 ∆𝐮𝒊 = 𝐩 − 𝐫 𝒊𝑘𝑠where, 𝐤 = [ 00𝜇 2 𝑘𝑛
status and 𝐤 = [ 0𝜇𝑘𝑛

(2.5)

(2.6)

Veand the element tangent stiffness matrix, ∂𝐫 𝑖 /𝜕𝐮, consists of twoparts as (2.7):
𝐤 𝑖𝑚 = ∫ (𝐁0 + 𝐀𝒊 𝐆)T 𝐃(𝐁0 + 𝐀𝒊 𝐆)dV + ∫ 𝐆T 𝐌𝑖 𝐆dVVe(2.9)

00 ] is used when Gap element is in stick𝑘𝑛0𝑘𝑠0𝜇𝑘𝑛
0 ] is used when Gap element is𝑘𝑛in sliding status; kn is an input large penalty constant for normal
direction of Gap element; ks is the penalty constant for sliding
directions and  is the static friction coefficient.
In sliding status, because the assumption of sliding direction is x
direction, the internal friction force 𝑓𝑠 = 𝑟𝑥𝐴 = −𝑟𝑥𝐵 and thenormal force 𝑓𝑛 = 𝑟𝑧𝐴 = −𝑟𝑧𝐵 obey the Mohr-Coulomb friction
theory (𝑓𝑠 = 𝜇𝑓𝑛 ) according to the equations (2.8)&(2.9); however
the increment of normal force 𝑓𝑛 contains the term, 𝜇𝑘𝑛 (∆𝑢𝑨 −∆𝑢𝑩 ). If this term is omitted from the calculation of the total normalforce, the contact analysis from equations (2.8)&(2.9) will produce asimilar result at the convergence state (∆𝑢𝑨 − ∆𝑢𝑩 ≅ 0).

here, the node force vector after ith iteration ri can be determined by(2.6):
𝐫 𝒊 = ∫ (𝐁0 + 𝑨𝒊 𝐆)T [𝐃(𝐁0 𝐮𝒊 + 0.5𝐀𝐢 𝛉𝐢 ) + 𝛔0 ]dV0𝑘𝑠
0−𝐤]𝐤(2.7)

Ve10here, 𝑟𝑥 = 𝑟𝑥𝐴 = −𝑟𝑥𝐵 and 𝑟𝑦 = 𝑟𝑦𝐴 = −𝑟𝑦𝐵 are the friction
forces in x and y directions and 𝑟𝑧 = 𝑟𝑧𝐴 = −𝑟𝑧𝐵 is the normal
force in z direction of Gap elements in stick condition, which aredetermined by Eq. (2.8)(g) Check the convergence for the incremental displacement ∆𝐔 𝒊(h) If step (g) converges, then output the results(i) If step (g) does not converge, then return to step (c) for the next
iteration.

2.3. Membrane structures include Gap elements
The principle of virtual work for the membrane structure include
Gap elements can be expressed as:(2.10)

V

where, 𝝐𝑚 and 𝝈𝑚 are the strain and stress of membraneelements; Rg is internal force of Gap elements; U is the whole
structure node displacement vector and P is the external node vector.
After eliminating δ𝐔 T and solving the global equation (2.10) by
Newton-Raphson method, the governing equation after the ith
iteration can be described as:(2.11)

here, 𝐊 𝒊𝒎 and 𝐑𝒊𝒎 are the tangent stiffness matrix and internal
force vector of membrane elements, which calculated by assembly
the elements stiffness in equation (2.7) and the element internal force
in equation (2.6) respectively; 𝐊 𝒊𝒈 and 𝐑𝒊𝒈 are the tangent stiffness
matrix and internal force vector of Gap elements, which calculated byassembly the elements stiffness in equation (2.9) and the element
internal force in equation (2.8) respectively depend on the status ofelements.
The displacements after (i+1)th iteration can be determined as:

1000

This updating of the node point displacements in the iteration iscontinue until the incremental displacements, ∆𝐔 𝒊 , are small ornearly zeroes. This convergence criteria will satisfy the condition ofusing the symmetric stiffness matrix of Gap elements in case ofsliding status.

1000
ETFE

3. Analysis procedures
The analysis procedures for membrane structures include Gapelements are listed below:(a) Input data(b) Assume all Gap elements are in stick condition for the first loop of
Newton-Raphson method(c) Determinate the structure stiffness matrix at step i: 𝐊 𝒊𝒎 +𝐊 𝒊𝒈(d) Determinate the internal force vector at step i: 𝐑𝒊𝒎 + 𝐑𝒊𝒈(e) Determinate the incremental displacement vector by Eq. (2.11)(f) Check the Gap elements status: the Gap elements will be assumedin stick status until equation (2.13) is reached. And it is assumed thatif Gap elements change into sliding status, this status will be remained
till the end of Newton-Raphson loop:√𝑟𝑥2 + 𝑟𝑦2 − μ|𝑟𝑧 | = 0

X

54
32
1

1000

(2.12)

1000

Y

2000

𝐔 𝒊+𝟏 = 𝐔 𝒊 + ∆𝐔 𝒊

4000
2000

1000

4000
2000

1000

Steelarch

OutsideframeH(𝐊 𝒊𝒎 + 𝐊 𝒊𝒈 )∆𝐔 𝒊 = 𝐏 − 𝐑𝒊𝒎 − 𝐑𝒊𝒈

4. Analysis examples
This numerical example is based on the experiments of new
method of form-finding for ETFE film structure.
4.1. Experiments overview
For the tension type of ETFE membrane structure, the films arepulled outward from the exterior to introduce membrane tension. In
those experiments, the stretching boundary of ETFE film with twoparallel steel arch model were took place. The dimension ofspecimens are described as Fig. 2. There are three type of arch height
H (200mm, 300mm and 400mm), which determined depend on thespan of outside frame.

4x200∫ (δ𝝐𝑚 )T 𝝈𝑚 dV + δ𝐔 T 𝐑𝑔 − δ𝐔 T 𝐏 = 𝟎ZXFig.2 Specimens overview (all dimensions in millimeter)First, the outside frames and steel arches were fixed first. After
that, the flatting ETFE film was put on the steel arches and outsideframe. After the ETFE film was fixed with outside frame, the forced
displacement was input in Z direction step by step until getting thetarget values. During the experiments, the movement of offset points
on ETFE films in X and Y direction, which shown in ¼ model ofspecimens in Fig. 2, were observed and were used for the forced
displacement condition in analysis. Those value can be referred in [6].
In addition, during the experiments, the sliding between the ETFE
films and steel arches were also observed at special points (from 1 to(2.13)

11

yields an oscillated numerical solution; however the small penalty
constant typically results in easy convergence, but the numericalresults might be less reliable. In this paper, the recommendation of S.H.
Lee [7] for those value were used. The normal stiffness penalty value
kn was get as three high order than ETFE films, and the lateral stiffness
penalty value ks was get as 0.1 kn.

Table 1. Experiments results on sliding observation
Sliding distances on arch direction (mm)200_Spec
300_Spec
400_Spec

Observed points
10.0
0.0
1.0

2
-1.0
6.0
8.0

3
2.5
11.0
15.0

4
4.5
16.5
21.0

5
6.5
21.0
27.0

1000

G2800

* : Fixed condition
 : Fixed in Y direction
 : Fixed in X direction
 : Gap elements

600

4.2. Input data of analysis examples
In this analysis examples, the triangular elements will be used tomodel the ETFE membrane. And the friction contact between ETFEfilms and steel arch was modeled by Gap elements. The formulationsof those elements were described in section 2 and section3. The
general view of this analysis will be shown in Fig. 3. However
because of the symmetric, ¼ analysis model is used to save the timeof calculation.

400

Y Direction

Specimens

31*5) in Fig. 2. Those experiments results were shown in Table 1 at thefinal step of experiments

G4G6G822*

200

00500

1*

1000

1500

X Direction(mm)

2000

21*
Z direction(mm)400300200
10000500

1000

1500

X direction (mm)2000

Fig. 4. Initial nodes coordinates for the modeling in case of300mm Arch model

zG1

yxG9Fig 3. The general view of analysis model
The Fig. 4 shows the dimension of ¼ analysis model, the positionof Gap elements and the boundary conditions. In this model, twohundreds triangular ETFE membrane elements were used. Thecontacts between ETFE membrane and steel arch were replaced bynine Gap elements at arch position. The input forced displacements atfixed nodes (1*31*) were determined from the experiments [6].
4.3. The assumption of local coordinate of Gap elements
In this analysis, the normal direction (z axis) of Gap element isassumed from the center of the arch to contact points. The sliding
direction (x axis) is assumed as the tangent direction of arch in point
contacts. The y direction is determined from x and z direction. This
local coordinate of Gap elements is assumed unchanged during thecalculation. The Fig. 5 illustrates the local coordinate of Gap elements.
4.4 Results and discussions
For penalty method, the large penalty constant which can achieve
convergence is strived to find. The very large penalty constant usuallyZYXFig. 5. The assumptions of local coordinate of Gap elements
The equivalent stresses on membrane elements and thedeformation shape, when static friction coefficient was chosen as 0.2
in case of 300mm arch model, were shown in Fig. 6. As it can be seenfrom Fig. 6, the sliding between the ETFE film and steel arch can be
only seen in arch direction. The reason can be explained as theassumption of local coordinate of Gap elements in section 4.3.
Although this results is different with experiment results, thisassumption still be used for the next investigation because the sliding

12corresponding with G2, G4, G6 and G8 Gap elements in analysis inFig. 4.

in arch direction is mainly observed in experiments.

Arch position

Fig. 8 The sliding results of experiments and analysis
As we can see in Fig. 8, the analysis results and experiment results
were almost the same in case of 200mm and 300mm arch models. In
case of 400mm arch model, the large difference between analysis andexperiment was observed. The reason can be explained as the sliding
in analysis was observed in x local direction of Gap elements as thetangent direction of arch, while in the experiment the sliding value
was measured on the curve of arch.

Fig. 6 The equivalent stress on membrane elements (N/m2)
A various values of static friction coefficients were used in order
to illustrate the accuracy and applicability of above proposedprocedure. In case of 300mm arch model, the figure 7 shows thesliding values between ETFE films and steel arch at 9 Gap elements
positions in Fig 2. As we can see in Fig. 7, the larger friction
coefficient the smaller sliding values. Those analysis results werecompletely coincide with the physically phenomenon.

G1G2G3G4G5G6G7G85. Conclusions
In this research, the friction contact of ETFE film and steel arch isboth investigated by experiments and analysis. The experiments show
the large effect of friction contact on the ETFE films during theextending boundary construction work. The finite element procedure,
which combined the nonlinear geometry of membrane elements andnonlinear boundary condition of friction contact, therefore, wasproposed. The node-to-node contact elements or Gap elements whose
symmetric stiffness matrix could be used to model the friction contact
in nonlinear geometry analysis. The comparison between theexperiments and analysis illustrated the accuracy and applicability ofthe proposed method. However, the assumption of local coordinate ofGap elements should be studied more in the future in order to have abetter evaluation for the friction contact.

G9

Acknowledgements-The support for the experiments from Taiyo
Kogyo Corporation is gratefully acknowledged.

Fig. 7 The sliding results in x direction of Gap elements300mm arch modelReferences[1] Kazuo Ishii: State-of-the-art report on the stress deformation
analysis of membrane structure, Research report on membranestructures 1990, Japan, pp69-105.
[2] B. Tabarrok, Z. Qin: Nonlinear analysis of tension structures,Computers & Structures Vol. 45, No. 5/6, pp. 973-984, 1992
[3] K.J. Bathe, A. Chaudhary: A solution method for planar andaxisymmetric contact problems, Int. J. Numer. Meth. Engng 21, 65-In many technical reports for the properties characteristics ofETFE films, the friction coefficient between ETFE films and steel
was described as 0.2 [8]. Therefore, this value was used for the next
investigation. The figure 8 shows the sliding values in arch direction
in both experiments and analyses. The dash line shows theexperiments results while the solid line shows the analysis results.
Four observed point (25) in experiments in Fig. 2 were13

[6] Eisuke Ishida: The study on form finding via extending the outerboundary of parallel arch model using ETFE films, Bachelor Thesis(in Japanese only)[7] S.H. Lee: Rudimentary considerations for adaptive Gap/Friction
element based on the penalty method, Computers & Structures, Vol47, No.6, pp. 1043-1056, 1[8] The technical report for ETFE films, Asahi Glass Co..

88, 1985[4] S.H. Ju, J.J. Stone, R.E. Rowlands: A New symmetric contact
element stiffness matrix for frictional contact problems, Computers &Structures, Vol 54, No.2, pp. 289-301, 1995
[5] S.H. Ju, R.E. Rowlands: A Three-dimensional frictional contact
element whose stiffness matrix is symmetric, Journal of Applied
Mechanics, Vol. 66, pp. 460-467, 1999

摩擦接触条件を考慮した ETFE 膜を用いた膜構造の非線形有限要素法
BUI Quang Hieu*1)
河端 昌也*2)
梗 概
フレームに支持された膜構造タイプでは、膜にプレストレスを導入する技術の一つが境界を延伸することである。この場合、膜と支持されたフレームの摩擦接触が観察された。したがって、摩擦接触を考慮した膜構造の解析手法が不可欠である。
本論文では、膜要素の非線形幾何と摩擦接触条件を含む非線形有限要素法について述べる。この分析では、膜要素とフレーム要素との間の摩擦接触を Gap 要素で置き換えた。提案手法を適用して、並行アーチでの筋アーチと ETFE フィルムとの間のすべりを強制変位条件下で調べる。提案された手順の精度と適用性は、実験によって確認された。

*1） 横浜国立大学大学院 都市イノベーション学府 大学院生
*2） 横浜国立大学大学院 都市イノベーション研究院 准教授

14