The precise finite strip method

PERGAMON

Computers and Structures 69 (1998) 773±783

The precise ®nite strip method W. X. Zhong a, Y. K. Cheung b, Y. Li b, * a

Research Institute of Engineering Mechanics, Dalian University of Technology, Dalian, People's Republic of China b Department of Civil and Structural Engineering of the University of Hong Kong, Hong Kong Received 20 December 1996; accepted 27 June 1997

Abstract This paper introduces a new development for the ®nite strip method. The precise integration method along the space coordinate is combined with the semi-analytical analysis of prismatic domain structures and then a novel precise ®nite strip method is proposed. By using such a novel semi-analytical method, the structure, which is as usual discretized along the transverse direction in a prismatic domain, will end up as a two point boundary value problem in the longitudinal direction. For the derived ODEs with constant coecient of the two point boundary value problems, the displacement function is computed by integration instead of being given by traditional series. The numerical result of the present method will be very precise, almost up to the computer precision. Thus, the precise ®nite strip method can deal with various local e€ect problems which cannot be solved by many traditional methods. # 1998 Published by Elsevier Science Ltd. All rights reserved. Keywords: Precise integration; Finite strip method; Two point boundary value problems

1. Introduction Y.K. Cheung [1] proposed the ®nite strip method for structural analysis of prismatic domain problems and in contrast to discretizing all the domain as commonly carried out in FEM, only the transverse crosssection is discretized. Along the longitudinal coordinate, however, the functions and their di€erential are still continuous and smooth. Thus, the method is regarded as a semi-analytical method [2]. Let y be assigned as the longitudinal coordinate and q(y) denotes the n-dimensional displacement vector at a nodal line to be solved, the ®nite strip di€erential equation can be derived as: _ 0, K22 q ‡ …K21 ÿ KT21 †_q ÿ K11 q ˆ ÿr0 ‡ m

…1†

where the dot (
) denotes a di€erential with respect to y and K12 ˆ KT21 , K22 K11 are n  n symmetric sti€ness matrices, which are y-independent; r0(y) and m0(y) are n-dimensional known vectors and correspond to the

* Corresponding author.

external driving forces, i.e. distribution force and distribution moment. The corresponding boundary conditions for Eq. (1) are given at the two ends of y0 and yf. For the domain y0
0045-7949/98/$ - see front matter # 1998 Published by Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 7 ) 0 0 1 0 5 - 3

774

W.X. Zhong et al. / Computers and Structures 69 (1998) 773±783

mal control [7, 8]. For constant coecient ODEs, the result of the present method may be very precise, almost up to the precision for computer expression. Thus, the novel method can deal with various local e€ects and can be e€ectively applied to the ODEs derived from other semi-analytical methods. The precise ®nite strip method also applies the method of mixed energy in the combination and elimination of intervals in the strip so that the numerical ill-conditioning problem is avoided when the intervals of integration are subdivided very small. Finally, the mixed energy of the whole strip can be transformed easily to the potential energy form and then the sti€ness matrix of the two end displacement vectors of the whole strip is obtained. Hence, it can be applied to a prismatic structure with complicated boundary conditions at the two ends. Based on the theory of multilevel substructural analysis, the strip can be regarded as a substructural module in a general purposed structural analysis system. Thus, the precise ®nite strip method can be merged into any multi-level substructure analysis system to extend its application area.

2. The Hamiltonian form of the ®nite strip equation The ®nite strip method is devoted to the analysis of prismatic structures with constant cross-sections. After discretization of the transverse cross-section, an ODE set along the longitudinal coordinate y can be derived. The displacement method is usually used in the discretization and the respective variational principle is derived as: … yf Y Y …q† ˆ L…q, q_ † dy, minimize , …2† y0

1 1 L…q, q_ † ˆ q_ T K22 q_ ‡ q_ T K21 q ‡ qT K11 q ÿ qT r0 ÿ q_ T m0 , 2 2 …3† where y = y0 and yf are the two ends of the strip and L is the Lagrangian function. Only the function q(y) is under variation. The above derivation is a purely displacement approach.

It is probably useful at this section to give the general form of Eqs. (1)±(3) for a beam-type strip. For the beam dynamic sti€ness matrix analysis based on the Timoshenko theory, the longitudinal coordinate is selected to be the strip axis and the nodal displacement vector q is composed of the transverse displacement W, the longitudinal displacement V and the rotation of the cross-section y, as shown in Fig. 1. Let the displacement function be written as q = {W, V, y}T, then the Lagrangian is: L…q, q_ † ˆ

1 _2 1 1 2 _ ‡ y†2 EJy ‡ EbtV_ ‡ mGbt…W 2 2 2 1 ÿ ro 2 ‰btW 2 ‡ btV 2 ‡ Jy2 Š, 2

