The Finite Strip Method (FSM) 1. Introduction This is the method of semi-numerical and semi-analytical semi-analytical nature. It is suitable for the analysis of rectangular plates and plane-stress elements or structures being the combination of both. Therefore, the following types of civil engineering structures can be dealt with: bridge slabs, box girders, roofs consisting of plane elements, etc. However, the method has one important drawback, which limits its versatility – the analysed elements must be rectangular and simply-supported at two opposite edges.
bridge slabs
roofs of industrial halls
2. Analysis of plates by FSM The plate is considered as thin, i.e. it is assumed that the variation of displacements across the plate thickness in negligible and the plate is replaced with a plane surface corresponding to the mid-surface of the real 3D body. Also the in-plane deformations are neglected. Hence, at each point of the mid-surface the deflection and its derivatives with respect to t o two in-plane co-ordinates x and and y are are enough to define the deformation state. 2 2–2 y q (x ,y )
1
1 q (x ,y )
I
2 x 1–1
q (x ,y )
The co-ordinate system is introduced with the axis x along along two opposite simply-supported edges. The two edges parallel to the axis y may may have any type of support. In the FSM the plate is divided into the finite number of strips running along the axis y, thus each of the strips is simply-supported. The strips are connected along the so called nodal lines. All the strips may have any width, not necessarily the same one. The starting point for the derivation of the method procedures is the approximation of displacement function w I (x ,y ) for the points on a single strip I . Let us denote the nodal lines along this strip as i and j and j .
i j
L w i
w j w i
b Nodal line
Finite strip
ϕ i
w j
ϕ j
Polynomial
Sine series
The deflection of the strip is approximated as a combination of the sine harmonic series in the longitudinal direction y (analytical aspect) and the polynomial function f m in the transverse direction x (numerical aspect) r
r m π y m π y w (x , y ) = ∑ f m (x ) sin = ∑ (Am + B m x + ...) sin L L m =1 m =1 I
where: r is the adopted number of harmonic functions in the series and Am , B m, etc. are the coefficients of the polynomial function f m, which must be found from the boundary conditions corresponding to the deflections and transverse slopes at the nodal lines i and j . The continuity conditions along the nodal lines require the adjacent strips to have the same deflections and slopes. For an arbitrary strip and the m -th harmonic function we have four prescribed nodal values w im , ϕ im , w jm , ϕ jm , so the required polynomial must be of the third order
f m (x ) = Am + B m x + C m x 2 + D m x 3 The nodal values of displacements and slopes are r
w i =
∑ w im sin m =1 r
ϕ i =
∑ ϕ im sin m =1
r
m π y L
w j =
m π y L
ϕ j =
w jm sin ∑ m =1
r
∑ ϕ jm sin m =1
m π y L m π y L
The unknown coefficients in the polynomial are obtained from the boundary conditions x = 0 ⇒
f m (x ) = w im
and
x = b ⇒
f m (x ) = w jm
and
∂f m (x ) ∂x ∂f m (x ) ∂x
= ϕ im
= ϕ jm
and after some reordering one gets:
f m (x ) = C 0 i w im + C 1i ϕ im + C 0 j w jm + C 1 j ϕ jm where the coefficients C are
x 2 x 3 C 0 i = 1 − 3 2 + 2 3 b b x 2 x 3 C 0 j = 3 2 − 2 3 b b
x 2 x 3 + C 1i = x − 2 b b 2 x 3 x 2 C 1 j = 2 − b b
Note, that these functions are identical with the shape functions N 2, N 3, N 5 and N 6 in the clampedclamped beam.
The deflection approximation can be given in the matrix form w I (x , y ) =
r
I I C b w bm sin ∑ m 1 =
m π y L
with the following definitions of the vectors
w im ϕ im I w bm = w jm ϕ jm
C b = [C 0i C 1i C 0 j C 1 j ] I
The superscript I denotes the I -th strip and the subscript b – the bending state. Note, that the nodal parameters w im , ϕ im , w jm , ϕ jm , are actually the amplitudes of m -th harmonic functions describing the deflections and transverse slopes along the nodal lines. Having specified the deflection approximation in terms of nodal parameters w im , ϕ im , w jm , ϕ jm , we can consider the energy and derive the equilibrium conditions for the strip. The total energy for one strip I consists of the strain energy and the energy of loading I
I
U I = U s + U q
The strain energy can be expressed in terms of bending and torsional moments, M x, M y, M xy and the corresponding curvatures L b
I U s
1 = 20
2 I 2 I 2 I ∂ w ∂ w ∂ w − M x dxdy − M + 2M xy y 2 2 ∂ ∂ x y ∂ ∂ x y 0
∫∫
while the loading energy involves the deflections L b
I U q
=−
∫ ∫ q (x , y )w dxdy I
0 0
In the matrix form the strain energy can be written as L b
T 1 = MI κ I dxdy 200
∫∫
I U s
where the vectors of moments and curvatures
∂ 2w I − 2 ∂2x I ∂ w I κ = − ∂y 2 2 I 2 ∂ w ∂x ∂y
M x MI = M y M xy
were introduced. After the substitution of deflection w I the curvatures vector takes the form κ
I
r
=
∑B m =1
where the matrix Bbm I is
I I bm w bm
− C 1i ' ' sin k m y − C 0 j ' ' sin k m y − C 1 j ' ' sin k m y − C 0i ' ' sin k m y 2 I 2 2 2 B bm = k m C 0 i sin k m y k m C 1i sin k m y k m C 0 j sin k m y k m C 1 j sin k m y 2k m C 0i ' cos k m y 2k m C 1i ' cos k m y 2k m C 0 j ' cos k m y 2k m C 1 j ' cos k m y
m π L Let us assume that the plate material is isotropic. Then from the theory of isotropic plates the moments can be found and k m =
2 I ∂ 2w I ∂ w M x = − D 2 + D 1 2 ∂y ∂x 2 I ∂ 2w I ∂ w M y = − D 2 + D 1 2 ∂ ∂ y x
M xy = 2D xy
2
w I ∂x ∂y
∂
where the plate stiffness, the coupled stiffness and the torsional stiffness Eh 3 D = 12(1 − ν 2 ) D 1 = ν D D xy =
1 − ν D 2
for the plate of the thickness h were introduced. The vector of moments can be now expressed in the matrix form I
MI = D b κ I
with the matrix of bending stiffness coefficients for the isotropic plate
D D 1 0 I D b = D 1 D 0 0 0 D xy The moment vector can also be expressed in terms of the nodal displacement parameters r
I
M =
∑D
I I I b B bm w bm
m =1
and its transpose is M
I T
r
=
∑
I T
I T
I
w bm B bm D b
m =1
Hence, the strain energy takes the form I U s
1 r r I T I T I I I w bn Bbn D b B bm dxdy w bm = 2 m =1 n =1 0 0
∑∑
L b
∫∫
The double integral will involve the following expressions
0 for m ≠ n = sin k y sin k ydy n ∫ m L 2 for m = n 0 L
L
0 for m ≠ n cos k m y cos k n ydy = L 2 for m = n 0
∫
Hence, the single sum is sufficient in this expression and we get I U s
1 r I T I T I I I w bm B bm Db B bm dxdy w bm = 2 m =1 0 0 L b
∑
∫∫
The loading energy takes the form I U q
r
=−
∑
I T w bm
m =1
L b
∫∫
I T
C bm q (x , y ) sin k m ydxdy
0 0
The energy of the entire plate U p is the sum of energies for all the N -strips N
U p =
∑ U
I
I =1
In the presented energy approach the theorem is used saying, that the system is in the equilibrium when the entire potential energy has the minimal value. This leads to the conditions of vanishing partial derivatives of U p calculated with respect to the nodal displacement parameters in wbm I for every m -th harmonic function ∂U p ∂w 0m
=
0,
∂U p ∂ϕ 0 m
=
0,
∂U p ∂w 1m
=
0,
∂U p ∂ϕ 1m
=
0, …
∂U p ∂w Nm
=
0,
After this differentiation the set of equations is obtained for every m N
∑
I I K bm w bm
I =1
N
=
∑P
I bm
I =1
where the stiffness matrix for a single strip I and m -th harmonic function L b
I K bm
=
∫∫
I T
I
I
Bbm Db B bm dxdy
0 0
and the loading vector for this strip and this harmonic function L b
I Pbm
=
∫∫
I T
C bm q (x , y ) sin k m ydxdy
0 0
were introduced. The explicit form of the strip stiffness matrix is
k b 1 k b 3 k b 2 I K bm = sym .
k b 4 k b 5 k b 1
− k b 5
k b 6
− k b 3
k b 2
with: k b 1 = 6
13 12 L 2 6L 2 L 4 D + Lbk m D + k m D xy + k m D 1 3 70 5 b 5 b b
∂U p ∂ϕ Nm
=
0
k b 2 = 2
L 1 4 2 4 2 2 D + Lb 3 k m D + Lbk m D xy + Lbk m D 1 b 210 15 15
k b 3 = 3
11 1 3 L 4 2 2 D + Lb 2 k m D + Lk m D xy + Lk m D 1 2 420 5 5 b
k b 4 = −6
L 9 12 L 2 6L 2 4 + − − D Lbk D k D k m D 1 m m xy 140 5 b 5 b b 3
k b 5 = −3
L 13 1 1 4 2 2 D + Lb 2 k m D − Lk m D xy − Lk m D 1 2 840 5 10 b
k b 6 =
L 1 1 1 4 2 2 D − Lb 3 k m D − Lbk m D xy − Lbk m D 1 280 5 30 b
The loading vectors depend on the type of loading. For instance: - the point load y P 0 L x 0 y 0 b
x
2 3 x 0 x 0 1 − 3 2 + 2 3 b b 2 3 x 0 x 0 x 0 − 2 b + b 2 I Pbm = P sin k m y 2 3 0 3 x 0 − 2 x 0 b 2 b 3 2 3 − x 0 + x 0 2 b b
- the patch load y q 0
L x 1 x 2
y 1 y 2
b
x 3 x 4 x − 2 + 3 b 2b 2 x 2 x 3 1 x 4 2 − 3 b + 4 b 2 I Pbm = q 0 c m 3 4 x x − b 2 2b 3 3 4 1 x 1 x + − 3 b 4 b 2
x c m =
1 (cos k m y 1 − cos k m y 2 ) k m n
n
x n = x 2 − x 1 for n = 1, 2, 3 and 4 The assembly of the stiffness matrix Kbm p and the loading vector Pbm p for the entire plate divided into N -strips is carried out according to the following scheme w 0m , ϕ 0m 0 1 w 1m , ϕ 1m 1 2 w 2m , ϕ 2m 2 3 w 3m , ϕ 3m 3
N
N
w Nm , ϕ Nm
Kbm p =
w 0m ϕ 0m w 1m ϕ 1m w 2m ϕ 2m w 3m ϕ 3m
Kbm 1 Kbm 2 Kbm Pbm p =
w 0m Pbm ϕ 0m w 1m Pbm 2 ϕ 1m w 2m ϕ 2m w 3m ϕ 3m Pbm 3
Kbm N
Pbm
w Nm ϕ Nm
w Nm ϕ Nm
Note, that the directions of axes of local co-ordinates in the strips coincide with global co-ordinates and no transformation is necessary. As for the boundary conditions, the simple supported opposite edges are already inherent in the system in the form of the sine series functions, which fulfil the conditions of vanishing displacements w and bending moments (second derivatives of w with respect to y ). On the other hand, the support conditions on two remaining edges must be introduced. If any of the edge parameters is vanishing, then the corresponding rows and columns in Kbm p and the corresponding element of Pbm p can be removed or replaced with zeros. For instance, if the edge along the 0-nodal line is simply supported, then the displacement w 0m vanishes and the first row and the first column in Kbm p must be modified as well as the first element of Pbm p . In this way the global set of equilibrium equations for the entire plate for the m -th harmonic function is obtained p
p
p
K bm w bm = Pbm
The solution of these equations provides the vector of amplitudes of sine functions for deflections and transverse slopes along all the nodal lines for the m -th harmonic function. The value of the displacement at an arbitrary point of the plate is obtained by a summation of results for all the assumed r -harmonic functions according to the formula w (x , y ) = I
r
∑ C b I w bm I sin m =1
m π y L
The method ensures the continuity of deflections and slopes between the strips along the nodal lines. However, due to the approximate form of the displacement function in the strips, the bending and torsional moments calculated using the second derivatives of displacements are not continuous. The approximate values of moments along the nodal lines can be obtained as mean values computed from the moments yielding from two adjacent strips. For instance, for the nodal line i lying between the strips I −1 and I we get
0 1
1 I −1 I Mright + Mleft 2 where: Mi =
1 I −1
i
I
I −1 Mright
N
r
=
∑D
)
I −1 I −1 b B bm
(x = b )w bm I −1
m =1 I
N
(
Mleft =
r
∑D m =1
I I b B bm
(x = 0 w ) bm I
and the appropriate matrices Bbm have the form
6 − b 2 sin k m y I −1 Bbm (x = b ) = 0 0
−
2 6 sin k m y sin k m y b b 2 2 0 k m sin k m y 0
4 6 sin sin k m y k y m b 2 b 2 I B bm (x = 0 ) = k m sin k m y 0 0 2k m cos k m y
0
−
−
4 sin k m y b 0
2k m cos k m y
6 2 sin sin k y k y m m b b 2 0 0 0
0
The differences between the moments obtained from the adjacent strips decrease with the increasing number of strips. Note, that generally the accuracy of the results obtained using the FSM depends on two parameters: the number of strips N and the number of the harmonic functions r . 3. Analysis of plane stress elements (plate-like elements loaded in their plane) by FSM We consider thin plane elements, i.e. it is assumed that the variation of displacements across the element thickness in negligible and the element is replaced with a plane surface corresponding to the mid-surface of the real 3D body. Due to the existence of only in-plane loading the displacements are also only in-plane. Thus, we have two displacement functions u and v . Likewise in the plate analysis, the element is divided into a finite number of strips, which span the entire length of the element between two simply-supported opposite edges. y u = 0, σ y = 0 p (x )
1
2
3
I
N
p (y ) L
v (x ,y ) u (x ,y ) u = 0, σ y = 0 x
0
1
2
Let us introduce the strain and stress vectors
∂u ε x ∂x ∂v ε = ε y = ∂y γ xy ∂u ∂v + ∂y ∂x
σ x σ = σ y τ xy
and write down the physical law for the plane stress element σ
= D p ε
i
j
N
where the plane stress stiffness matrix has the general form
D p =
E x ν x E x 1 E y ν E (1 − ν x ν y ) y y 0 0
0 (1 − ν x ν y )E xy 0
For the isotropic material E x = E y = E , ν x = ν y = ν , E xy =
E = G 2(1 + ν )
The boundary conditions for the simply supported opposite edges y = 0 and y = L are u = 0 , σ y = 0 The approximation of displacements in a single strip has the form m π y m π y r r sin L sin L u (E + F m x + ...) f m (x ) = v = m π y m m =1 cos m π y m =1 cos L L
∑
∑
where E m, F m, etc. are the unknown coefficients. The assumed form of displacement functions fulfils the boundary conditions at the simply supported edges. The presence of the function sin(m πy L) in the displacement u ensures the fulfilment of the displacement boundary conditions. / As for the stress boundary conditions we have
σ y = ηε x + ξε y = η
∂u ∂x
+ ξ
∂v ∂y
The differentiation with respect to x in the first term does not change the function sin(m πy / L), while the differentiation with respect to y in the second term transforms the function cos(m πy L) into / / sin(m πy L). Hence, both the terms contain the desired function sin(m πy / L) and the conditions for the vanishing σ y are also automatically fulfilled. The calculation of displacements for the strip I limited by the nodal lines i and j involves the following nodal parameters: u im and u jm for u as well as v im and v jm for v . Thus, each of the two approximations of displacements will require a linear polynomial function with two coefficients E m and F m. They are found from the boundary conditions, which for instance for u read: x = 0 ⇒
f m (x ) = u im
x = b ⇒
f m (x ) = u jm
This leads to the following matrix relation r u I I C pm w pm v = m =1
∑
where the matrix of coefficients is
x m π y x m π y − 1 sin 0 sin 0 b L b L I C pm = m π y x m π y x 0 0 cos 1 − cos L b L b and the nodal parameters for the single strip I are assembled into the vector
u im v im I w pm = u jm jm Again, similarly as in the plate analysis, to derive the strip stiffness matrix and the loading vector the energy approach is used. The strain energy for a strip I in the plane stress element can be expressed as L b
I U s
h = 20
∫∫
σ
L b
h dxdy = 20
I T I ε
0
∫ ∫ (σ ε
x x
+ σ y ε y + τ xy γ xy
)dxdy
0
The potential energy of the loading is L b
I U q
=
∫ ∫ (p (x )u + p (y )v )dxdy 0 0
The vector of strains can be expressed as r
ε
=
∑B
I I pm w pm
m =1
with the matrix of derivatives
1 − sin k m y b I B pm = 0 x k m 1 − cos k m y b
x x − k m 1 − sin k m y − k m sin k m y 0 b b 1 x 1 cos k m y − cos k m y k m cos k m y b b b 1 sin k m y b
0
0
m π . L With this in hand the physical law can be expressed in the matrix form and k m =
r
σ
=
∑D
I I I p B pm w pm
m =1
Thus, the energy parts can be put as: I
h r I T I T I I I w pm ∫ ∫ B pm D p B pm dxdy w pm ∑ 2 m =1 0 0 L b
U s =
I U q
r
=−
∑ m =1
L b I p (x ) C dxdy pm 0 0 p (y )
I T w pm
∫∫
Considering the fact, that the energy for the entire plane stress element consists of the energies for all the strips N
U e =
∑ U
I
I =1
and using the similar conditions of equilibrium as in the plate case
∂U e ∂w pm
=
0
we get the set of equilibrium equations for the entire element for each m -th harmonic function N
∑
I I K pm w pm
N
=
I =1
∑P
I pm
I =1
The stiffness matrix for the single strip I is obtained from L b
I K pm
=
∫∫
h
I T
I
I
B pm D p B pm dxdy
0 0
and its explicit form is
k p 1 k p 3 k p 4 k p 2 k p 5 I K pm = k p 1 sym .
− k p 5
k p 6 − k p 3
k p 2
where: 2
L Lbk m k p 1 = α + δ , 2b 6
2
Lbk m L k p 2 = β + δ , 6 2b
2
k p 4 = −
L Lbk m α + δ , 2b 6
k p 5 =
Lk m (γ + δ ) , 4
k p 3 =
Lk m (γ − δ ) 4 2
k p 6 =
Lbk m L β − δ 12 2b
and
α =
hE y hE x , β = , γ = h ν x α = h ν y β , δ = hE xy 1 − ν x ν y 1 − ν x ν y
The following steps of the solution are analogous to the plate analysis. 4. Analysis of compound structures Let us consider a box girder structure and its discretisation into the finite strips. The cross-section looks like this
Such structures can be considered as compound of rectangular elements, which are simultaneously subjected to bending and plane stress action. The stiffness matrix of such an element (strip) can be obtained as an appropriate assembly of stiffness matrices for strips in the bending state and in the plane stress state. The equilibrium for a single strip in its local co-ordinates can be expressed by the matrix equation ~ I ~ I ~ I K m w m = Pm Let us represent the strip matrices in bending and plane stress using the (2 ×2)-submatrices
k b 1 k b 3 k b 4 k b 2 k b 5 ~ I K bm = k b 1 sym .
− k b 5
k p 1 k p 3 k p 4 k p 2 k p 5 ~ I K pm = k p 1 sym .
− k p 5
k b 6 K bmii K bmij = I I − k b 3 K bmji K bmjj k b 2 I
I
k p 6 K pmii I K pmij I = I I − k p 3 K pmji K pmjj k p 2
Now the stiffness matrix for a strip in a compound structure can be given in the following form I K pmii I 0 K pmij 0 I I K bmii 0 K bmij ~ I K m = I K pmjj 0 I sym . K bmjj
This matrix corresponds to the following vectors of nodal displacement parameters (amplitudes) and nodal forces for the I -th strip and the m -th harmonic function ~ X ~im Y im ~ Z ~im ~ I M Pm = ~ im jm X ~ Y jm ~ w im ~ M jm
~ u im ~ v im ~ w im ~ ~ I = ϕ im w m ~ u ~ jm v jm w ~ im ~ ϕ jm
Each nodal line has four nodal parameters – three for displacements and one for transverse slope. ~ y, y
Displacements and forces in local co-ordinates I
i
~ z z
~ , ϕ jm ~ M jm
α ~ z
x ~ z z
~ x
~ y, y
ϕ jm , M jm
x ~ x
z
j
, jm v jm Y
~ ~ , X u jm jm ~ ~ , Z w jm jm
I
i j
~ ~ Y v jm , jm
~ y, y
Displacements and forces in global co-ordinates
u jm , X jm w jm , Z jm
~ x
x
The relations between the forces are ~ X jm = X jm cos α + Z jm sin α ~ Y jm = Y jm ~ Z jm = − X jm sin α + Z jm cos α ~ M jm = M jm what leads to the following transformation rule
cos α 0 − sin α 0
~ X ~ jm 1 Y jm =~ 0 cos α 0 Z jm Z ~ jm 0 0 1 M jm M jm 0
sin α 0 X jm 0 0 Y jm
The direction cosine matrix can be denoted as
cos α 0 C= − sin α 0
0
sin α 0
1
0
0
0 cos α 0 0
0
1
and the transformation matrix is
C 0 T= 0 C Now the transformation of vectors of the displacement and force nodal parameters can be given as ~ I I ~ I = Tw I w Pm = TPm m m ~ w m = TT w m I
I
~ I I Pm = TT Pm
and for the strip stiffness matrix we have ~ I I K m = TT K m T The following steps of the solution are analogous to the plate analysis. 5. Analysis of multi-span and column supported plates by the FSM The plates with inter-span supports, both point type or knife-edge type, can be solved using a combination of the FSM and the flexibility method. Reactions in the “additional” inter-span supports are considered as redundant forces and the single span plate for which the solution can be found using the pure FSM is taken as a modified (determinate) system. Thus the scheme of the flexibility method is set. The continuous reaction in the knife-edge support can be approximated as a set of point loads, for instance at the points where the support line intersects the nodal lines of the FSM discretisation.
P X 1 X 3 X N X 2 The identity of the modified system with removed supports and the original system with supports present is ensured by the kinematical conditions requiring the deflections at the removed supports to be zero
δ 1 = 0 δ 2
=
0
...
δ N = 0 After consideration of the basic states: X 1 = 1, X 2 = 1, X N = 1 and P and the application of the superposition rule the canonical equations of the flexibility method are obtained
δ 11 X 1 + δ 12 X 2
+ ... + δ 1N X N + δ 1P =
δ 21 X 1 + δ 22 X 2
+ ... + δ 2N X N + δ 2P =
0 0
...
δ N 1 X 1 + δ N 2 X 2
+ ... + δ NN X N + δ NP = 0
The flexibility coefficients δ ik , which are the appropriate negative deflections (the redundant forces vectors have upwards orientation, opposite to the axis z ) at the points where forces X i are applied are obtained from the FSM analysis of the modified one-span plate loaded by an appropriate loading state X k = 1. For instance, the set of displacements δ 1i , δ 2i ,… δ Ni follows from the state X i = 1.
δ 1i δ 3i δ ii
δ Ni
δ 2i 1 X i = Thus, N basic states and the state P must be solved to formulate the set of canonical equations. These equations are then solved and the values of redundant forces are found. Then the last stage of the solution follows, where the modified one-span plate is loaded simultaneously with the external loading and all the redundant forces. Alternatively, the superposition rule can be used but this requires the computer storage of the complete results for all the basic states, including deflections and moments.
