close

Вход

Забыли?

вход по аккаунту

?

512

код для вставкиСкачать
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.
Документ
Категория
Без категории
Просмотров
2
Размер файла
329 Кб
Теги
512
1/--страниц
Пожаловаться на содержимое документа