in which t and b are the thickness and width of the strip, respectively. Putting F = bt, it can be seen that: 2 3 mGF 0 0 6 7 K22 ˆ4 0 EF 0 5, 0 0 EJ 2 3 2 0 0 ÿro F 6 7 K11 ˆ4 0 ÿro 2 F 0 5, 2

0

0 0 6 K21 ˆ4 0 0 0 0

0 3

mGF ÿ ro 2 J

mGF 7 0 5, 0

and the external force vectors r0=m0=0. The Eq. (1)type ODE cannot be solved directly in Lagrangian system because there are three terms of q, q_ and q_ which make it impossible to use the method of separation of variables. Applying the analogy theory in optimal control, the equations can be transformed into the Hamiltonian system. Let us take the Legendre's transformation, but corresponding to the station space of coordinate y instead of time space of t used in the Lagrangian system. Let us de®ne: p ˆ ÿ@ L=@ q_ ˆ ÿ…K22 q_ ‡ K21 q† ‡ m0 ,

…4†

where the minus sign is only for convenience and p can be interpreted as the internal force. A set of the new dual equations can be derived in the following

Fig. 1. A typical strip.

W.X. Zhong et al. / Computers and Structures 69 (1998) 773±783

form: q_ ˆ Aq ÿ Dp ÿ m,

…5a†

p_ ˆ ÿBq ÿ AT p ‡ r,

…5b†

in which ÿ1 ÿ1 D ˆKÿ1 22 , A ˆ ÿK22 K21 , B ˆ K11 ÿ K12 K22 K21 ,

r ˆr0 ‡ AT m0 , m ˆ ÿDm0 :

…6†

Evidently, B and D are symmetric matrices and the matrices A, B and D are termed the system matrices. The respective boundary condition for Eqs. (5a) and (b) is the same at the two ends y0 and yf. Thus, the original problem with n second-order ODEs shown in Eq. (1) has been simpli®ed into the 2n ®rst-order ODEs and the integration solution can then be carried out directly with the feature of such a two point boundary value problem. Citing the Timoshenko beam dynamic sti€ness matrix as an example, the corresponding system matrices A, B and D are: 2 3 2 3 0 0 1 ÿro 2 F 0 0 6 7 6 7 A ˆ4 0 0 0 5 , B ˆ 4 0 ÿro 2 F 0 5, 0 0 0 0 0 ÿro 2 J 2 3 ÿ1 …mGF † 0 0 6 7 ÿ1 D ˆ4 0 …EF † 0 5: 0

0

…EF †ÿ1

The corresponding variational Eqs. (5a)±(b) is given as … yf d …pT q_ ÿ H…p, q†† dy ˆ 0,

principle

y0

for …7†

where the Hamilton function is 1 1 H…p, q† ˆ pT Aq ‡ qT Bq ÿ pT Dp ÿ qT r ÿ pT m: 2 2

…8†

It is easily veri®ed that: @ H=@ q ˆ Aq ÿ Dp ÿ m ˆ q_ @ H=@ q ˆ AT p ‡ Bq ‡ r ˆ ÿ_p: Thus, the new dual Eqs. (5a)±(b) can also be written as: q_ ˆ @ H=@ p, p_ ˆ ÿ@ H=@ q:

…9†

As is well known, the above dual di€erential equation is not convenient for numerical computation, so the equation for a ®nite length Dy must be derived instead. Although the ®nite di€erence method can be used to do step by step integration, it is, however, not suciently precise, so instead the precise integration method [6±9] is applied here.

775

3. The solution of the dual equations and the interval mixed energy matrix function Let any two stations ya and yb>ya be selected arbitrarily within the domain (y0, yf) of the strip and let them make up an interval (ya, yb). Evidently, if the vectors qa and pb are given at the two ends of the interval respectively, the solution q and p in the interval (ya, yb) is ®xed, where qa means the displacement vector at station ya and pb is the force vector at station yb. The mixed energy of the interval (ya, yb) can be de®ned as [8]: … yb …10† V…qa ,pb † ˆ pTb qb ÿ …pT q_ ÿ H…p, q†† dy: ya

The importance of the above equation is that by applying Eq. (9), its variational equation can be derived as: dV…qa , pb † ˆ qTb dpb ‡ pTa dqa i:e: @ V=@ qa ˆ pa , @ V=@ pb ˆ qb :

…11†

From Eq. (10), the mixed energy V is obviously a quadratic functional of qa and pb. The general form of a quadratic form is: 1 V…qa , pb † ˆ pTb Fqa ‡ qTa Qqa 2 1 T ÿ pb Gpb ‡ hTq qa ‡ hTp pa ‡ Const, 2

