INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 3689—3702 (1997) A CHEBYSHEV COLLOCATION MULTIDOMAIN METHOD TO SOLVE THE REISSNER—MINDLIN EQUATIONS FOR THE TRANSIENT RESPONSE OF AN ANISOTROPIC PLATE SUBJECTED TO IMPACT BO KJELLMERT* ¸ulea> ºniversity of ¹echnology, Division of Physics, S-97187 ¸ulea> , Sweden ABSTRACT The transient response of an anisotropic rectangular plate subjected to impact is described through a Chebyshev collocation multidomain discretization of the Reissner—Mindlin plate equations. The trapezoidal rule is used for time-integration. The spatial collocation derivative operators are represented by matrices, and the subdomains are patched by natural and essential conditions. At each time level the resulting governing matrix equation is reduced by two consecutive block Gaussian eliminations, so that an equation for the variables at the subdomain corners has to be solved. Back-substitution gives the variables at all other collocation points. The time history as represented by computed contour plots has been compared with analytical results and with photos produced by holographic interferometry. The agreements are satisfactory. ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) No. of Figures: 6. No. of Tables: 2. No. of References: 14. KEY WORDS: multidomain Chebyshev collocation; partitioning; transient response; Reissner—Mindlin plate 1. INTRODUCTION Many experimental and mathematical techniques have been developed to study the response of plates subjected to impact, and the problem has many different aspects. This paper examines only a few of the techniques which have inspired the present work. The transient displacement field of a plate may be visualized experimentally using holographic interferometry. A series of snapshots is produced which shows how the out-of-plane component of the field evolves from the beginning of the impact and forward.1 This optical method gives much information on the displacement field and may serve as a test-bench for plate models. The Reissner—Mindlin plate model is often adopted to analyse the three coupled waves which together form the displacement field. Recently, El-Rehab investigated the waves in an isotropic disk soon after impact using an approximate analytical method.2 To find analytic expressions for waves in an anisotropic plate one can rely on the classical Kirchhoff plate model and assume that the impact pressure is delta shaped in space and time. Transform methods then yield expressions for the displacement field of the plate, expressions which are valid before the waves reach the boundary of the plate.3 * Correspondence to: B. Kjellmert, Division of Physics, Luleas University of Technology, S-97187 Luleas , Sweden. E-mail: Bo.Kjellmert @ mt.luth.se CCC 0029—5981/97/203689—14$17.50 ( 1997 John Wiley & Sons, Ltd. Received 21 October 1995 Revised 19 December 1996 3690 B. KJELLMERT In general, only numerical methods are available for the analysis of the transient response of a plate. This paper describes a method which gives solutions of engineering accuracy for an anisotropic rectangular Reissner—Mindlin plate. Standard methods are combined as follows. The trapezoidal rule is used for time integration. At each time-level the boundary value problem is discretized by a multidomain Chebyshev collocation method. All subdomains are rectangular, have the same size and contain the same set of collocation points. In principle, the discretization, and the matrix expressions for the collocation derivatives4,5 result in a large matrix equation for all unknowns at all collocation points. This equation is solved in steps. Indeed, a domain partitioning method leads to two consecutive block Gaussian eliminations which transform the large equation into an equation for the variables at all subdomain corners. Here one has taken advantage of the fact that the governing differential equations and the discretizations are identical in each subdomain, so that many blocks of the large matrix equation also become identical. In all, three matrices have to be decomposed and some help matrices must be computed. Then a small number of matrix multiplications are sufficient to proceed from one time-level to another. The partitioning technique makes the solver efficient. Four facts have been relevant to the choice of method. (1) The anticipated solution field is smooth. (2) The initial-boundary value problem is stiff. (This property comes from the Reissner— Mindlin theory.6) (3) The plate is rectangular. (4) The impact pressure is concentrated to a small part of the plate. The first two points are important for all thin plates, while the last two take into account the specific experimental situation to be modelled. Point (1) indicates that a spectral discretization in space is more efficient than a finite difference or finite element discretization.7 Due to point (3) Legendre or Chebyshev polynomial expansions form expedient approximations.5 Points (2) and (4) imply that any discretization in space must be fine. The author has chosen a Chebyshev spectral method in space in spite of the difficulties encountered because of the concentrated impact. Being a collocation method the boundary conditions are enforced, both the external ones and the internal ones between subdomains. In other words, continuity and flux continuity are enforced. By physical arguments this strong subdomain coupling is important for an accurate modelling of waves. Here the present method has an advantage in comparison with finite element or spectral element methods, where flux continuity is ensured only as part of the convergence process. The plate and impact models are given in Section 2 and the discretization is presented in Section 3. Some details of the block Gaussian eliminations are also presented in Section 3. The out-of-plane deformation field as computed by the present method is compared in Section 4 with analytical and experimental results. The discussion is concentrated to one anisotropic rectangular plate hit by a concentrated laser pulse. 2. MATHEMATICAL MODEL 2.1. Plate and impact models An anisotropic rectangular plate is considered. Its side lengths are 2a and 2b and its thickness is h. Cartesian coordinates are introduced such that the plate is defined in the region G"M(x, y, z)3R3 D(x, y)3)1, z3[!h/2, h/2 ]N Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) ( 1997 John Wiley & Sons, Ltd. A CHEBYSHEV COLLOCATION MULTIDOMAIN METHOD 3691 where )1"M(x, y)3R2D x3[!a, a], y3[!b, b]N is the plate surface. The boundary ! of )1 is looked upon as a union of four sides and four corner %95 points, and the open region ) is defined such that )1")X! . %95 The dynamics is described by the Reissner—Mindlin plate theory using notations for the kinematics taken from a textbook by Hughes.8 Thus, at time t the displacements in x-, y- and z-directions are u "!zh (x, y, t) 1 1 u "!zh (x, y, t) (1) 2 2 u "w (x, y, t) 3 Note the signs chosen for the rotations h and h ; slope is the sum of fibre rotation and shear 1 2 strain. Indices 1, 2 and 3 are alternative representations for indices x, y and z, respectively. Note also that a comma means differentiation; f "Lf/Lx, f "Lf/Ly and f "Lf/Lt for any function f. ,1 ,2 ,t The plate moments m , m ("m ), and m and shear forces q and q are related to h , h and 11 12 21 22 1 2 1 2 w by constitutive equations: m "!D h !D h !D h !D h 11 11 1,1 16 1,2 16 2,1 12 2,2 m "!D h !D h !D h !D h 22 12 1,1 26 1,2 26 2,1 22 2,2 m "!2D h !D h !D h !2D h (2) 12 16 1,1 66 1,2 66 2,1 26 2,2 q "!A h #A w 1 55 1 55 ,1 q "!A h #A w 2 44 2 44 ,2 D , D , D , D , D and D are the flexural moduli, while A and A are the in-plane 11 12 16 22 26 66 44 55 moduli. As can be seen from the paper by Mindlin9 A and A "(12/h2) D for an isotropic 44 55 66 plate, if no shear correction factors are included when defining A and A . For the plates 44 55 investigated A "0, so that the corresponding terms in the expressions for q and q have been 45 1 2 omitted in equation (2). The impact is modelled by a transverse force per unit area, q. No other forces or couples are assumed to exist; the plate is free. Let q"S (x, y) ¹ (t ) (3) where x2#y2 exp s x2#y2!r2 S (x, y)" 0, CA BD , x2#y2(r2 and (x, y) 3 )1 x2#y2 *r2 and (x, y) 3 )1 and G ¹(t)" ( 1997 John Wiley & Sons, Ltd. q sin (nt/q ), 0 i 0, 0)t)q i q (t i Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) 3692 B. KJELLMERT The constants s and r determine the function S, which is smooth (3C=) and vanishes outside a circle with radius r. By choosing s and r one may simulate the spatial dependence of the actual impact pressure. Sometimes one also knows from experiments the impact time q and the total i momentum mv transferred to the plate by the impact. It is then easy to obtain q from mv by 0 integrating q analytically over the impact time and numerically over the plate area. 2.2. Initial-boundary value problem Initially, when the plate is at rest (t)0) h "0, 1 h "0, w"0, 2 (x, y) 3 )1 (4) When t'0 the dynamics in the plate is described by8 m !q "( oh3/12) h ab,b a a,tt , (x, y ) 3 ) q #q"( oh) w a,a , tt (5) and the dynamics on the free boundary by m n "0, ab b q n "0, a a (x, y)3! (6) %95 The plate surface is considered to be composed of a number of rectangular subdomains. On the boundary between two subdomains, e.g. subdomain u and v, essential and natural conditions have to be satisfied.8 They read (a"1, 2 and b"1, 2): hu!hv"0, wu!wv"0, a a (x, y) 3 ! */5 mu nu !mv nv "0, qu nu!qv nv"0, a a a a ab b ab b (7a) (7b) Here ! is the union of all boundaries between subdomains, and a repeated index a or b means */5 summation. n is the unit vector, outward to the considered subdomain. n is not unique at corner a a points. However, by assumption, n at a corner forms the angle 3n/4 with the sides meeting at the a corner. There is one exception to this prescription. If a subdomain corner (x, y) is situated on the side of the plate ((x, y) 3 ! W ! ) the vector n is outward and normal to this side. %95 */5 a Equations (2)—(7) yield an initial-boundary value problem adapted to a multidomain decomposition. 2.3. Subdomains and dimensionless variables First, dimensionless lengths and times are introduced by the replacements xPx/a, yPy/b, tPt/s (8a) where s"1 second. The plate surface )1 is hence mapped on a biunit square. This is covered by P · Q rectangular subdomains of equal size, Q subdomains in x-direction and P in y-direction, and P"Q"1 or Q"3, 5, . . . and P"3, 5, . . . . The subdomains are labelled in ascending order from right to left and from top to bottom. See Figure 1. Each rectangle is then mapped on a biunit Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) ( 1997 John Wiley & Sons, Ltd. A CHEBYSHEV COLLOCATION MULTIDOMAIN METHOD 3693 Figure 1. Enumeration of subdomains (Q"P"3) and definition of x- and y-directions computational square, and the first derivative operators change L 1 L L 1 L P , P Lx Q Lx Ly P Ly (8b) when local coordinate systems are introduced; the second derivative operators change in an analogous way. (The origo of each local xy-system is at the centre of the computational biunit square.) Secondly, dimensionless dependent variables are defined: ºX"h , 1 º½"h , 2 ºZ"w/(ab)1@2 (9) 3. NUMERICAL PROCEDURE A Chebyshev collocation multidomain method has been developed to solve the problem defined above by equations (2)—(9). The solver will be described in Sections 3.1—3.4. 3.1. Space discretization and patching Each subdomain is discretized using the Gauss—Lobatto—Chebyshev collocation points. Let there in each subdomain be M#1 points in x-direction and N#1 points in y-direction. M and N are even numbers. The subdomain collocation points are given in Figure 2. When patching the subdomains one must ascertain that the number of conditions equals the number of unknowns. Furthermore, it is appropriate to treat separately side points and corner points. Some side points are on the free edge of the plate (pu 3 ! ), others are shared by two ij %95 subdomains (pu is physically the same point as pv for the subdomains u and v, and some index kl ij pairs (i, j ) and (k, l), and pu 3 ! ). If the point is on the free boundary (pu 3 ! ) equation (6) ij %95 ij */5 applies and if the point is on the internal boundary (pu 3 ! ) equations (7a) and (7b) apply. There ij */5 are three kinds of corner points: (a) (b) (c) pu 3 ! , external corner point ij %95 pu and pv 3 ! W! , corner point at the side of the plate kl %95 */5 ij pu and pv 3 ! , internal corner point kl */5 ij ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) 3694 B. KJELLMERT Figure 2. Pseudocodes for creation and enumeration of the subdomain collocation points p ("(x , y ), local co-ordinate ij i j system) Remember that in cases (a) and (c) the unit normal vector is directed diagonally and outwards, while it is normal to the edge of the plate in case (b). Let us now define the patching at corner points. (a) At an external corner point equation (6) applies. (b) Consider first a plate with P"Q"3 (see Figure 1). At two physically coinciding points of type (b) equation (7a) is satisfied. At one of them equation (6) is also satisfied. By assumption, subdomains 1, 3, 7 and 9 have two points of type (b) where equation (6) applies, and none of subdomains 2, 4, 6 and 8 have points of this kind. In other words, along the edges of the plate every second subdomain has two points of type (b) where equation (6) applies, and among them are the subdomains situated at the corners of the plate. Obviously this prescription also works generally provided P and Q are odd. (When P"Q"1, no points of type (b) exist.) (c) Four subdomain points coincide physically at an internal corner point. The subdomains meeting there are looked upon as two diagonally coupled pairs, coupled at physically the same point. Each pair has together two points of type (c) where equations (7a) and (7b) apply. This prescription is consistent with the directions of the unit normal vector, and implies for each pair a coupling of all points on the four subdomain edges meeting at the corner. The present patching method has been chosen since it fulfils symmetry requirements and yields a rather simple code. Further information on this subject can be found in the monograph by Canuto et al.5 The dependent variables are classified like the collocation points as interior (B), side (S), and corner (C) ones. Let º be the vector of all variables at all collocation points. Let º"(ºB, ºS, ºC)T, where ºB"(ºB1, . . . , ºBPQ)T, ºS"(ºS1, . . . , ºSPQ)T, ºC"(ºC1, . . . , ºCPQ)T. The values of º in the interior of subdomain u ("1, . . . , PQ) are ºBu"(ºXu, º½ u, ºZu)T where the vectors ºXu, º½ u, ºZu are ordered in the same way as the interior collocation points. Similarly, ºSu"(ºXu, º½ u, ºZu)T and ºCu"(ºXu, º½ u, ºZu)T for all values at side points and corner points, where the vectors ºXu, º½ u, ºZu each are ordered in the same way as the corresponding collocation points. Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) ( 1997 John Wiley & Sons, Ltd. A CHEBYSHEV COLLOCATION MULTIDOMAIN METHOD 3695 3.2. Matrix formulation The initial-boundary value problem is time discretized using the trapezoidal rule. Let the time step be *¹. At time t"N¹ · *¹ the time discretized versions of equations (2)—(7) modified by (8)—(9) form a boundary value problem where the input depends on the solution at the previous time level. N¹ is an integer less than approx. 2000. This problem is space discretized using the grid and patching described in Section 3.1. For each subdomain the differential operators L/Lx, L/Ly and L2/Lx2, L2/LxLy, L2/Ly2 are transferred into matrices by matrix representations for the collocation derivatives.4 The matrix equation for the vector º at time N¹ · *¹ takes the form A ¸BB ¸BS ¸BC ¸SB ¸SS ¸SC BA B A B ºB RB ºS " RS ¸CB ¸CS ¸CC ºC (10) RC where the vector RB on the right side depends on the external pressure and the previous vector º. RS and RC are both null-vectors, the boundary conditions being time-independent and independent of pressure. ¸BB, ¸BS and ¸BC correspond to equation (5) for the interior of the subdomains. Since interior points of two subdomains are not connected, ¸BB, ¸BS and ¸BC are diagonal block matrices. Let the blocks of ¸BB, ¸BS, and ¸BC be denoted MBB, MBS and MBC, respectively. MBB, MBS and MBC are the same for all subdomains. MBB is a 3 · (M!1)· (N!1)]3 · (M!1) · (N!1) matrix. Moreover, the condition number of MBB is rather small and its inverse, MBB~1, can be computed. ¸SB, ¸SS and ¸SC describe the patching conditions along the internal sides (see equation (7)) and the free boundary conditions at the external sides (see equation (6)). Similarly, ¸CB, ¸CS and ¸CC describe the corner conditions. One observes that ¸CB has only zero elements, since the differential equations corresponding to the boundary conditions contain no L2/LxLy-operator. Observe that the patching conditions given above are direct translations of the strong formulations in equations (6) and (7). No flux conditions in the sense of Macaraeg and Street10 are used. ºB is eliminated from equation (10) using block Gaussian elimination. One finds A GSS GSC ¸CS ¸CC BA B A B ºS QS " ºC RC (11) The Gauss transforms GSS ("¸SS!¸SB · ¸BB~1 · ¸BS), GSC ("¸SC!¸SB · ¸BB~1 · ¸BC) and QS ("RS!¸SB · ¸BB~1 · RC) are easily computed, the matrix ¸BB~1 being block diagonal with blocks MBB~1. GSS is a rather well-conditioned matrix of size 6PQ (M#N!2) ]6PQ (M#N!2). The structure of GSS is illustrated in Figure 3. Another block elimination reduces equation (11) to GCC · ºC"QC (12) Here GCC ("¸CC!¸CS · GSS~1 · GSC) is a full 12PQ by 12PQ matrix. It is rather ill-conditioned. The right side QC ("RC!¸CS · GSS~1 · QS) contains information coming from RB defined in equation (10). ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) 3696 B. KJELLMERT Figure 3. MATLAB spy plot of GSS. Non-zero elements are black. M"N"6, Q"P"3 3.3. Matrix properties The inverses of the matrices MBB, GSS and GCC are parts of the Gauss transforms which form the basis of the present method. It is hence of interest to study these matrices. Some results for their condition numbers, calculated by a standard ¸º-solver,11 are shown in Table I. (The corresponding plate data are given in Section 4.) The matrix MBB is a sum of a truncated stiffness matrix corresponding to the left-hand side of equation (5) and a mass matrix corresponding to the right-hand side of the same equation. The stiffness matrix has a rather big condition number, which is proportional to N4. The mass matrix is well conditioned. It is diagonal and proportional to *¹ ~2. Therefore, the sum MBB becomes better conditioned as the time step *¹ decreases (see Table I). However, letting the norm of MBB be equal to 1, the stiffness elements of MBB are of the order of the machine’s floating-point precision if *¹ is too close to zero. The choice of *¹ is further discussed in Section 3.4. The matrix GSS ("¸SS!¸SB · ¸BB~1 · ¸BS) describes the side conditions for all subdomains. In addition, GSS takes into account the bulk equations and depends on *¹ through the matrix MBB~1. Also GSS becomes better conditioned as *¹ decreases; GSS is singular when *¹ ~1 equals zero. The matrix GSS is sparse and banded (see Figure 3). Although GSS~1 is also needed at each time-level, one does not compute it, but instead uses the ¸º-factors of GSS to obtain QC ("!¸CS · GSS~1 · QS) and GCC ("¸CS · GSS~1 · ¸BS). The sparse solver of MA¹¸AB12 and the banded solver of Numerical Recipes13 turn out to be less effective for computing the ¸º-factors than a modified version of a standard ¸º-solver.11 This solver has been altered so that multiplying by some zeros is avoided. (The computations are run in double precision on a SºN-SPARC4.) Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) ( 1997 John Wiley & Sons, Ltd. 3697 A CHEBYSHEV COLLOCATION MULTIDOMAIN METHOD Table I. Condition numbers of matrices MBB, GSS and GCC · P"Q"3. M"N N 6 12 18 *¹ 0·50!5 0·50!6 0·50!7 0·50!5 0·50!6 0·50!7 0·50!5 0·50!6 0·50!7 Cond (MBB) 0·98#4 0·29#2 0·12#1 0·29#5 0·10#3 0·20#1 0·88#5 0·27#3 0·36#1 Cond (GSS) 0·12#5 0·55#3 0·59#2 0·10#6 0·19#5 0·23#3 0·18#6 0·52#5 0·50#3 Cond (GCC) 0·23#10 0·21#10 0·25#10 0·84#10 0·82#10 0·85#10 0·18#11 0·20#11 0·20#11 Table II. Sizes of MBB, GSS and GCC Size (MBB) Size (GSS) Size (GCC) Q"P"7 M"N"8 Q"P"5 M"N"10 Q"P"3 M"N"18 Q"P"1 M"N"54 147]147 4116]4116 588]588 243]243 2700]2700 300]300 867]867 1836]1836 108]108 8427]8427 636]636 12]12 The matrix GCC describes the corner conditions for all subdomains, taking into account the side and bulk equations. For realistic values of *¹, GCC is worse conditioned than the matrices MBB and GSS. However, as a rule GCC is much smaller than MBB and GSS (see Table II). Cond (GCC) is of the order 1010, but the number of significant figures of the solution ºC ("GCC~1 · QC) can be recovered by one step of iterative improvement. According to the definitions of GCC and GSS, GCC"¸CC!¸CS · (¸SS!¸SB · ¸BB~1 · ¸BS)~1 · GSC, where ¸BB~1 is formed by MBB~1 blocks which depend on *¹. This expression suggests that the matrix GCC is more related to the boundary than to the bulk. A study of Table I confirms these ideas. Indeed, Table I shows that the condition number of GCC, Cond (GCC), depends weakly on *¹. It is also clear from Table I that Cond (GCC) is proportional to N2. This behaviour is representative of first-derivative operators. In conclusion, the time-independent, first-order equation (6) for the boundary is more important for Cond (GCC) than the time-dependent, second-order equation (5) for the bulk. The Reissner—Mindlin theory is the origin of the large condition numbers. The theory gives rise to a stiff system of equations (5) and (6), provided the plate is thin (h2;ab). The coefficients in equation (5) are then of different orders of magnitude, while the coefficients in equation (6) are of the same order, and large condition numbers are generated. This drawback of the Reissner— Mindlin theory has recently been discussed6 for a plate modelled by spectral elements. It was shown, among other things, that the stiffness matrix becomes worse conditioned as the plate thickness diminishes. 3.4. Choices of M, N, P, Q and *¹ The governing parameters are chosen starting out from characteristic lengths and times which originate from the descriptions of the plate and the impact. ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) 3698 B. KJELLMERT The impact pressure is concentrated over a little spot at the centre of the plate. Let the characteristic impact length be j . For the plate investigated, h(j (a, b, and waves with * * wavelengths of the order j are to be resolved. The distance between the collocation points must * then be less than j /3. This is secured by large enough values of Q, P and/or M, N. Large values of * these parameters give, however, large matrices. Table II shows the sizes of MBB, GSS, and GCC for four choices of Q, P and M, N, choices which yield almost the same space resolution (+0·1 a+0·1 b). As a rule of thumb, the present solver works well when the sizes of MBB and GSS are of the same order of magnitude, and GCC is rather small. The best alternative in Table II is Q"P"3 and M"N"18. Two characteristic times influence the choice of the time step *¹, viz., an impact time and a cut-off time. The impact time q is determined by the duration of the laser pulse and of the i heating and recoiling of particles in the spot hit by the pulse. To model the impact one requires that *¹;q . The cut-off time comes from the differential equations describing the plate * dynamics, (2) and (5). These equations are three coupled wave equations. Therefore, in an arbitrary direction, three plane waves may propagate in the Reissner—Mindlin plate; cf. e.g. Reference 14. The corresponding dispersion relation, u"u (k), has three branches, i.e. given a wave number k three angular frequencies u are possible. A bending mode exists for all angular frequencies, and two thickness shear modes exist for angular frequencies above a cut-off value, u . (Here it is assumed that A44"A55, so that only one cut-off frequency appears; #0 u2 "12·A44/(oh3) and u2 "12·A55/(oh3).) Thus 2n/u is a characteristic time, and to #02 #0 #01 account for the shear modes one requires: *¹;2n/u . #0 The trapezoidal rule is used for time integration. This rule implies no algorithmic damping because of the time-discretization, but relative period errors (rpe) appear. They increase with *¹. According to the test oscillation problem described in the textbook by Hughes,8 *¹ (u/2n))0·125 for rpe)0·05. A reasonable choice of *¹ is thus *¹+0·125·2n/u #0 (+0·44]10~6 for the anisotropic plate defined in Section 4.) As can be seen from Section 3.3 this value of *¹ is so large that one takes into account the stiffness elements of matrix MBB, and it is so small that the condition numbers of MBB, GSS and GCC are acceptable. 4. RESULTS Computations have been run to simulate the behaviour of a glass-fibre reinforced epoxy plate impacted by a laser pulse. The pulse hits the plate at the centre of its upper surface, and the impact time and length are known approximately: q +10 ks and j "5—10 mm. This paper examines * * a specific plate with the following data (SI-units): a"0·15, b"0·10, h"0·0041, D "173·0, 11 D "77·0, D "27·0, D "!18·9, D "3·0, D "51·0, A "A "12D /h2 , and 22 12 16 26 66 44 55 66 A "0. The bending moduli have been determined through the Kirchhoff plate theory, while the 45 in-plane moduli have been arrived at through the assumption that only the isotropic bulk material contributes to them. No shear correction factors are introduced. The contour plots in Figure 4 show the out-of-plane displacement field (ºZ) and the corresponding velocity field (»Z). According to the trapezoidal rule, at time t"N¹ · *¹, &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&" Figure 4. Contour plots of displacement field (ºZ, left) and displacement velocity field (»Z, right) in out-of-plane direction at times 0, 40, 80, 120 ks after impact. Time axis downward. q "10]10~6, s"100, r"b. *¹"0·5]10~6. * M"N"18, Q"P"3 Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) ( 1997 John Wiley & Sons, Ltd. A CHEBYSHEV COLLOCATION MULTIDOMAIN METHOD ( 1997 John Wiley & Sons, Ltd. 3699 Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) 3700 B. KJELLMERT (*¹/2)»Z (N¹ )"(ºZ (N¹ )!ºZ (N¹!1))!(ºZ (N¹!1)!ºZ(N¹!2))#· · · ! (!1)NT(ºZ(1)!ºZ(0)), so that the right column of the plots also gives information on the field ºZ. At t"0 ks, just after impact, »Z is similar to ºZ and both fields are near pictures of the impact pressure. (Note the new definition of time t"0; the impact starts at t"!q and ends at * t"0.) Then the displacement waves propagate outwards, but at t"40 ks they have not yet reached the boundary of the plate. It can be seen, however, that small amplitude waves move ahead of the visible displacement, the velocity disturbance being spread over a larger portion of the plate than the displacement disturbance. At t"80 ks the waves have reached most of the boundary, and at t"120 ks the waves have reached the whole boundary and some fast reflected waves have returned to the centre of the plate. When producing Figure 4 the preparatory computations require 2260 s on a SUN-SPARC4, using double precision fortran code. The following computations require 24·8 s per time-level. The Chebyshev coefficients of ºZ and »Z are determined and stored at desired time-levels. Then a simple MATLAB program yields the contour plots starting from the Chebyshev coefficients. The accuracy of the algorithm is estimated by some indirect methods. First, the mean out-of-plane velocity should be constant after impact. In fact, the plate is at rest initially. Therefore, the total linear impulse on the plate, applied to it by the pulse, equals the linear momentum of the plate just after impact. This linear momentum is conserved, the plate being free. Hence, also the average out-of-plane velocity, »ZM, is conserved. To find »ZM numerically Gauss-Lobatto integrations of »Z are performed in x- and y-directions over each subdomain. »ZM is then calculated as the mean of all subdomain means of »Z. A time-history of »ZM is given in Figure 5. The linear impulse on the plate has been normalized so that »ZM should equal 1 after impact. Figure 5. Time history of the mean of the out-of-plane velocity, »ZM. For other data, see Figure 4 Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) ( 1997 John Wiley & Sons, Ltd. A CHEBYSHEV COLLOCATION MULTIDOMAIN METHOD 3701 Figure 6. (a) Contour plot of the displacement field at t"40 ks according to the Kirchhoff plate theory and an analytical transform solution. (b) Photo of the plate considered at t"40 ks after impact. The photo was produced by holographic interferometry Choosing *¹"0.5]10~7 yields a »ZM-curve almost identical to the one in Figure 5, at least for the first 80 ks after impact. In conclusion, the time integration algorithm seems to be stable, and the time step *¹"0·5]10~6 is satisfactory for the considered case. Secondly, the deformation of the plate has also been calculated from the Kirchhoff plate theory using an analytic transform method.3 Here it is assumed that the laser pulse is a delta-shaped pressure in space and time and that the boundary of the plate can be ignored. However, the initial ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) 3702 B. KJELLMERT development of the deformation field should be modelled satisfactorily by this procedure. The contours obtained by the Kirchoff model at t"40 ks are shown in Figure 6(a). This picture agrees rather well with the corresponding plot of the ºZ-field in Figure 4. Thirdly, direct observation of the deformation field of the plate is also possible, using holographic interferometry.1 Of course, comparing experimentally produced photos with the computed contour plots of the displacement field gives an ultimate test of both the basic Reissner—Mindlin plate theory and its numerical implementation. In Figure 6(b) is shown the experimentally obtained photo at t"40 ks, a photo which is consistent with the contour plot of ºZ at t"40 ks in Figure 4. The agreement is perhaps better than it should be; the laser spot is somewhat less than that defined by the parameters s"100 and r"b of the simulation, and the impact time is longer than 10 ks according to some experimental data. In addition, one side of the real plate is partly clamped, while in the model all sides are free. This discrepancy means little soon after impact (t"40 ks), and surprisingly little when the waves have reached the sides (t"80 ks). Several other plates have also been studied using the present solver, a transform solution of the transient development, and the hologram interferometry technique. In particular orthotropic plates seem to be modelled very accurately. The in-plane moduli A and A cannot, however, be 44 55 determined by this way of working. What is worse, computations show that the values of A and 44 A do influence the simulation of the behaviour of the transient waves. 55 ACKNOWLEDGEMENTS Karl-Evert Fällström and Anders Wa> hlin of Lulea> University, Div. of Experimental Mechanics, have provided Figures 6(a) and 6(b), respectively. The author has discussed this work with Karl-Evert Fällström. The reviewers’ comments have improved this paper. REFERENCES 1. K.-E. Fällström, H. Gustavsson, N.-E. Molin and A. Wa> hlin, ‘Transient bending waves in plates studied by hologram interferometry’, Exp. Mech., 29, 378—387 (1989). 2. M. El-Raheb, ‘Flexural waves in a disk soon after impact’, J. Acoust. Soc. Am., 96, 221—234 (1994). 3. K.-E. Fällström, O. Lindblom, R. Näslund, L.-E. Persson, ’A study of vibrations in infinite plates’, (1995), submitted. 4. H. L. Ku and D. Hatziavramidis, ‘Chebyshev expansion methods for the solution of the extended Graetz problem’, J. Comput. Phys., 56, 495—512 (1984). 5. C. Canuto, M. Y. Hussaini, A. Quarteroni and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer, New York, 1988. 6. U. Zrahia and P. Bar-Yoseph, ‘Plate spectral elements based upon Reissner—Mindlin theory’, Int. J. Numer. Meth. Engng., 38, 1341—1360 (1995). 7. C. A. J. Fletcher, Computational Galerkin Methods, Springer, New York, 1984. 8. T. J. R. Hughes, ¹he Finite Element Method, Prentice-Hall, Englewood Cliffs, N.J., 1987. 9. R. D. Mindlin, ‘Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates’, J. Appl. Mech., 18, 336—343 (1951). 10. M. G. Macaraeg and C. L. Street, ‘Improvements in spectral collocation through a multiple domain technique’, Appl. Numer. Math., 2, 95—108 (1986). 11. G. E. Forsythe, M. A. Malcom and C. B. Moler, Computer Methods for Mathematical Computations, Prentice-Hall, Englewood Cliffs, N.J., 1977. 12. MATLAB, »ersion 4.2a, 1994. MATLAB is a trademark of The Math Works, Inc. 13. W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes, 2nd edn., Cambridge University Press, New York, 1992. 14. T. S. Chow, ‘On the propagation of flexural waves in an orthotropic laminated plate and its response to an impulsive load’, J. Comput. Mater., 5, 306-319 (1971). . Int. J. Numer. Meth. Engng., 40, 3689—3702 (1997) ( 1997 John Wiley & Sons, Ltd.