…12†

where the Const term has no e€ect on variation and will be neglected later. In Eq. (12), the matrices Q, G, F are functions of ya and yb and are of dimension n  n. Q and G are symmetric matrices and for the static analysis problem, G is usually positive de®nite. Substituting Eq. (12) into Eq. (11) yields:        hp qb F ÿG qa ˆ ‡ …13a; b† hq pa pb q FT where the nonhomogeneous terms hp and hq are from external forces. Evidently the physical interpretation for Eqs. (13a) and (b) is that G is the ¯exibility matrix at end b with end a clamped; Q is the sti€ness matrix at end a with end b free; and F is the transfer of displacement from end a to end b. The model is always clamped at end a, but free at end b, which is the `cantilever beam' type boundary condition corresponding to the mixed energy formulation. Eqs. (13a) and (b) is derived from the dual Eqs. (5a) and (b) or Eq. (9) and is another form of Eqs. (5a) and (b). The mixed energy matrix functions Q, G and F are certainly dependent on the matrices A, B and D and the vectors hp and hq also depend on r and m. Their interrelationship can be expressed by di€erential equations, which are derived as follows.

776

W.X. Zhong et al. / Computers and Structures 69 (1998) 773±783

Let yb be increased to yb+Dyb with Dyb in®nitesimally small, while ya and qa remained ®xed and let pb be increased according to the di€erential Eq. (5b) at yb, i.e. Dpb=_pb Dyb. The result is that vector pa must not be changed, i.e. Dpa=0 and the vector qb must be increased to qb+Dqb, with Dqb=_qbDyb. Based on that observation, it is possible to derive the following di€erential equations from Eqs. (13a) and (b): @ hq @Q @ FT qa ‡ p ‡ FT p_ b ‡ , @ yb @ yb b @ yb @ hp @F @G q_ b ˆ q ÿ p ÿ G_pb ‡ : @ yb a @ yb b @ yb 0ˆ

Substituting Eq. (4) into the above equation to eliminate p_ b and q_ a and noting that the three vectors qa, pb and qb are not mutually independent, it is also possible to eliminate qb further by using Eqs. (13a) and (b), thus yielding     @ hp @F @G ÿ …A ÿ GB†F Qa ‡ ÿ ‡ D ‡ AG ‡ GAT ÿ GBG pb ‡ ÿ …A ÿ GB†hp ÿ Grb ‡ mb ˆ 0 @ yb @ yb @ yb 

  T  @ hq @q @F ÿ FT BF qa ‡ ÿ FT …AT ÿ BG†F pb ÿ FT …Bhb ÿ rb † ‡ ˆ 0: @ yb @ yb @ yb

Because the two end vectors qa and pb in the above equations are arbitrary, hence: @F @G @Q ˆ …A ÿ GB†F, ˆ D ‡ AG ‡ GAT ÿ GBG, ˆ FT BF, @ yb @ yb @ yb @ hp @ hq ˆ …A ÿ GB†hp ‡ Grb ÿ mb , ˆ FT …Bhp ÿ rb †: @ yb @ yb

…14† …15a; b†

Repeating the same procedure with yb, pb ®xed and letting ya, pa vary, it is possible to arrive at: @F @G @Q ˆ ÿF…A ÿ DQ†, ˆ ÿFDFT , ˆ ÿB ÿ AT Q ÿ QA ‡ QDQ, @ ya @ ya @ ya @ hp @ hq ˆ F…Dhq ‡ ma †, ˆ ra ÿ …AT ÿ QD†hq : @ ya @ ya

…16† …17a; b†

The boundary conditions for these ®rst-order di€erential equations are: Q ˆ G ˆ 0, F ˆ In , hq ˆ hp ˆ 0 when ya ÿ4yb ,

…18†

where In is the unit matrix of dimension n. For the normal case, A, B and D are independent of y because of the uniform cross-section; therefore, Q, G, F, hq and hp are dependent only on the interval length t. Let: t ˆ yb ÿ ya

…19†

then Q…za , zb † ˆ Q…t†, @ Q=@ yb ˆ dQ=dt, @ Q=@ ya ˆ ÿdQ=dt, etc:

…20†

and Eqs. (14), (15a) and (15b) are then reduced to dq=dt ˆ FT BF, dG=dt ˆ D ‡ AG ‡ GAT ÿ GBG, dF=dt ˆ …A ÿ GB†F,

…21†

dhp =dt ˆ …A ÿ GB†hp ‡ Gr ÿ m, dhq =dt ˆ FT …Bhp ÿ r†:

…22†

The above di€erential equations give the interrelationship between the system matrices and the interval matrices Q, G and F. However, the same reduction from Eqs. (16), (17a) and (17b) gives again: dQ=dt ˆ B ‡ AT Q ‡ QA ÿ QDQ, dG=dt ˆ FDFT , dF=dt ˆ F…A ÿ DQ†,

…23†

dhp =dt ˆ ÿF…Dhq ‡ m†, dhq =dt ˆ …AT ÿ QD†hq ÿ r:

…24†

W.X. Zhong et al. / Computers and Structures 69 (1998) 773±783

777

The respective boundary conditions are [see Eq. (18)]: Q ˆ G ˆ 0, F ˆ In , hq ˆ hp ˆ 0 when t ˆ 0, Eqs. (21) and (22) appear to be quite di€erent from Eqs. (23) and (24), although they represent the derivatives (forward and backward) of Q, G, F and hq and hp. Having derived the di€erential equations and boundary conditions for the mixed energy matrix functions, it is still a challenging problem to obtain a solution because of nonlinearity. However, when t is a very small interval, the Taylor series expansion method of the solution can be applied. Let: Q…t† ˆ y1 t ‡ y2 t2 ‡ y3 t3 ‡ y4 t4 ‡ O…t5 †,

…25a†

G…t† ˆ g1 t ‡ g2 t2 ‡ g3 t3 ‡ g4 t4 ‡ O…t5 †,

…25b†

F 0 …t† ˆ f1 t ‡ f2 t2 ‡ f3 t3 ‡ f4 t4 ‡ O…t5 †, F…t† ˆ In ‡ F 0 …t†,

…25c†

hp …t† ˆ p1 t ‡ p2 t2 ‡ p3 t3 ‡ p4 t4 ‡ O…t5 †,

…25d†

hq …t† ˆ c1 t ‡ c2 t2 ‡ c3 t3 ‡ c4 t4 ‡ O…t5 †,

…25e†

where yi, gi and fi (i = 1, 2, 3, 4) are n  n coecient matrices and pj and cj are n-dimensional vectors to be determined. Substituting Eq. (25b) into Eq. (21) and comparing the coecients of di€erent powers of t, gives: g1 ˆ D, g2 ˆ …Ag1 ‡ g1 AT †=2, g3 ˆ …Ag2 ‡ g2 AT ÿ g1 Bg1 †=3, g4 ˆ …Ag3 ‡ g3 AT ÿ g2 Bg1 ÿ g1 Bg2 †=4:

…26a†

Applying similar procedures to the other equations in Eqs. (25a)±(e), the coecient can be obtained successively: f1 ˆ A, f2 ˆ …Af1 ÿ g1 B†=2, f3 ˆ …Af2 ÿ g2 B ÿ g1 Bf1 †=3, f4 ˆ …Af3 ÿ g3 B ÿ g2 Bf1 ÿ g1 Bf2 †=4,

…26b†

for F, and then y1 ˆ B, y2 ˆ …fT1 B ‡ Bf1 †=2, y3 ˆ …fT2 B ‡ Bf2 ‡ fT1 Bf1 †=3, y4 ˆ …fT3 B ‡ Bf3 ‡ fT2 Bf1 ‡ fT1 Bf2 †=4,

…26c†

for Q. After the coecient matrices have been calculated, the coecient vectors for hq and hp are given as: c1 ˆ ÿ r, c2 ˆ …fT1 c1 ‡ Bp1 †=2, c3 ˆ …fT2 c1 ‡ Bp2 ‡ fT1 Bp1 †=3, c4 ˆ …fT3 c1 ‡ Bp3 ‡ fT1 Bp2 ‡ fT2 Bp1 †=4, p1 ˆ ÿm, p2 ˆ …f1 p1 ÿ Dc1 †=2, p3 ˆ …f2 p1 ÿ Dc2 ÿ f1 Dc1 †=3, p4 ˆ …f2 p1 ÿ Dc3 ÿ f1 Dc2 ÿ f2 Dc1 †=4:

…27†

Taylor series expansion up to the fourth-order means, when t is small enough, Q, G, F, hq and hp can be calculated very precisely. Such properties enable us to use the 2N type algorithm, i.e. for N = 20, a given length of interval will be subdivided into 220=1 048 576 very small intervals of equal length, which implies t being divided by 106 and the neglected term in the Taylor series expansion is t with a multiplier of 10ÿ24 compared with the ®rst term. Such precision is in fact beyond the computer double precision of 10ÿ16 for real numbers. The computational techniques are given below. 4. Computational techniques for precise integration 4.1. Interval combination and the submatrix Eqs. (13a) and (b) has given the relationship of the state vectors at the two ends of ya and yb. For two given adjacent intervals (ya, yb) and (yb, yc) connected at yb as shown in Fig. 2, we can eliminate the interior state vector de®ned at yb so as to form a new interval and the corresponding dual equations will then connect the state vectors de®ned at ya and yc, respectively. Applying Eqs. (13a) and (b) to the intervals 1 and 2 gives: qb ˆ F1 qa ÿ G1 pb ‡ hp1

…28a†

778

W.X. Zhong et al. / Computers and Structures 69 (1998) 773±783

Fig. 2. Combining intervals 1 and 2 to form a new interval a.

pa ˆ Q1 qa ‡ FT1 pb ‡ hq1 ,

…28b†

qc ˆ F2 qb ÿ G2 pc ‡ hp2

…29a†

pb ˆ Q2 qb ‡ FT2 pc ‡ hq2 ,

…29b†

based on which the equation of the combined interval (ya, yc) can be derived and it should still have the form as follows qc ˆ Fc qa ÿ Gc pc ‡ hpc

…30a†

pa ˆ Qc qa ‡ FTc pc ‡ hqc ,

…30b†

where the subscript c represents the corresponding quantities. Solving the vectors qb and pb from Eqs. (28a) and (b) yields qb ˆ …I ‡ Q2 G1 †ÿT ‰…F1 qa ‡ hp1 † ÿ G1 …FT2 pc ‡ hq2 †Š …31a†

tures, only one needs to be analyzed for its sti€ness and the data can be transferred to all other identical substructures. Let us divide the strip along the longitudinal direction by using the 2N-type algorithm. First, assign the length of the integration domain to be ysRyfÿy0, which can be further divided uniformly into 2N very small subintervals of length Dy. With N = 20, the length of a subinterval will be Dy = ys/1 048 576 and all these subintervals are identical and periodic. Therefore, each step integration will reduce the number of subintervals by half. After N = 20 steps, all 1 048 576 subintervals would be integrated into an interval of length ys as required. It is further suggested to subdivide the whole strip domain with length L = yfÿy0 into 2m integration lengths ys, i.e. ys=L/2m, so that in order to obtain the solution for the whole strip one can also apply the 2Ntype algorithm (with N = m) to combine the 2m integration lengths ys. Using the above division method, each interval has identical periodicity and therefore the mixed energy submatrices functions are of the same form, i.e.: Qc ˆ Q ‡ FT …Qÿ1 ‡ G†ÿ1 F Gc ˆ G ‡ F…Gÿ1 ‡ Q†ÿ1 FT Fc ˆ F…I ‡ QG†ÿT F

pb ˆ …I ‡ Q2 G1 †ÿ1 ‰…FT2 pc ‡ hq2 † ‡ Q2 …F1 qa ‡ hp1 †Š: …31b† Finally substituting Eqs. (31a) and (b) back to Eqs. (28a), (28b), (29a) and (29b), and comparing with Eqs. (30a) and (b), will ®nally result in the desired interval combination equations: ÿ1 Qc ˆ Q1 ‡ FT1 …Qÿ1 2 ‡ G1 † F 1 ,

…32a†

ÿ1 T Gc ˆ G2 ‡ F2 …Gÿ1 1 ‡ Q2 † F 2 ,

…32b†

Fc ˆ F2 …I ‡ Q2 G1 †ÿT F1 ,

…32c†

hqc ˆ hq1 ‡ FT1 …I ‡ Q2 G1 †ÿ1 …hq2 ‡ Q2 hp1 †,

…33a†

hpc ˆ hp2 ‡ F2 …I ‡ Q2 G1 †ÿT …hp1 ÿ G1 hq2 †:

…33b†

Note that due to Eqs. (30a) and (b), Qc, Gc and Fc are the mixed energy submatrices combining from station a to c. 4.2. The precise division and integration for a given interval in the strip In structural mechanics, the substructure technique had been widely used to improve computational eciency. For a structure with identical periodic substruc-

…34†

and hqc ˆ hq ‡ FT …I ‡ QG†ÿ1 …hq ‡ Qhp † hpc ˆ hp ‡ F…I ‡ QG†ÿT …hp ÿ Qhq †:

…35†

For the above equations, the basic calculation of this 2N type algorithm is the repeated execution N times and each time, the computed Qc, Gc, Fc, hpc and hqc are put into the right-hand side of the above equations to calculate the next round's Qc, Gc, Fc, hpc and hqc of the doubled length subinterval. The advantage of this division and integration technique based on interval combination is higher precision and less calculation steps compared with the step by step integration commonly used. Note that all the derivations above are exact except for the Taylor series expansion. This error can be controlled by using very small Dy, which is why the integration domain ys is further divided into 2N smaller intervals. However, the use of 2N type algorithm changes the order of integration of the original di€erential Eqs. (23) and (24), which is usually considered step by step forward. As shown in Fig. 3, there are three subintervals numbered 1±3 which will be combined to generate a new subinterval c. Obviously, it can be performed in two ways:

W.X. Zhong et al. / Computers and Structures 69 (1998) 773±783

779

4. Sum the matrices for the interval ys, i.e.: Q… ys † ˆQc ; G… ys † ˆ Gc ; F 0 … ys † ˆIn ‡ Fc 0 ; hp … ys † ˆ hpc ; hq … ys † ˆ hqc ;

…36†

Fig. 3. Two procedures for interval combination: (1) 1 + 2 4 a, a + 3 4 C; (2) 2 + 3 4 b, b + 1 4 C.

1. combine subintervals 1 and 2 to get a new subinterval a ®rst, then combine a and 3 to get the ®nal subinterval c; 2. combine subintervals 2 and 3 ®rst to get the subinterval b then combine 1 and b to get the ®nal subinterval c. Based on the matrix inversion lemma [10, 11] and Eqs. (23) and (24), it is not dicult to prove that the results of the two combinations are identical. Therefore, the use of 2N algorithm is justi®ed. Care must be taken in the computation process. Since the matrix F is the sum of an unit matrix and a very small matrix F 0 , the addition should not be executed when the interval is very small so as to avoid unnecessary numerical error. Only F 0 (t) is generated and stored in the computer memory instead of F(t) and Eq. (34) should be replaced by: Qc ˆ Q ‡ …In ‡ F 0 †T …Qÿ1 ‡ G†ÿ1 …In ‡ F 0 †,

…34 0 a†

Gc ˆ G ‡ …In ‡ F 0 †…Gÿ1 ‡ Q†ÿ1 …In ‡ F 0 †T ,

…34 0 b†

Fc 0 ˆ …F 0 ÿ GQ=2†…I ‡ QG†ÿ1 ‡ …I ‡ QG†ÿ1 …F 0 ÿ GQ=2† ‡ F 0 …I ‡ QG†ÿ1 F 0 :

…34 0 c†

It is easily seen from Eqs. (34 0 a)±(b) that originally Q and G are non-negative symmetric matrices and the additional terms are also symmetric and non-negative de®nite, so that the newly formed Qc and Gc must be no less than the original Q and G, which dictates that the combination based on Eqs. (34 0 a)±(c) will give ever increasing matrices Q and G. The following programming process gives the computational steps of the precise integration for the submatrix functions of the mixed energy in a given length ys interval. 1. t = ys/2N. 2. Generate Q(t), G(t), F 0 (t), hp and hq from Eqs. (25a)±(e), (26a) and (b) and (27). 3. Integrate the mixed energy submatrices N times to obtain the submatrices of the interval ys by the execution of Eqs. (34 0 a)±(c) and (35), i.e.: 0

0

‰Q ˆ Qc ; G ˆ Gc ; F ˆ Fc ; hp ˆ hpc and hq ˆ hqc ;Š

4.3. The potential energy formulation of an interval by the matrix functions Since an interval domain can also be treated as a substructure, potential energy formulation for intervals should also be possible. The reason why the above mentioned methods are seldom used in FEM. is because that the length of the subinterval t for the 2N algorithm is very small and the submatrices will tend to become in®nitesimal when t approaches zero. Hence, in practice serious numerical ill-conditioning will prevent the widespread application of the displacement method using the 2N algorithm. On the other hand, the mixed energy sub-matrices Q, G and F 0 can be seen to be of the order of t when t approaches zero. Thus, the numerical ill-conditioning is nearly completely eliminated and superior performance can be expected. However, for a ®nite length interval, the displacement representation will have no such ill-conditioned problem as mentioned above, so that after having compute the interval matrices Q(ys), G(ys), F(ys) and vectors hp(ys) and hq(ys) for the integration domain, the representation can be transformed back to the displacement form. The interval sti€ness matrix can be obtained by transforming Eqs. (13a) and (b) with displacement vectors at the right-hand side: pa ˆ …Q ‡ FT Gÿ1 F†qa ÿ FT Gÿ1 qb ‡ FT Gÿ1 hp ‡ hq , …37a† pb ˆ Gÿ1 Fqa ÿ Gÿ1 qb ‡ Gÿ1 hp

…37b†

Comparing with the standard form of displacement representation: pa ˆ Kaa qa ‡ Kab qb ÿ fa ,

…38a†

pb ˆ ÿKba qa ÿ Kbb qb ‡ fb ,

…38b†

it can be seen that Kbb ˆ Gÿ1 , Kba ˆ ÿKbb F, Kaa ˆ Q ‡ FT Kbb F

…39†

and fa ˆ ÿ…FT Kbb hp ‡ hq †, fb ˆ Kbb hp :

…40†

780

W.X. Zhong et al. / Computers and Structures 69 (1998) 773±783

The corresponding potential energy is: 1 T 1 q Kaa qa ‡ qTb Kbb qb ‡ qTa Kab qb ÿ qTa fa ÿ qTb fb , 2 a 2

…41†

where fa, fb are the equivalent external force vectors at the stations ya and yb. If ya and yb are selected as y0 and yf of the strip, then the sti€ness and equivalent forces of the whole strip are obtained. The solution will give the two end displacement vectors qa and qb (i.e. q0 and qf) and then they are substituted into Eqs. (38a) and (b) to ®nd pa and pb so that the two end state vectors are obtained. 4.4. State vector solving After the two end state vectors have been solved from the whole structure, the control will be transferred back to the precise strip analysis program module with the two computed end state vectors. The internal state vectors can be solved by reversing the interval combination process. First of all the middle station between y0 and yf is computed from Eqs. (31a) and (b), then the 1/4 and 3/4 stations, then 1/8, 3/8, 5/ 8, 7/8 stations, etc.

5. Numerical examples

Fig. 4. Rectangular plates with simply supported and clamped boundaries: (a) a uniformly distributed load p; (b) a central concentrated load P.

Some numerical examples will be studied to demonstrate the applications in structural analysis and the features of the present method. The numerical results

Table 1 Accuracy test on simply support plate with uniformly distributed load The method of ®nite strip

W1(max)

Mx1 My1

Mxy4

8  8 strip by Shi$ 5  4 strip by Cheung$ 8  8 spline strips$ 8 strip precise integration The exact solution Multiplier

0.406236 0.406 0.406236 0.4074 0.406236 qa4/D  10ÿ2

0.47950 0.481, 0.478 0.47951, 0.47950 0.4778 0.47886 qa210ÿ1

0.32498 0.32499 0.32032 0.32482 qa210ÿ1

$ Referenced from [12].

Table 2 Accuracy test on plate with four clamped sides and under a central concentrated load The method of ®nite strip

W1(max)

8  8 strip by Z.C. Shi$ HO2 strip by Cheung$ 8  8 spline strips$ 8 strip precise integration The exact solution Multiplier

0.560175 0.555 0.56027 0.565 0.56120 Pa2/D  10ÿ2

My2

P

0.124182 0.12654 0.12429 0.125936 0.12577

Mx3

P

0.124182 0.11276 0.12427 0.125936 0.12577

W.X. Zhong et al. / Computers and Structures 69 (1998) 773±783

Fig. 5. A plate with two simply supported sides and two free edges under a concentrated load P at point B.

obtained are compared with those obtained by other researchers. From these comparisons, the accuracy, convergence and eciency of the precise ®nite strip method are found to be very good.

5.1. Rectangular plates A typical application for the ®nite strip method is the analysis of plate problems. Rectangular plates with di€erent boundary conditions and loads are shown in Fig. 4. Due to symmetry conditions, a quarter of the plate is divided into eight strips. For the traditional ®nite strip method, the longitudinal trial function takes eight terms of the series or eight spline functions. For the precise ®nite strip method, the longitudinal function is computed from the precise integration in 2n steps (n = 4). The numerical results are compared with other solutions shown in Tables 1 and 2.

781

For the traditional ®nite strip method, the general trial functions are restricted to analyze the concentrated load e€ect exactly, but the precise integration can describe any local change in detail, and in Table 2 it can be seen that the precise integration strip method gives the most accurate moments under a concentrated load. Consider further a plate with two sides simply supported and under a concentrated load P at point B shown in Fig. 5. The numerical results are shown in Table 3 and Fig. 6. Compared with the solutions by Coull, FEM, the traditional ®nite strip methods, the computing by the precise ®nite strip method needs less time for convergence. 5.2. Architectural sandwich panels A typical sandwich panel section is divided into a number of layers which are clamped at one end and free at the other end. The core of high density has a Poisson's ratio of 0.25 and a Young's modulus of 5175 kN/m2 and the core of low density has a Poisson's ratio of 0.177 and a Young's modulus of 4685 kN/m2. Poisson's ratio and Young's modulus of each steel face are 0.3 and 207 GPa, respectively. The thicknesses of the face and core are 0.508 and 12.7 mm respectively, as shown in Fig. 7. The ¯exural behavior of sandwich panels subject to bending has been investigated analytically, experimentally and numerically [13] and the boundary e€ect cannot be described in more detail by any common method. Consider a sandwich panel with one end clamped and the other end free and subjected a uniformly distributed load of 4790 N/m2. By the precise

Table 3 Convergence test on the moment My(y = 0) of a plate with simple supports at two ends and under a concentrate load P at point B The node x = ÿ b/2 x = ÿ 7b/16 x = ÿ 6b/16 x = ÿ 5b/16 x = ÿ b/4 x = ÿ 3b/16 x = ÿ b/8 x = ÿ b/16 x=0 x = b/16 x = b/8 x = 3b/16 x = b/4 x = 5b/16 x = 6b/16 x = 7b/16 x = b/2 $ Referenced from [12].

Typical strip$ (16  35) 0.1391 0.1443 0.1509 0.1602 0.1725 0.1885 0.2092 0.2369 0.2759 0.3366 0.4865 0.3466 0.2962 0.2677 0.2513 0.2424 0.2389

Typical strip$ (16  99) 0.1391 0.1440 0.1509 0.1602 0.1725 0.1885 0.2092 0.2366 0.2747 0.3315 0.5286 0.3415 0.2950 0.2675 0.2512 0.2424 0.2389

Spline strip$ (16  8) 0.1396 0.1445 0.1515 0.1609 0.1734 0.1896 0.2108 0.2394 0.2812 0.3516 0.4763 0.3617 0.3016 0.2704 0.2503 0.2473 0.2400

Precise strip 0.14051 0.14405 0.15102 0.16046 0.17285 0.18890 0.20981 0.23785 0.27800 0.35479 0.48594 0.36482 0.29827 0.26873 0.25184 0.24268 0.23944

782

W.X. Zhong et al. / Computers and Structures 69 (1998) 773±783

Fig. 6. Moment Mx and shearing force Qz along central section.

®nite strip method, the boundary layer e€ect can be studied in great detail and as an illustration the numerical solutions of the shear forces are shown in Fig. 8.

6. Concluding remarks The solution of second-order matrix di€erential equations is of great importance in engineering appli-

Fig. 7. A sandwich panel with variable core densities.

W.X. Zhong et al. / Computers and Structures 69 (1998) 773±783

783

Fig. 8. Shear force variation of sandwich layers under a uniformly distributed load.

cations as the semi-analytical ®nite strip method and the wave propagation problems always end up as two point boundary value problems. For structural vibration, the precise ®nite strip method has been developed for initial value problems and can give the numerical solution very precisely since its error only comes from the arithmetic rounding o€ error of the computer. The method subdivides each interval into a very dense mesh along the longitudinal coordinate; thus, it is especially suited to calculate the boundary or local e€ects of structures. With the precise integration form, the novel method also has no restriction in various loading and complicated boundary conditions at two ends.

Acknowledgements The appointment of the senior author as Honorary Professor of the University of Hong Kong is gratefully acknowledgment.

References [1] Cheung YK. Finite strip method in structural analysis, 2nd edn. Pergamon Press, Oxford, 1976. [2] Zienkiewicz OC, Taylor RL. The Finite Element Method, Vol. 2, 4th edn. McGraw-Hill Press, U.K, 1991.

[3] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes. Cambridge University Press, U.K., 1992. [4] Cheung YK, Fan SC, Wu CQ. Spline ®nite strip method in structural analysis. In Proceedings of the International Conference on Finite Element Methods, Shanghai, 1982, pp. 36±46. [5] Kong JS. Analysis of plate-type structures by ®nite strip, ®nite prism and ®nite layer methods. Ph.D. thesis, The University of Hong Kong, Hong Kong, 1994. [6] Zhong WX, Williams FW. A precise time step integration method. Procedings of International Mechanics Engineering, Part C 1994;208:427±30. [7] Zhong W, Zhong X. Computational structural mechanics optimal control and semi-analytical method for PDE. Computers and Structures 1990;37:993±1004. [8] Zhong W. Computational Structural Mechanics and Optimal Control. Dalian University of Technology Press, China, 1993 (in Chinese). [9] Angel E, Bellman R. Dynamic Programming and Partial Di€erential Equations. Academic Press, New York, 1972. [10] Stengel RF. Stochastic Optimal Control. Wiley, New York, 1986. [11] Howson WP, Banerjee JR, Williams FW. Concise equations and program for exact eigen-solutions of plane frames including member shear. Advanced Engineering Software 1983;5:137±41. [12] Wu CQ, Cheung YK, Fan SC. Spline Finite Strip Method in Structural Analysis. Guangdong Academic Press, Cuangzhou PR. China, 1986 (in Chinese). [13] Boswell LF. (ed. Bull JW), Analysis of Thin-walled Box Beam Structures, Vol. 6. Elsevier, Oxford, 1990, pp. 188± 192.