Initial and Boundary Conditions in ... - Semantic Scholar
1648
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
Initial and Boundary Conditions in Multidimensional Wave Digital Filter Algorithms for Plate Vibration Chien-Hsun Tseng and Stuart Lawson, Senior Member, IEEE
Abstract—Recently, multidimensional wave digital filter (MDWDF) structures have been proposed for the modeling of plate vibration problems. In this paper, we discuss how initial and boundary conditions may be properly embedded into such an algorithm in terms of the state quantities that are an integral part of the algorithm. Due to the essential feature of fully-local interconnectivity in the MDWDF model, different types of boundary conditions can be easily satisfied in a very simple and efficient manner. Instead of remodifying the whole algorithm, usually required by finite elements based methods, boundary conditions in terms of state outputs are simply attached to the model. This feature is especially useful when dealing with the mixed-edges boundary conditions frequently encountered in practice. Graphical results obtained from implementing the MDWDF algorithm are given to further demonstrate the capacities of the method in efficiently handling a fourth-order Mindlin plate vibration system with various types of boundary conditions. Index Terms—Courant–Friedrichs–Levy bound, golden section search method, Kirchhoff circuit, losslessness, Mindlin system, multidimensional wave digital filters (MDWDF), partial differential equations (PDEs), passivity, stability, trapezoidal rule.
I. INTRODUCTION
T
HE modeling of plate vibration in mechanics is of great importance to applications in engineering and science, especially in the areas of aerospace, marine and construction. The vibration aspects in micorelectromechanical systems applications such as electronic packaging, micro-robots, etc. are increasingly important because of their enhanced sensitivities to vibration. The most popular kind of models for mechanical vibration systems can be represented by sets of linear and/or nonlinear partial differential equations (PDEs) with properly imposed initial and boundary conditions. In most cases of complex vibration-influenced design, the equations can not be solved analytically and so numerical approaches are sought. As far as practical interest is concerned, it is, therefore, desirable to numerically simulate the behavior of such PDE systems as cheaply and rapidly as possible, especially when dealing with complicated boundary conditions. Numerical techniques used to solve the PDEs such as finite elements and finite differences have had some success. The finite-element methods permit easy inclusion of local grid refinement and handling of complex geometries. These methods, Manuscript received August 12, 2002; revised May 18, 2003. This work was supported by EPSRC(UK) under Grant GR/N22298. This paper was recommended by Associate Editor M. Nakhla. The authors are with the School of Engineering, University of Warwick, Coventry CV4 7AL, U.K. (e-mail: [email protected]; [email protected]; [email protected]) Digital Object Identifier 10.1109/TCSI.2004.832796
however, are computationally expensive and harder to correctly set up the simulation plane than the frequently used finite difference methods. The latter approaches, on the other hand, have difficulties in handling irregular boundaries, and need local grid refinement to increase the accuracy. A remarkable alternative approach to the modeling and simulation of a system of PDEs builds on the wave digital filter (WDF) paradigm and analogies with electrical networks, an area rich in mathematical structure [1], [2]. By using such a paradigm and properties of electrical networks, concepts well known to circuit theorists such as passivity, losslessness, stability, finite propagation velocity, etc., can be translated to the actually considered physical problems so as to preserve important relationships between variables. The technique involves first finding a multidimensional lumped electrical network which represents the behavior of the system. From this network, a discrete-time equivalent is developed that is a multidimensional WDF (MDWDF) [3]–[5]. Proceeding in this way, one can substantially achieve for the resulting algorithm: • massive parallelism because each point in the -D grid can be updated simultaneously if sufficient computing resources are available [3], [6]; • high accuracy because of making use of the WDF structure known to have low roundoff noise characteristics and suppression of parasitic oscillations [2], [7], [8]; • multidimensional passivity (reflecting the passivity of the system to be modeled) that guarantees robustness and basically all numerical stability properties required of an accurate numerical integration method [9], [10]; • full local interconnectivity due to the generation of a second-order difference equation for which its behavior at any point in space is directly influenced only by the points in its nearest neighborhood [4], [11]. This local interconnectivity is particularly helpful to the establishment of the boundary conditions in a very simple manner. Other second-order methods do not offer this property [12]. The technique of MDWDF originally introduced for the implementation of digital filters [1], [2], [13] has been successfully applied to the numerical integration of PDEs [3], [5], [14]–[16], etc. These PDEs systems usually involve low-order theories for the governing equations. As a result, the systems can be directly handled by a standard PDE solver such as the MATLAB PDE Toolbox [17]. In this paper we discuss in particular the crucial issue of the setting of initial and boundary conditions with regard to a recently proposed MDWDF algorithm which simulates plate vibration. We begin by considering the thick plate theory due to Mindlin [18]–[20], the first approximate theory for dealing especially
1057-7122/04$20.00 © 2004 IEEE
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
1649
with flexural waves in relatively thick plates, although thin plates can also be dealt with. We note that a plate is typically considered to be thin if the ratio of its thickness to representative lateral dimension (e.g., circular plate diameter, square plate side length) is 1/20 or less. Modeling of the Mindlin system had been achieved using many different techniques such as the finite element, the boundary element, the finite difference, the differential quadrature, the collocation, the Galerkin and the Ritz methods [21], [22]. The idea of using a digital network to model Mindlin plate vibration was first reported in [16] where the networks used were based on: 1) digital wave guides (DWN) and 2) MDWDF using the space transformation method [3]. Good results were obtained using the DWN model only, the MDWDF model did not yield results due to problems in the formulation of the initial and boundary conditions [16]. This paper begins with the modeling of the Mindlin system in detail. In Sections II and II, we present a variant of the model discussed by Bilbao [16]. Following a direct space–time analysis methodology [4], we start from the system of PDEs and derive the equivalent analogue network description. This results in a MDWDF algorithm. A useful feature of MDWDF modeling is the ability at this stage to satisfy the system’s stability requirement in order to ensure the passivity of the circuit obtained. The next step is to properly incorporate the given initial and boundary conditions into the open-field MDWDF algorithm obtained. It turns out that due to the advantage of fully local interconnectivity of the MDWDF approach, a straightforward and efficient way to accommodate various boundary conditions without remodifying the whole MDWDF algorithm, is provided. Other finite element-based methods do not have this advantage. Finally, two cases of plate vibration are considered to demonstrate capacities of the method in handling the high-order governing equations. II. PROBLEM STATEMENT The stress equations of motion of a two-dimensional (2-D) Mindlin plate can be written as a system of PDEs [18], [20]
(2.1)
where is the plate thickness, are the bending rotations of a transverse normal to the and axes, respectively, are the transverse shear forces per unit length of the are the bending moments per unit plate, and length of plate. The plate material is assumed to have the density , Young’s modulus/modulus of elasticity in tension and com. pression , as well as Poisson’s ratio , such that Assuming a plane stress distribution in accordance with Hooke’s law, the plate stress-displacement relations are then given by
(2.2) is a shear correction factor being applied to compenwhere sate for an error in assuming a constant shear stress through is the modulus of elasthe plate thickness, is the flexural rigidity ticity in shear, and of the plate. When the stress-displacement equations in (2.2) are incorporated in the system defined by (2.1), we obtain the governing equations of motion described as a system of eight , , and PDEs as follows by defining : (2.3)
(2.4)
In view of (2.3) and (2.4), it was noted in [16] that the complete Mindlin system can be decomposed into two separate subsystems where the corresponding system variables are written and . Furthermore, the as differentiable equations of motion can be reduced to a single equation for the displacement by firstly substituting the plate stress-displacement equations of (2.2) into (2.1). Followed up by differentiating these equations with respect to and then summing, one finally obtains a fourth-order Mindlin plate vibration system in time and space
(2.5) In order to apply the circuit theory for representing the behavior of the system in (2.3)–(2.4), the same physical dimension must be kept for all plate variables. As a result, scaled dependent combined with positive constants variables , 2, 3 are introduced as (2.6) where the transverse velocity , and velocities of transverse and are treated as the graphical voltage bending rotation denoted by corresponding to the graphical current denoted by . It was noted in [16] that the plate system in (2.3)–(2.4) can be replaced by the following form written as symmetric hyperbolic PDEs: (2.7)
1650
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
(2.8)
where the operators , denote the partial derivatives with respect to temporal and spatial coordinates in and axes. III. GRAPHICAL NETWORK DESCRIPTION AND MULTIDIMENSIONAL WAVE DIGITAL FILTERING ALGORITHM In this section, the derivation of a three-dimensional (3-D) Kirchhoff circuit is presented. Moreover, the linear discretization using a generalized trapezoidal rule is applied to each component of the Kirchhoff circuit for achieving a structure that can be numerically implemented in terms of difference equations.
Fig. 1. Multidimensional Kirchhoff circuit representation for Mindlin systems (2.3)–(2.4).
with respect In addition, the partial derivative operators to the two-port Jaumann structure [25] are given by defining arbitrary positive constants
A. Network Description The key success in the derivation of a MDWDF algorithm for the numerical solution of Mindlin system is the derivation of a MD Kirchhoff circuit from the underlying PDEs. By representing these equations (2.7)–(2.8) in graphical form, one can express mathematical manipulations as network operations known as the Kirchhoff circuit theory in the following:
Remark 3.1: In view of Fig. 1, we note that by applying the techniques of circuit theory, inductors of a MD network description resulted from the plate stress-displacement PDEs in (2.3)–(2.4) can not be realized by real-world components. Instead, the network description is a graphical representation of time and space differentiation operators for the system of (2.3)–(2.4). Furthermore, the graphical circuit representation in Fig. 1 has resulted in a convenient way of understanding the passivity concept and grid generation for MD physical systems [4], [9], [10]. B. MDWDF Algorithm
where , , and . Following the Kirchhoff’s laws, the above equations can be immediately interpreted in the form of a 3D Kirchhoff circuit represented in Fig. 1, see also [16, Fig. 5.16]. In view of the MDKC of Fig. 1, the conservation of energy holds as a result of Kirchhoff’s laws. Furthermore, the graphical circuit contains three main loops on the left, five main loops on the right, and two gyrators coupling in the middle with the gyration coefficients , and . The self-inductors with inductances are given by
To arrive at an approximation of each mesh voltage across the corresponding inductance in the 3D Kirchhoff circuit of Fig. 1, a generalized trapezoidal rule [3], [4] is applied to the integration of each mesh voltage involved in two dimensions where , take only discrete values with the temporal step size and the spatial and in and coordinates, respectively. As a step sizes, result, we obtain the port resistances and for the two-port Jaumann filters and self-inductors, respectively, given by
(3.2)
When the condition for the passivity of inductances is held in the MDKC of Fig. 1, i.e., , clearly the simplest choice for these is given by
(3.1)
(3.3)
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
1651
Fig. 2. MDWDF algorithm for numerical simulation of Mindlin system (2.3)–(2.4).
Substituting in
chosen by (3.3) into (3.1)–(3.2), it simply results , , , 3, 6, 7, and
now transform the MDKC of Fig. 1 into the MDWDF algorithm depicted in Fig. 2 (see also [16, Fig. 5.16]) where the shift operators are obtained by
(3.4)
where the plus sign “ ” occurs when . Following from (3.2)–(3.3), we are able to determine the resistances of port adaptors in Fig. 2 as shown in (3.5) at the bottom of the page. To guarantee that the gyrator couplings between the loops with , and in the MDKC of Fig. 1 are memorycurrents less, reflection-free ports are necessary [16] which consequently determined by leads to the port resistances
Defining and following from 1, 2, 3 in (3.1), we obtain a lower bound of a ratio between the density of sampling in space and that of sampling where is the in time, i.e., phase velocity of a wave [18] traveling freely across the whole plate according to the fourth-order system (2.5). This is the well-known Courant–Friedrichs–Levy (CFL) bound [26], a stability criterion for many approaches to the numerical solution of the wave equation. For a simplified reference circuit where , the are determined by (3.3) and hence resulting in the minimal ratio can be obtained by the golden section search method. Thus, it imposes the least restriction on the density of the sampling in time for a given density of the sampling in space so that the numerical stability of the Mindlin system can be guaranteed when simulated by a passive MDWDF circuit. and are minimized to Furthermore, the inductances achieve the losslessness of the system. Assuming that the voltage wave representation is adopted and making use of properties of the two-port lattice structure which is equivalent to the two-port Jaumann structure [25], we may
As a result, a delay-free directed loop is not created leading from and the input to the output terminal. Let be the currents corresponding to the wave quantities and , respectively, in the MDWDF algorithm of Fig. 2. By synthesizing various ports in the MDKC of Fig. 1, using the series form and given below in terms of currents
(3.6)
(3.5)
1652
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
the wave quantities of these series adaptors in the MDWDF of Fig. 2 can be described as
Referring back to the MDWDF algorithm of Fig. 2, the relations of wave and state quantities can be seen as
for
for
(3.7)
Of particular interest are the reflection-free ports in the gyrator couplings. According to (3.7), the then become independent . Let the vector of voltage waves [5] be of for and where the vector given by , and the vector voltage and current are defined waves by , , , and
(3.8) Thus, the coupled two-port series adaptor described by vector waves is obtained by (3.9) (3.10) where the vector waves and , and are and the coupling matrix of the port resistances given in (3.5). It is to be noted that vector wave variables were first used by Nitsche in his development of a MDWDF model for Timoshenko’s equations for beam bending [5]. They have also been used for the Mindlin plate model in [16]. Finally, for , 2 the two-port series adaptors (in the middle of Fig. 2) due to the gyrator couplings are given by (3.11) (3.12) where
(3.20) Furthermore, the state input-output relations are given by
(3.21) where the plus sign “ ” occurs when and and outputs.
, are vectors of state
IV. INITIAL AND BOUNDARY CONDITIONS In this section, we discuss how to accommodate the properly imposed initial and boundary conditions into an open-field MDWDF algorithm depicted in Fig. 2 in terms of state equations given by (3.21). A. Initial Conditions Let initial conditions of the Mindlin plate variables given in (2.6) be imposed as follows:
(4.1)
where , and . We first consider the first derivative initial conditions . Replaced the plate variables in (3.13)–(3.19) by their initial values in (4.1) and following from the results obacross each selftained in [15], we have the forward waves in the MDKC of Fig. 1: inductor with inductance
are given by
We further consider the vector of voltage a coupled inductor with inductance series connections of and , we have
(4.2) given in (3.8) across . Due to , by defining
The desired plate variables following from (3.5)–(3.12) can now and (indices be expressed in the form of forward waves are omitted) as (3.13) (3.14) (3.15) (3.16) (3.17) (3.18) (3.19)
Replaced by (3.10) and substituting by (3.18) with their , 7, it holds that initial values (4.3) Taking sums and differences of the plate variables replaced by their initial vales in (3.14)–(3.17), and then combining together given in (4.2), with (3.13), (3.19), (4.3), and
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
Fig. 3.
1653
Decomposed rectangular Mindlin plates with boundaries calculated by state output on edges and vertices.
we obtain two sets of state output equations, respectively, corresponding to initial conditions of the plate variables in two sepand (inarate subsystems dices and are omitted)
(4.4)
(4.5)
Let , , , and , where are given function of . Solving (4.4)–(4.5), we thus obtain the remaining state outputs shown in the equation at the bottom of the page. B. Boundary Conditions There are four common types of boundary conditions applied to the plate edge of Mindlin plates [23], [22], [20]: Free edge Hard-type simply-supported edge Soft-type simply-supported edge Clamped edge.
Considering a rectangular Mindlin plate whose grid points are calculated by state outputs (3.21) which involve both the spatial shifts and time delay, it seems clear that the plate can be decomposed into three separate subplates depicted in Fig. 3. Starting with the straight-line plate boundaries parallel with either the or axis illustrated in Fig. 3, clearly one of four neighborhood points is outside the boundary of each subplate. This has resulted in four combinations of state in, , , and , puts respectively, corresponding to the southern boundary (SB), western boundary (WB), northern boundary (NB), and eastern is one spatial step boundary (EB) where each state input size beyond the boundary of the rectangular plate. Instead of being obliged to determine sets of state inputs which are not computable, our aim is to find expressions for the corresponding set of state outputs that accommodate the given boundary conditions. Consequently, a state is enough for each equation determining each state output straight-line edge of the plate. This, however, will require us to impose at most one additional boundary condition. Similar to the straight-line edge grid points, two of four grid for each subplate. As a result, points are outside the vertex at most two additional boundary conditions have to be imposed for providing vertices of these subplates. Thus, each set of veris determined by six state output equations. tices By assembling these three subplates into one complete Mindlin plate, we, therefore, can obtain vertices of the rectangular plate , , , , which are included in the directions of SB, WB, NB, and EB, respectively. In the following, we will discuss how to find these additional conditions that associate with various boundary conditions given above. Substituting (3.20) into the plate variables described by the state equations in (3.13)–(3.19) and replaced the given boundary values, we then obtain the following state
1654
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
output equations directly corresponding to various boundary conditions
that, gyrator coupling between the loops with currents and and , taking respectively, determine the plate variables sums and differences of (4.8) and (4.10) yields
(4.6) (4.15) (4.7)
(4.8)
(4.9)
Combining the boundary conditions (4.11) and (4.12) with two additional conditions which result in (4.7) and (4.15), one is able to assemble these three subplates into one complete Mindlin plate so that the state outputs appearing at for the plate are obtained as foleach set of vertices lows:
(4.10) (4.11) (4.12) We note that (4.11) is a result of mutual inductors occurring in a coupled two-port series adaptor that is described by the plate and in the MDWDF algorithm of Fig. 2. variables 1) Free Edges (F-F-F-F): Following the conditions from (4.8), and (4.11)–(4.12), the state outputs outside the straight-line edge exclusive the vertices are obtained for each subplate by providing the additional conditions as follows:
(4.16)
(4.13) (4.17) and as shown in (4.14) at the bottom of the page. It is worth noting that the additional condition is given to provide the state outputs and for EB and WB, respectively, in the subplate I of Fig. 3. We now consider vertices of these three subplates by following the above results obtained. Due to the
Remark 4.1: We note that (4.15) implies and , which results in the backward waves and by following (3.7). Resulting from (3.12), we thus . We further consider the voltages of the gyrator obtain coupling between the loops with currents and . Due to the
(4.14)
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
condition at the voltage
1655
, we have a short circuit that occurs , i.e.,
Following from the discussion above, it is clear that by comwith the additional conbining the existing condition , the lower gyrator coupling between the corredition , sponding series adaptors vanishes, i.e., 4. We now consider another possible set of additional conditions being made for determining vertices of the plate: which correspond to the upper gyrator coupling between the loops with currents and . Taking sums and differences of (4.7) and (4.9), we immediately obtain
(4.18)
Based on the results obtained in Remark 4.1, equations in (4.18) have resulted in dropping the entire upper gyrator coupling with from the MDWDF circuit in Fig. 2, i.e., coefficient , 2. Since the MDWDF circuit is designed in a symmetric form, such additional conditions resulting in the removal of either gyrator coupling simply yield the same results. 2) Hard-Type Simply-Supported Edges (HS-HS-HSHS): Combining the additional condition given by (4.12) with the boundary conditions given by (4.6), (4.9), and (4.11), respectively, it is seen that the following state outputs outside each straight edge are readily is added to EB obtained. In particular, the condition and WB of the subplate III in Fig. 3 for obtaining the state and outputs
which re(4.11) with the additional conditions sults in (4.12) and (4.18), we simply obtain vertices of the rectangular plate in terms of state outputs as follows:
(4.21)
(4.22)
(4.23)
(4.24) 3) Soft-Type Simply-Supported Edges (SS-SS-SS-SS): The given boundary conditions resulting from (4.6), (4.11)–(4.12) simply yield the following state outputs outside boundaries of each subplate shown in Fig. 3:
(4.25)
(4.26) (4.19)
Considering two additional conditions that yield (4.15) and wipe out the lower gyrator coupling, and then combining with the existing conditions given by (4.6), and (4.11)–(4.12), the vertices of the rectangular plate are immediately obtained as follows:
(4.20)
Similar to the case of the free-edges boundary conditions for obtaining vertices of these three subplates, the second additional condition is applied to associate with the existing con. This is due to the gyrator coupling between the dition loops with currents and . We thus obtain the state equations (4.18). Combining the existing boundary conditions (4.6) and
(4.27)
(4.28)
1656
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
and (4.32) and (4.33) shown at the bottom of the page, and
(4.29)
(4.30) We note that as mentioned in the case of the free-edges boundary conditions, vertices of the type SS-SS-SS-SS can be also ob. tained by employing the additional conditions This has resulted in the upper gyrator coupling vanished and the resulting vertices are available in (4.21)–(4.24). 4) Clamped Edges (C-C-C-C): The conditions can be set immediately by following the state equations in (4.6), (4.9), and (4.10) so that the state outputs outside each boundary are obtained for these three subplates as follows:
(4.31)
(4.34) To determine vertices of the rectangular plate, two additional conditions either or may be imposed. More specifically, as mentioned in the case of combined with free-edges boundary condition, (respectively, combined with ) has resulted in the termination of the upper gyrator connection (respectively, the lower gyrator connection). Thus, we describe only the verare tices for which additional conditions implemented as shown in (4.35)–(4.38) at the bottom of the page. 5) Mixed Edges: Based on the boundary conditions considered above, clearly different combinations of edge conditions can be achieved in a straightforward manner. In order to make use of the results obtained previously, we confine our attention to the mixed-edges boundary conditions for which two opposite edges have the same type of edge condition. Clearly, there are six such combinations.
(4.32) (4.33)
(4.35)
(4.36)
(4.37)
(4.38)
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
5a)
5b)
5c)
5d)
5e)
5f)
C-F-C-F: Clamped edges on the western and eastern boundaries, while free edges on the northern and southern boundaries. SS-F-SS-F: Soft-type simply-supported edges on the western and eastern boundaries, while free edges on the northern and southern boundaries. F-HS-F-HS: Free edges on the western and eastern boundaries, while hard-type simply-supported edges on the northern and southern boundaries. SS-HS-SS-HS: Soft-type simply-supported edges on the western and eastern boundaries, while hard-type simply-supported edges on the northern and southern boundaries. SS-C-SS-C: Soft-type simply-supported edges on the western and eastern boundaries, while clamped edges on the northern and southern boundaries. C-HS-C-HS: Clamped edges on the western and eastern boundaries, while hard-type simply-supported edges on the northern and southern boundaries.
1657
5c)
5d)
In the following, we describe the state outputs for all six combinations of boundary conditions. 5a)
5b)
C-F-C-F: For this type of mixed boundaries, no additional boundary condition is needed for straight-line edges as well as vertices. As a result, the state outputs outside SB and NB are simply attached by (4.13), while the state outputs outside WB and EB are given by (4.33) and (4.34). The state outputs appeared at each of the rectangular plate are set of vertices also obtained by making use of the state output equations given by (4.6), (4.9), (4.11) and (4.12), and (4.15) as shown in the equation at the bottom of the page. SS-F-SS-F: No additional condition is required to establish this type of mixed boundary for the straight-line edges of the rectangular plate. We thus have the state
5e)
5f)
outputs outside SB and NB that are given by (4.13) and the state outputs outside WB and EB that are given by (4.26). Combining the additional condition with the existing condition which results in the state equations (4.15), the state outputs appeared of the plate are then at each set of vertices given by (4.27) and (4.30). F-HS-F-HS: Imposing the additional condition on WB and EB, the state outputs on their straight-line edges are given by (4.14), while the state outputs on SB and NB are simply attached by (4.19). With the additional condition, four sets of vertices for the plate are also obtained as shown in the first equation at the bottom of the next page. SS-HS-SS-HS: No additional condition is needed to build up such a combination of boundary conditions for the straight-line edges of the plate. Thus, the state outputs outside the straight-line edges are simply attached by (4.19) and (4.26) for (SB,NB) and (WB,EB), respectively. Furthermore, the sets of vertices are also settled by (4.21) and (4.24), provided the adholds. ditional condition SS-C-SS-C: Operating this type of boundary condition does not require any additional condition, as the combination of existing conditions take place at each straight-line edge and vertex. As a result, the state outputs outside the straight-line edges are set immediately, following from the state outputs given by (4.31)–(4.32) for SB and NB and the state outputs given by (4.26) for are also estabWB and EB. The vertices lished as shown in the second equation at the bottom of the next page. C-HS-C-HS: The state outputs outside each straight-line edge can be set immediately by making use of (4.33)–(4.34) and (4.19) directly attached to
1658
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
(WB,EB) and (SB,NB), respectively. Provided the adholds, four sets of vertices ditional condition are obtained as for the case of SS-C-SS-C. Tables I and II briefly summarize the common edges and mixed edges boundary conditions for which additional conditions built on the basic conditions are imposed on the Mindlin model.
V. NUMERICAL RESULTS We present two types of numerical example for simulating the vibrations of isotropic Mindlin plates, implemented by the MDWDF algorithm depicted in Fig. 2 with the choices given by (3.3): isotropic thin plate and isotropic thick of plates. In the first case, the material body involved is an isotropic thin square plate for the purpose of obtaining the
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
1659
p
TABLE I COMMON-EDGE BOUNDARY CONDITIONS FOR MINDLIN PLATES. “ ” : BASIC CONDITION; “ ” : ADDITIONAL CONDITION
2
p
TABLE II MIXED EDGE BOUNDARY CONDITIONS FOR MINDLIN PLATES. “ ” : BASIC CONDITION; “ ” : ADDITIONAL CONDITION
2
TABLE III MINIMAL RATIO AND INDUCTANCES OBTAINED BY THE GOLDEN SECTION SEARCH METHOD
well-known Chladni figures [24]. For the second case, three material bodies concerning isotropic thick plates are presented. is taken to be 5/6, In both cases, the shear correction factor 4 cm, and the material paramthe grid spacing is set to eters used in the algorithm are shown in Table III for either brass or steel. To visibly show the deformation of the plate surface, some simulation results are illustrated in terms of the displacement by taking the integral of the transverse velocity at a period of time . As a result, we have
At
the initial value of displacement is given by for both cases.
A. Isotropic Thin Plate This example uses a thin square plate with the side length 1 m for obtaining Chladni’s figures [24]. The square plate is assumed to be made of brass with sheeting at mm in thickness. The initial equation for this thin plate is given by
quantities outside the plate edges are given by (4.13)–(4.14) and (4.16)–(4.17). Using the golden section search method, the minis obtained that also results in the minimal inimal ratio and . These results are illustrated in Table III. ductances It is worth noting that the minimal ratio obtained has indicated the least restriction on the density of the sampling in time for a given density of the sampling in space according to the CFL given in (3.4), these results shown bound. Combined with in Table III have essentially not only guaranteed the stability of Mindlin system but also resulted in a lossless case of the systems. The nodal patterns obtained for different values of and are shown in Fig. 4. They are very similar to those appearing in [24] which were obtained experimentally. They are known as Chladni figures and the lines indicate parts of the plate that , the have zero displacement. Along the main diagonal, Chladni figures in the upper triangular part are generated using the negative sign of (5.1) whereas the figures in the lower triangular part are generated using the positive sign of the same equation. This graphical arrangement recorded as far as the 7 7 ” and “ ,” respecmode is indicated in Fig. 4 by “ tively. B. Isotropic Thick Plates
for nodal lines where modes variables
(5.1)
is the amplitude at the point ; the number of provided . Since the plate and move toward the same trajectory, we may set when . We note that the initial equation (5.1) is particularly designed to combine with the free edges boundary condition where the state ouput
In this subsection, the implementation of the MDWDF algorithm is used to demonstrate several aspects of thick plate vibrations. These aspects are: plane-wave propagation with free edges, plate deformation with hard-type simply-supported edges, and a combination of plate deformation and plane-wave propagation with mixed edges. 1) Plane-Wave Propagation: This example demonstrates the plane-wave propagation in terms of the transverse velocity
1660
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
Fig. 4. Diagram of Chladni’s figures for elastic square thin plates made of brass with free edges.
by using a square Mindlin plate with side length 1 m. cm The square plate body has a uniform thickness of and is made of steel. The MDWDF algorithm is implemented with the free edges boundary condition and initialized only and with a transverse velocity that is centred at takes the form
The minimal values of the ratio and inductances obtained by the golden section search method are shown in Table III. The resulting wave propagation is illustrated in Fig. 5 where four subsequent time-points are recorded. More specifically, the reflection effect of plane waves across the boundaries is visibly shown in the last two plots of Fig. 5. Similar results to those depicted in Fig. 5 are also to be found in [16] where the digital waveguide network approach was used. 2) Plate Deformation: This example concerns the symmetric plate deformation of an isotropic square plate made of 1 m. To visibly demonstrate brass with side lengths the plate deformation in terms of upward and downward surface displacement, the transverse velocity is initialized by
(5.2) where are positive integers representing the number of directions, respectively, and are signed modes in the
and is a positive scale factor. integers such that . In We note that the plus sign “ ” in (5.2) is only when , , this example, we take and . The boundary conditions involved are set by the hard-type simply-supported edges for which the state outputs outside the plate edges are determined by (4.19)–(4.24). The minimal ratio and inductances for the plate thickness 0.1 m are illustrated in Table III. Plots depicted in Fig. 6 at four 0.1 m. time-points show the plate deformation in the case , , By taking the -norm of the plate variables , , the results of detailed comparison made for the two plate thick0.1 m and 0.2 m with a common ground initial nesses transverse velocity are given in Fig. 7 at different time points for which . In view of Fig. 7, we make the following observations which are commonly obtained through material experiments: i) Due to the energy reflection and transmission from and across the boundaries, the rate of the total energy, sum of the kinetic and strain, for the elastic plate body is gradually decreased with increasing time. This can be shown in Fig. 7(b)-(d) where plate velocities and the plate stress are related to the kinetic and strain energies, respectively. 0.1 m, the ii) As compared to the plate with thickness 0.2 m was able to susthicker plate with thickness tain vibration longer before material failure occurred according to Fig. 7(a). However, the thicker plate generates a larger strain energy of the material body according to Fig. 7(c)-(d). As a result, a significant energy jump can be clearly seen in Fig. 7(d) which indicates that the thicker plate is approaching the status of material failure.
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
Fig. 5.
1661
Plane-wave propagation in terms of the transverse velocity for an elastic thick square plate made of steel with free edges at different time points.
Fig. 6. Symmetric plate deformation for an elastic thick square plate made of brass with hard-type simply-supported edges at different time points.
Observations i) and ii) are consistent with those results obtained in [18]. 3) Combination of Plate Deformation and Plane-Wave Propagation: To demonstrate the effects of a combination of plate deformation and plane-wave propagation, an isotropic rectan2m gular elastic plate body made of steel with side lengths
1 m is used. Boundary conditions are of the mixed and edges SS-F-SS-F where state outputs outside boundaries are given by (4.13), (4.26), and (4.27)–(4.30). The initial value of the transverse velocity is given by
1662
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
Fig. 7. ` -norm quantity comparisons of the surface displacement w , transverse velocity v , and plate stress Q and M for two different plate thicknesses of an 0.1 m; “ :” : h 0.2 m. elastic square plate made of brass. “ ” : h
0
=
=
Fig. 8. Combination of plate deformation and plane-wave propagation for an elastic thick rectangular plate made of steel with mixed edges SS-F-SS-F at different time points.
while the remaining plate variables are set to 0. For the case 0.1 m with the scale paramof the plate thickness set by , the results of minimal ratio and inductances eter shown in Table III, are applied to the MDWDF algorithm for the
optimal performance of maintaining the passivity and stability of the Mindlin system. Considering different time slots , the graphic illustrations shown in Fig. 8 have given a clear picture of the plate deformation and plane-wave propagation including
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
energy reflection and transmission from and across the boundaries in terms of the to and fro motion of the constantly varying displacement. VI. CONCLUSION We have discussed predominantly the implementation of initial and boundary conditions in a MDWDF algorithm for Mindlin plate vibration, described by a fourth-order PDE system. The MDWDF algorithm, implemented by state output quantities, requires only local interconnections for which the behavior at any point in space is directly related to the points in its nearest neighborhood. This property is thus very suited to serve as a basis for the construction of different types of common boundary conditions applied to the area of mechanics. Based on the establishment of these edge conditions, several types of mixed edge boundary conditions encountered in practice can be easily established in a very simple and efficient manner without remodifying the whole algorithm. Two cases with different types of boundary conditions were investigated for demonstrating the capacities of the MDWDF approach: a thin isotropic plate and a thick isotropic plate. We were able to demonstrate, for the thin elastic plate, that the Chladni figures we have produced were similar to those produced experimentally by Waller. Results obtained from the simulations of the thick elastic plate bodies have clearly shown the plate deformation and plane-wave propagation including energy reflection and transmission from and across the boundaries. The fourth-order system we have investigated is of sufficient complexity to enable the effective consideration of a suitable parallel implementation. This is the subject of ongoing research. REFERENCES [1] A. Fettweis, “Digital filters related to classical filter networks,” AEÜ, vol. 25, pp. 79–89, 1971. , “Wave digital filters: Theory and practice,” Proc. IEEE, vol. 74, [2] pp. 270–327, 1986. [3] A. Fettweis and G. Nitsche, “Transformation approach to numerical integrating PDEs by means of WDF principles,” Multidim. Syst. Signal Processing, vol. 2, pp. 127–159, 1991. [4] , “Numerical integration of partial differential equations using principles of multidimensional wave digital filters,” J. VLSI Signal Processing, vol. 3, pp. 7–24, 1991. [5] G. Nitsche, “Numerische Lösung Partieller Differentialgleichungen mit Hilfe von Wellendigitalfiltern,” Ph.D. dissertation, Fakultät Electrotech. Info., Ruhr-Universität, Bochum, Germany, 1993. [6] X. Liu and A. Fettweis, “Multidimensional digital filtering by using parallel algorithms based on diagonal processing,” Multidim. Syst. Signal Processing, vol. 1, pp. 51–66, Mar 1990. [7] A. Fettweis and K. Meerkötter, “On parasitic oscillations in digital filters under looped conditions,” IEEE Trans. Circuits Syst., vol. CAS-24, pp. 475–481, Sept. 1977. [8] C. T. Mullis and R. A. Roberts, “Synthesis of minimum roundoff noise fixed point digital filters,” IEEE Trans. Circuits Syst., vol. CAS-23, pp. 551–562, Sept. 1976. [9] A. Fettweis, “Discrete passive modeling of physical systems described by PDEs,” in Proc. EUSIPCO-92, 6th Eur. Signal Processing Conf., vol. 1, Brussels, Belgium, 1992, pp. 55–62. [10] G. Hemetsberger, “Stability verification of multidimensional Kirchoff circuits by suitable energy functions,” in Proc. IEEE Int. Conf. Acoustics, Speech, Signal Processing, vol. 6, Adelaide, Australia, 1994, pp. 13–16. [11] C. Q. Xu, S. C. Bass, and X. Wang, “Accommodating linear and nonlinear boundary conditions in wave digital simulations of PDE systems,” J. Circuits Syst. Comput., vol. 7, pp. 563–597, 1997.
1663
[12] C. A. J. Fletcher, Computational Techniques for Fluid Dynamics.-1: Fundamental and General Techniques, 2nd ed. New York: SpringerVerlag, 1991. [13] A. Fettweis, H. Levin, and A. Sedlmeyer, “Wave digital lattice filters,” Int. J. Circuit Theory Applicat., vol. 2, pp. 203–211, 1974. [14] H. Krauss, R. Rabenstein, and M. Gerken, “Simulation of wave propagation by multidimensional digital filters,” J. Simulation Practice Theory, vol. 4, pp. 361–382, 1996. [15] C. H. Tseng and S. S. Lawson, “Incorporating initial and boundary conditions into the multidimensional wave digital filter,” in Proc. IEEE Int. Conf. Acoustics, Speech, Signal Processing, Orlando, FL, May 2002. [16] S. D. Bilbao, “Wave and scattering methods for the numerical integration of partial differential equations,” Ph.D dissertation, Dept. Elect. Eng., Stanford Univ., Stanford, CA, 2001. [17] Computer Solutions Europe AB, Partial Differential Equation Toolbox, The Math Works, Natick, MA, 1996. [18] K. F. Graff, Wave Motion in Elastic Solids. New York: Dover, 1975. [19] R. D. Mindlin, “Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates,” Trans. ASME J. Appl. Mech., vol. 18, pp. 31–38, 1951. [20] K. Liew, C. Wang, Y. Xiang, and S. Kitipornchai, Vibration of Mindlin Plates: Programming the p-Version Ritz Method, first ed. Amsterdam, The Netherlands: Elsevier, 1998. [21] E. Hinton and B. S. Al-Janabi, Numerical Methods and Software for Dynamic Analysis of Plates and Shells, E. Hinton, Ed. Swansea, U.K.: Pineridge Press, 1987. [22] H. Huang, Static and Dynamic Analyzes of Plates and Shells. London, U.K.: Springer-Verlag, 1989. [23] S. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells, 2nd ed. New York: McGraw-Hill, 1959. [24] M. D. Waller, Chladni Figures: A Study in Symmetry. London, U.K.: G. Bell and Sons, 1961. [25] R. Nouta, “The use of Jaumann structures in wave digital filters,” Int. J. Circuit Theory Applicat., vol. 2, pp. 163–174, 1974. [26] A. A. Samarskij, Theorie der Differenzenverfahren (in German). Leipzig: Akademicshe Verlagsgesellschaft, 1984.
Chien-Hsun Tseng was born in Taiwan in 1968. He received the B.Sc. degree in applied mathematics from National Cheng Kung University, Tainan, Taiwan, the M.Math. degree from University of New South Wales, Sydney, Australia, and the Ph.D. degree in electrical engineering from Curtin University of Technology, Perth, Australia, in 1995, 1997, and 2000, respectively. From 2000 to 2001, he was a Research Scientist at Weierstrass Institute for Applied Analysis and Stochastics (WIAS), Berlin, Germany. Since June 2001, he has been with the University of Warwick, Warwick, U.K., where he is a Research Fellow in the School of Engineering. His current research interests include optimal digital and statistical filters design through optimization techniques, modeling of physical systems using digital signal processing techniques, and stochastic and numerical optimizations.
Stuart Lawson (M’80–SM’91) received the B.Sc. degree in mathematics (with first class honors) from Woolwich Polytechnic and the D.I.C., U.K., and the Ph.D. degree in electrical engineering from Imperial College, University of London, London, U.K., in 1971 and 1975, respectively. From 1971 to 1978, he was a Research Engineer with Plessey Telecommunications, U.K. He taught and researched at City University, London, U.K., from 1978 to 1987. Since 1988, he has been at the University of Warwick where he is now a Reader in the School of Engineering. His current research interests include digital filters and filter banks, image compression, modeling physical systems using digital signal processing techniques and sonar signal processing. Dr. Lawson is a Fellow of Institution Electrical Engineers and a member of the Advisory Panel for the IEE Signal Processing Professional Network.
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
Initial and Boundary Conditions in Multidimensional Wave Digital Filter Algorithms for Plate Vibration Chien-Hsun Tseng and Stuart Lawson, Senior Member, IEEE
Abstract—Recently, multidimensional wave digital filter (MDWDF) structures have been proposed for the modeling of plate vibration problems. In this paper, we discuss how initial and boundary conditions may be properly embedded into such an algorithm in terms of the state quantities that are an integral part of the algorithm. Due to the essential feature of fully-local interconnectivity in the MDWDF model, different types of boundary conditions can be easily satisfied in a very simple and efficient manner. Instead of remodifying the whole algorithm, usually required by finite elements based methods, boundary conditions in terms of state outputs are simply attached to the model. This feature is especially useful when dealing with the mixed-edges boundary conditions frequently encountered in practice. Graphical results obtained from implementing the MDWDF algorithm are given to further demonstrate the capacities of the method in efficiently handling a fourth-order Mindlin plate vibration system with various types of boundary conditions. Index Terms—Courant–Friedrichs–Levy bound, golden section search method, Kirchhoff circuit, losslessness, Mindlin system, multidimensional wave digital filters (MDWDF), partial differential equations (PDEs), passivity, stability, trapezoidal rule.
I. INTRODUCTION
T
HE modeling of plate vibration in mechanics is of great importance to applications in engineering and science, especially in the areas of aerospace, marine and construction. The vibration aspects in micorelectromechanical systems applications such as electronic packaging, micro-robots, etc. are increasingly important because of their enhanced sensitivities to vibration. The most popular kind of models for mechanical vibration systems can be represented by sets of linear and/or nonlinear partial differential equations (PDEs) with properly imposed initial and boundary conditions. In most cases of complex vibration-influenced design, the equations can not be solved analytically and so numerical approaches are sought. As far as practical interest is concerned, it is, therefore, desirable to numerically simulate the behavior of such PDE systems as cheaply and rapidly as possible, especially when dealing with complicated boundary conditions. Numerical techniques used to solve the PDEs such as finite elements and finite differences have had some success. The finite-element methods permit easy inclusion of local grid refinement and handling of complex geometries. These methods, Manuscript received August 12, 2002; revised May 18, 2003. This work was supported by EPSRC(UK) under Grant GR/N22298. This paper was recommended by Associate Editor M. Nakhla. The authors are with the School of Engineering, University of Warwick, Coventry CV4 7AL, U.K. (e-mail: [email protected]; [email protected]; [email protected]) Digital Object Identifier 10.1109/TCSI.2004.832796
however, are computationally expensive and harder to correctly set up the simulation plane than the frequently used finite difference methods. The latter approaches, on the other hand, have difficulties in handling irregular boundaries, and need local grid refinement to increase the accuracy. A remarkable alternative approach to the modeling and simulation of a system of PDEs builds on the wave digital filter (WDF) paradigm and analogies with electrical networks, an area rich in mathematical structure [1], [2]. By using such a paradigm and properties of electrical networks, concepts well known to circuit theorists such as passivity, losslessness, stability, finite propagation velocity, etc., can be translated to the actually considered physical problems so as to preserve important relationships between variables. The technique involves first finding a multidimensional lumped electrical network which represents the behavior of the system. From this network, a discrete-time equivalent is developed that is a multidimensional WDF (MDWDF) [3]–[5]. Proceeding in this way, one can substantially achieve for the resulting algorithm: • massive parallelism because each point in the -D grid can be updated simultaneously if sufficient computing resources are available [3], [6]; • high accuracy because of making use of the WDF structure known to have low roundoff noise characteristics and suppression of parasitic oscillations [2], [7], [8]; • multidimensional passivity (reflecting the passivity of the system to be modeled) that guarantees robustness and basically all numerical stability properties required of an accurate numerical integration method [9], [10]; • full local interconnectivity due to the generation of a second-order difference equation for which its behavior at any point in space is directly influenced only by the points in its nearest neighborhood [4], [11]. This local interconnectivity is particularly helpful to the establishment of the boundary conditions in a very simple manner. Other second-order methods do not offer this property [12]. The technique of MDWDF originally introduced for the implementation of digital filters [1], [2], [13] has been successfully applied to the numerical integration of PDEs [3], [5], [14]–[16], etc. These PDEs systems usually involve low-order theories for the governing equations. As a result, the systems can be directly handled by a standard PDE solver such as the MATLAB PDE Toolbox [17]. In this paper we discuss in particular the crucial issue of the setting of initial and boundary conditions with regard to a recently proposed MDWDF algorithm which simulates plate vibration. We begin by considering the thick plate theory due to Mindlin [18]–[20], the first approximate theory for dealing especially
1057-7122/04$20.00 © 2004 IEEE
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
1649
with flexural waves in relatively thick plates, although thin plates can also be dealt with. We note that a plate is typically considered to be thin if the ratio of its thickness to representative lateral dimension (e.g., circular plate diameter, square plate side length) is 1/20 or less. Modeling of the Mindlin system had been achieved using many different techniques such as the finite element, the boundary element, the finite difference, the differential quadrature, the collocation, the Galerkin and the Ritz methods [21], [22]. The idea of using a digital network to model Mindlin plate vibration was first reported in [16] where the networks used were based on: 1) digital wave guides (DWN) and 2) MDWDF using the space transformation method [3]. Good results were obtained using the DWN model only, the MDWDF model did not yield results due to problems in the formulation of the initial and boundary conditions [16]. This paper begins with the modeling of the Mindlin system in detail. In Sections II and II, we present a variant of the model discussed by Bilbao [16]. Following a direct space–time analysis methodology [4], we start from the system of PDEs and derive the equivalent analogue network description. This results in a MDWDF algorithm. A useful feature of MDWDF modeling is the ability at this stage to satisfy the system’s stability requirement in order to ensure the passivity of the circuit obtained. The next step is to properly incorporate the given initial and boundary conditions into the open-field MDWDF algorithm obtained. It turns out that due to the advantage of fully local interconnectivity of the MDWDF approach, a straightforward and efficient way to accommodate various boundary conditions without remodifying the whole MDWDF algorithm, is provided. Other finite element-based methods do not have this advantage. Finally, two cases of plate vibration are considered to demonstrate capacities of the method in handling the high-order governing equations. II. PROBLEM STATEMENT The stress equations of motion of a two-dimensional (2-D) Mindlin plate can be written as a system of PDEs [18], [20]
(2.1)
where is the plate thickness, are the bending rotations of a transverse normal to the and axes, respectively, are the transverse shear forces per unit length of the are the bending moments per unit plate, and length of plate. The plate material is assumed to have the density , Young’s modulus/modulus of elasticity in tension and com. pression , as well as Poisson’s ratio , such that Assuming a plane stress distribution in accordance with Hooke’s law, the plate stress-displacement relations are then given by
(2.2) is a shear correction factor being applied to compenwhere sate for an error in assuming a constant shear stress through is the modulus of elasthe plate thickness, is the flexural rigidity ticity in shear, and of the plate. When the stress-displacement equations in (2.2) are incorporated in the system defined by (2.1), we obtain the governing equations of motion described as a system of eight , , and PDEs as follows by defining : (2.3)
(2.4)
In view of (2.3) and (2.4), it was noted in [16] that the complete Mindlin system can be decomposed into two separate subsystems where the corresponding system variables are written and . Furthermore, the as differentiable equations of motion can be reduced to a single equation for the displacement by firstly substituting the plate stress-displacement equations of (2.2) into (2.1). Followed up by differentiating these equations with respect to and then summing, one finally obtains a fourth-order Mindlin plate vibration system in time and space
(2.5) In order to apply the circuit theory for representing the behavior of the system in (2.3)–(2.4), the same physical dimension must be kept for all plate variables. As a result, scaled dependent combined with positive constants variables , 2, 3 are introduced as (2.6) where the transverse velocity , and velocities of transverse and are treated as the graphical voltage bending rotation denoted by corresponding to the graphical current denoted by . It was noted in [16] that the plate system in (2.3)–(2.4) can be replaced by the following form written as symmetric hyperbolic PDEs: (2.7)
1650
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
(2.8)
where the operators , denote the partial derivatives with respect to temporal and spatial coordinates in and axes. III. GRAPHICAL NETWORK DESCRIPTION AND MULTIDIMENSIONAL WAVE DIGITAL FILTERING ALGORITHM In this section, the derivation of a three-dimensional (3-D) Kirchhoff circuit is presented. Moreover, the linear discretization using a generalized trapezoidal rule is applied to each component of the Kirchhoff circuit for achieving a structure that can be numerically implemented in terms of difference equations.
Fig. 1. Multidimensional Kirchhoff circuit representation for Mindlin systems (2.3)–(2.4).
with respect In addition, the partial derivative operators to the two-port Jaumann structure [25] are given by defining arbitrary positive constants
A. Network Description The key success in the derivation of a MDWDF algorithm for the numerical solution of Mindlin system is the derivation of a MD Kirchhoff circuit from the underlying PDEs. By representing these equations (2.7)–(2.8) in graphical form, one can express mathematical manipulations as network operations known as the Kirchhoff circuit theory in the following:
Remark 3.1: In view of Fig. 1, we note that by applying the techniques of circuit theory, inductors of a MD network description resulted from the plate stress-displacement PDEs in (2.3)–(2.4) can not be realized by real-world components. Instead, the network description is a graphical representation of time and space differentiation operators for the system of (2.3)–(2.4). Furthermore, the graphical circuit representation in Fig. 1 has resulted in a convenient way of understanding the passivity concept and grid generation for MD physical systems [4], [9], [10]. B. MDWDF Algorithm
where , , and . Following the Kirchhoff’s laws, the above equations can be immediately interpreted in the form of a 3D Kirchhoff circuit represented in Fig. 1, see also [16, Fig. 5.16]. In view of the MDKC of Fig. 1, the conservation of energy holds as a result of Kirchhoff’s laws. Furthermore, the graphical circuit contains three main loops on the left, five main loops on the right, and two gyrators coupling in the middle with the gyration coefficients , and . The self-inductors with inductances are given by
To arrive at an approximation of each mesh voltage across the corresponding inductance in the 3D Kirchhoff circuit of Fig. 1, a generalized trapezoidal rule [3], [4] is applied to the integration of each mesh voltage involved in two dimensions where , take only discrete values with the temporal step size and the spatial and in and coordinates, respectively. As a step sizes, result, we obtain the port resistances and for the two-port Jaumann filters and self-inductors, respectively, given by
(3.2)
When the condition for the passivity of inductances is held in the MDKC of Fig. 1, i.e., , clearly the simplest choice for these is given by
(3.1)
(3.3)
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
1651
Fig. 2. MDWDF algorithm for numerical simulation of Mindlin system (2.3)–(2.4).
Substituting in
chosen by (3.3) into (3.1)–(3.2), it simply results , , , 3, 6, 7, and
now transform the MDKC of Fig. 1 into the MDWDF algorithm depicted in Fig. 2 (see also [16, Fig. 5.16]) where the shift operators are obtained by
(3.4)
where the plus sign “ ” occurs when . Following from (3.2)–(3.3), we are able to determine the resistances of port adaptors in Fig. 2 as shown in (3.5) at the bottom of the page. To guarantee that the gyrator couplings between the loops with , and in the MDKC of Fig. 1 are memorycurrents less, reflection-free ports are necessary [16] which consequently determined by leads to the port resistances
Defining and following from 1, 2, 3 in (3.1), we obtain a lower bound of a ratio between the density of sampling in space and that of sampling where is the in time, i.e., phase velocity of a wave [18] traveling freely across the whole plate according to the fourth-order system (2.5). This is the well-known Courant–Friedrichs–Levy (CFL) bound [26], a stability criterion for many approaches to the numerical solution of the wave equation. For a simplified reference circuit where , the are determined by (3.3) and hence resulting in the minimal ratio can be obtained by the golden section search method. Thus, it imposes the least restriction on the density of the sampling in time for a given density of the sampling in space so that the numerical stability of the Mindlin system can be guaranteed when simulated by a passive MDWDF circuit. and are minimized to Furthermore, the inductances achieve the losslessness of the system. Assuming that the voltage wave representation is adopted and making use of properties of the two-port lattice structure which is equivalent to the two-port Jaumann structure [25], we may
As a result, a delay-free directed loop is not created leading from and the input to the output terminal. Let be the currents corresponding to the wave quantities and , respectively, in the MDWDF algorithm of Fig. 2. By synthesizing various ports in the MDKC of Fig. 1, using the series form and given below in terms of currents
(3.6)
(3.5)
1652
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
the wave quantities of these series adaptors in the MDWDF of Fig. 2 can be described as
Referring back to the MDWDF algorithm of Fig. 2, the relations of wave and state quantities can be seen as
for
for
(3.7)
Of particular interest are the reflection-free ports in the gyrator couplings. According to (3.7), the then become independent . Let the vector of voltage waves [5] be of for and where the vector given by , and the vector voltage and current are defined waves by , , , and
(3.8) Thus, the coupled two-port series adaptor described by vector waves is obtained by (3.9) (3.10) where the vector waves and , and are and the coupling matrix of the port resistances given in (3.5). It is to be noted that vector wave variables were first used by Nitsche in his development of a MDWDF model for Timoshenko’s equations for beam bending [5]. They have also been used for the Mindlin plate model in [16]. Finally, for , 2 the two-port series adaptors (in the middle of Fig. 2) due to the gyrator couplings are given by (3.11) (3.12) where
(3.20) Furthermore, the state input-output relations are given by
(3.21) where the plus sign “ ” occurs when and and outputs.
, are vectors of state
IV. INITIAL AND BOUNDARY CONDITIONS In this section, we discuss how to accommodate the properly imposed initial and boundary conditions into an open-field MDWDF algorithm depicted in Fig. 2 in terms of state equations given by (3.21). A. Initial Conditions Let initial conditions of the Mindlin plate variables given in (2.6) be imposed as follows:
(4.1)
where , and . We first consider the first derivative initial conditions . Replaced the plate variables in (3.13)–(3.19) by their initial values in (4.1) and following from the results obacross each selftained in [15], we have the forward waves in the MDKC of Fig. 1: inductor with inductance
are given by
We further consider the vector of voltage a coupled inductor with inductance series connections of and , we have
(4.2) given in (3.8) across . Due to , by defining
The desired plate variables following from (3.5)–(3.12) can now and (indices be expressed in the form of forward waves are omitted) as (3.13) (3.14) (3.15) (3.16) (3.17) (3.18) (3.19)
Replaced by (3.10) and substituting by (3.18) with their , 7, it holds that initial values (4.3) Taking sums and differences of the plate variables replaced by their initial vales in (3.14)–(3.17), and then combining together given in (4.2), with (3.13), (3.19), (4.3), and
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
Fig. 3.
1653
Decomposed rectangular Mindlin plates with boundaries calculated by state output on edges and vertices.
we obtain two sets of state output equations, respectively, corresponding to initial conditions of the plate variables in two sepand (inarate subsystems dices and are omitted)
(4.4)
(4.5)
Let , , , and , where are given function of . Solving (4.4)–(4.5), we thus obtain the remaining state outputs shown in the equation at the bottom of the page. B. Boundary Conditions There are four common types of boundary conditions applied to the plate edge of Mindlin plates [23], [22], [20]: Free edge Hard-type simply-supported edge Soft-type simply-supported edge Clamped edge.
Considering a rectangular Mindlin plate whose grid points are calculated by state outputs (3.21) which involve both the spatial shifts and time delay, it seems clear that the plate can be decomposed into three separate subplates depicted in Fig. 3. Starting with the straight-line plate boundaries parallel with either the or axis illustrated in Fig. 3, clearly one of four neighborhood points is outside the boundary of each subplate. This has resulted in four combinations of state in, , , and , puts respectively, corresponding to the southern boundary (SB), western boundary (WB), northern boundary (NB), and eastern is one spatial step boundary (EB) where each state input size beyond the boundary of the rectangular plate. Instead of being obliged to determine sets of state inputs which are not computable, our aim is to find expressions for the corresponding set of state outputs that accommodate the given boundary conditions. Consequently, a state is enough for each equation determining each state output straight-line edge of the plate. This, however, will require us to impose at most one additional boundary condition. Similar to the straight-line edge grid points, two of four grid for each subplate. As a result, points are outside the vertex at most two additional boundary conditions have to be imposed for providing vertices of these subplates. Thus, each set of veris determined by six state output equations. tices By assembling these three subplates into one complete Mindlin plate, we, therefore, can obtain vertices of the rectangular plate , , , , which are included in the directions of SB, WB, NB, and EB, respectively. In the following, we will discuss how to find these additional conditions that associate with various boundary conditions given above. Substituting (3.20) into the plate variables described by the state equations in (3.13)–(3.19) and replaced the given boundary values, we then obtain the following state
1654
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
output equations directly corresponding to various boundary conditions
that, gyrator coupling between the loops with currents and and , taking respectively, determine the plate variables sums and differences of (4.8) and (4.10) yields
(4.6) (4.15) (4.7)
(4.8)
(4.9)
Combining the boundary conditions (4.11) and (4.12) with two additional conditions which result in (4.7) and (4.15), one is able to assemble these three subplates into one complete Mindlin plate so that the state outputs appearing at for the plate are obtained as foleach set of vertices lows:
(4.10) (4.11) (4.12) We note that (4.11) is a result of mutual inductors occurring in a coupled two-port series adaptor that is described by the plate and in the MDWDF algorithm of Fig. 2. variables 1) Free Edges (F-F-F-F): Following the conditions from (4.8), and (4.11)–(4.12), the state outputs outside the straight-line edge exclusive the vertices are obtained for each subplate by providing the additional conditions as follows:
(4.16)
(4.13) (4.17) and as shown in (4.14) at the bottom of the page. It is worth noting that the additional condition is given to provide the state outputs and for EB and WB, respectively, in the subplate I of Fig. 3. We now consider vertices of these three subplates by following the above results obtained. Due to the
Remark 4.1: We note that (4.15) implies and , which results in the backward waves and by following (3.7). Resulting from (3.12), we thus . We further consider the voltages of the gyrator obtain coupling between the loops with currents and . Due to the
(4.14)
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
condition at the voltage
1655
, we have a short circuit that occurs , i.e.,
Following from the discussion above, it is clear that by comwith the additional conbining the existing condition , the lower gyrator coupling between the corredition , sponding series adaptors vanishes, i.e., 4. We now consider another possible set of additional conditions being made for determining vertices of the plate: which correspond to the upper gyrator coupling between the loops with currents and . Taking sums and differences of (4.7) and (4.9), we immediately obtain
(4.18)
Based on the results obtained in Remark 4.1, equations in (4.18) have resulted in dropping the entire upper gyrator coupling with from the MDWDF circuit in Fig. 2, i.e., coefficient , 2. Since the MDWDF circuit is designed in a symmetric form, such additional conditions resulting in the removal of either gyrator coupling simply yield the same results. 2) Hard-Type Simply-Supported Edges (HS-HS-HSHS): Combining the additional condition given by (4.12) with the boundary conditions given by (4.6), (4.9), and (4.11), respectively, it is seen that the following state outputs outside each straight edge are readily is added to EB obtained. In particular, the condition and WB of the subplate III in Fig. 3 for obtaining the state and outputs
which re(4.11) with the additional conditions sults in (4.12) and (4.18), we simply obtain vertices of the rectangular plate in terms of state outputs as follows:
(4.21)
(4.22)
(4.23)
(4.24) 3) Soft-Type Simply-Supported Edges (SS-SS-SS-SS): The given boundary conditions resulting from (4.6), (4.11)–(4.12) simply yield the following state outputs outside boundaries of each subplate shown in Fig. 3:
(4.25)
(4.26) (4.19)
Considering two additional conditions that yield (4.15) and wipe out the lower gyrator coupling, and then combining with the existing conditions given by (4.6), and (4.11)–(4.12), the vertices of the rectangular plate are immediately obtained as follows:
(4.20)
Similar to the case of the free-edges boundary conditions for obtaining vertices of these three subplates, the second additional condition is applied to associate with the existing con. This is due to the gyrator coupling between the dition loops with currents and . We thus obtain the state equations (4.18). Combining the existing boundary conditions (4.6) and
(4.27)
(4.28)
1656
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
and (4.32) and (4.33) shown at the bottom of the page, and
(4.29)
(4.30) We note that as mentioned in the case of the free-edges boundary conditions, vertices of the type SS-SS-SS-SS can be also ob. tained by employing the additional conditions This has resulted in the upper gyrator coupling vanished and the resulting vertices are available in (4.21)–(4.24). 4) Clamped Edges (C-C-C-C): The conditions can be set immediately by following the state equations in (4.6), (4.9), and (4.10) so that the state outputs outside each boundary are obtained for these three subplates as follows:
(4.31)
(4.34) To determine vertices of the rectangular plate, two additional conditions either or may be imposed. More specifically, as mentioned in the case of combined with free-edges boundary condition, (respectively, combined with ) has resulted in the termination of the upper gyrator connection (respectively, the lower gyrator connection). Thus, we describe only the verare tices for which additional conditions implemented as shown in (4.35)–(4.38) at the bottom of the page. 5) Mixed Edges: Based on the boundary conditions considered above, clearly different combinations of edge conditions can be achieved in a straightforward manner. In order to make use of the results obtained previously, we confine our attention to the mixed-edges boundary conditions for which two opposite edges have the same type of edge condition. Clearly, there are six such combinations.
(4.32) (4.33)
(4.35)
(4.36)
(4.37)
(4.38)
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
5a)
5b)
5c)
5d)
5e)
5f)
C-F-C-F: Clamped edges on the western and eastern boundaries, while free edges on the northern and southern boundaries. SS-F-SS-F: Soft-type simply-supported edges on the western and eastern boundaries, while free edges on the northern and southern boundaries. F-HS-F-HS: Free edges on the western and eastern boundaries, while hard-type simply-supported edges on the northern and southern boundaries. SS-HS-SS-HS: Soft-type simply-supported edges on the western and eastern boundaries, while hard-type simply-supported edges on the northern and southern boundaries. SS-C-SS-C: Soft-type simply-supported edges on the western and eastern boundaries, while clamped edges on the northern and southern boundaries. C-HS-C-HS: Clamped edges on the western and eastern boundaries, while hard-type simply-supported edges on the northern and southern boundaries.
1657
5c)
5d)
In the following, we describe the state outputs for all six combinations of boundary conditions. 5a)
5b)
C-F-C-F: For this type of mixed boundaries, no additional boundary condition is needed for straight-line edges as well as vertices. As a result, the state outputs outside SB and NB are simply attached by (4.13), while the state outputs outside WB and EB are given by (4.33) and (4.34). The state outputs appeared at each of the rectangular plate are set of vertices also obtained by making use of the state output equations given by (4.6), (4.9), (4.11) and (4.12), and (4.15) as shown in the equation at the bottom of the page. SS-F-SS-F: No additional condition is required to establish this type of mixed boundary for the straight-line edges of the rectangular plate. We thus have the state
5e)
5f)
outputs outside SB and NB that are given by (4.13) and the state outputs outside WB and EB that are given by (4.26). Combining the additional condition with the existing condition which results in the state equations (4.15), the state outputs appeared of the plate are then at each set of vertices given by (4.27) and (4.30). F-HS-F-HS: Imposing the additional condition on WB and EB, the state outputs on their straight-line edges are given by (4.14), while the state outputs on SB and NB are simply attached by (4.19). With the additional condition, four sets of vertices for the plate are also obtained as shown in the first equation at the bottom of the next page. SS-HS-SS-HS: No additional condition is needed to build up such a combination of boundary conditions for the straight-line edges of the plate. Thus, the state outputs outside the straight-line edges are simply attached by (4.19) and (4.26) for (SB,NB) and (WB,EB), respectively. Furthermore, the sets of vertices are also settled by (4.21) and (4.24), provided the adholds. ditional condition SS-C-SS-C: Operating this type of boundary condition does not require any additional condition, as the combination of existing conditions take place at each straight-line edge and vertex. As a result, the state outputs outside the straight-line edges are set immediately, following from the state outputs given by (4.31)–(4.32) for SB and NB and the state outputs given by (4.26) for are also estabWB and EB. The vertices lished as shown in the second equation at the bottom of the next page. C-HS-C-HS: The state outputs outside each straight-line edge can be set immediately by making use of (4.33)–(4.34) and (4.19) directly attached to
1658
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
(WB,EB) and (SB,NB), respectively. Provided the adholds, four sets of vertices ditional condition are obtained as for the case of SS-C-SS-C. Tables I and II briefly summarize the common edges and mixed edges boundary conditions for which additional conditions built on the basic conditions are imposed on the Mindlin model.
V. NUMERICAL RESULTS We present two types of numerical example for simulating the vibrations of isotropic Mindlin plates, implemented by the MDWDF algorithm depicted in Fig. 2 with the choices given by (3.3): isotropic thin plate and isotropic thick of plates. In the first case, the material body involved is an isotropic thin square plate for the purpose of obtaining the
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
1659
p
TABLE I COMMON-EDGE BOUNDARY CONDITIONS FOR MINDLIN PLATES. “ ” : BASIC CONDITION; “ ” : ADDITIONAL CONDITION
2
p
TABLE II MIXED EDGE BOUNDARY CONDITIONS FOR MINDLIN PLATES. “ ” : BASIC CONDITION; “ ” : ADDITIONAL CONDITION
2
TABLE III MINIMAL RATIO AND INDUCTANCES OBTAINED BY THE GOLDEN SECTION SEARCH METHOD
well-known Chladni figures [24]. For the second case, three material bodies concerning isotropic thick plates are presented. is taken to be 5/6, In both cases, the shear correction factor 4 cm, and the material paramthe grid spacing is set to eters used in the algorithm are shown in Table III for either brass or steel. To visibly show the deformation of the plate surface, some simulation results are illustrated in terms of the displacement by taking the integral of the transverse velocity at a period of time . As a result, we have
At
the initial value of displacement is given by for both cases.
A. Isotropic Thin Plate This example uses a thin square plate with the side length 1 m for obtaining Chladni’s figures [24]. The square plate is assumed to be made of brass with sheeting at mm in thickness. The initial equation for this thin plate is given by
quantities outside the plate edges are given by (4.13)–(4.14) and (4.16)–(4.17). Using the golden section search method, the minis obtained that also results in the minimal inimal ratio and . These results are illustrated in Table III. ductances It is worth noting that the minimal ratio obtained has indicated the least restriction on the density of the sampling in time for a given density of the sampling in space according to the CFL given in (3.4), these results shown bound. Combined with in Table III have essentially not only guaranteed the stability of Mindlin system but also resulted in a lossless case of the systems. The nodal patterns obtained for different values of and are shown in Fig. 4. They are very similar to those appearing in [24] which were obtained experimentally. They are known as Chladni figures and the lines indicate parts of the plate that , the have zero displacement. Along the main diagonal, Chladni figures in the upper triangular part are generated using the negative sign of (5.1) whereas the figures in the lower triangular part are generated using the positive sign of the same equation. This graphical arrangement recorded as far as the 7 7 ” and “ ,” respecmode is indicated in Fig. 4 by “ tively. B. Isotropic Thick Plates
for nodal lines where modes variables
(5.1)
is the amplitude at the point ; the number of provided . Since the plate and move toward the same trajectory, we may set when . We note that the initial equation (5.1) is particularly designed to combine with the free edges boundary condition where the state ouput
In this subsection, the implementation of the MDWDF algorithm is used to demonstrate several aspects of thick plate vibrations. These aspects are: plane-wave propagation with free edges, plate deformation with hard-type simply-supported edges, and a combination of plate deformation and plane-wave propagation with mixed edges. 1) Plane-Wave Propagation: This example demonstrates the plane-wave propagation in terms of the transverse velocity
1660
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
Fig. 4. Diagram of Chladni’s figures for elastic square thin plates made of brass with free edges.
by using a square Mindlin plate with side length 1 m. cm The square plate body has a uniform thickness of and is made of steel. The MDWDF algorithm is implemented with the free edges boundary condition and initialized only and with a transverse velocity that is centred at takes the form
The minimal values of the ratio and inductances obtained by the golden section search method are shown in Table III. The resulting wave propagation is illustrated in Fig. 5 where four subsequent time-points are recorded. More specifically, the reflection effect of plane waves across the boundaries is visibly shown in the last two plots of Fig. 5. Similar results to those depicted in Fig. 5 are also to be found in [16] where the digital waveguide network approach was used. 2) Plate Deformation: This example concerns the symmetric plate deformation of an isotropic square plate made of 1 m. To visibly demonstrate brass with side lengths the plate deformation in terms of upward and downward surface displacement, the transverse velocity is initialized by
(5.2) where are positive integers representing the number of directions, respectively, and are signed modes in the
and is a positive scale factor. integers such that . In We note that the plus sign “ ” in (5.2) is only when , , this example, we take and . The boundary conditions involved are set by the hard-type simply-supported edges for which the state outputs outside the plate edges are determined by (4.19)–(4.24). The minimal ratio and inductances for the plate thickness 0.1 m are illustrated in Table III. Plots depicted in Fig. 6 at four 0.1 m. time-points show the plate deformation in the case , , By taking the -norm of the plate variables , , the results of detailed comparison made for the two plate thick0.1 m and 0.2 m with a common ground initial nesses transverse velocity are given in Fig. 7 at different time points for which . In view of Fig. 7, we make the following observations which are commonly obtained through material experiments: i) Due to the energy reflection and transmission from and across the boundaries, the rate of the total energy, sum of the kinetic and strain, for the elastic plate body is gradually decreased with increasing time. This can be shown in Fig. 7(b)-(d) where plate velocities and the plate stress are related to the kinetic and strain energies, respectively. 0.1 m, the ii) As compared to the plate with thickness 0.2 m was able to susthicker plate with thickness tain vibration longer before material failure occurred according to Fig. 7(a). However, the thicker plate generates a larger strain energy of the material body according to Fig. 7(c)-(d). As a result, a significant energy jump can be clearly seen in Fig. 7(d) which indicates that the thicker plate is approaching the status of material failure.
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
Fig. 5.
1661
Plane-wave propagation in terms of the transverse velocity for an elastic thick square plate made of steel with free edges at different time points.
Fig. 6. Symmetric plate deformation for an elastic thick square plate made of brass with hard-type simply-supported edges at different time points.
Observations i) and ii) are consistent with those results obtained in [18]. 3) Combination of Plate Deformation and Plane-Wave Propagation: To demonstrate the effects of a combination of plate deformation and plane-wave propagation, an isotropic rectan2m gular elastic plate body made of steel with side lengths
1 m is used. Boundary conditions are of the mixed and edges SS-F-SS-F where state outputs outside boundaries are given by (4.13), (4.26), and (4.27)–(4.30). The initial value of the transverse velocity is given by
1662
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 51, NO. 8, AUGUST 2004
Fig. 7. ` -norm quantity comparisons of the surface displacement w , transverse velocity v , and plate stress Q and M for two different plate thicknesses of an 0.1 m; “ :” : h 0.2 m. elastic square plate made of brass. “ ” : h
0
=
=
Fig. 8. Combination of plate deformation and plane-wave propagation for an elastic thick rectangular plate made of steel with mixed edges SS-F-SS-F at different time points.
while the remaining plate variables are set to 0. For the case 0.1 m with the scale paramof the plate thickness set by , the results of minimal ratio and inductances eter shown in Table III, are applied to the MDWDF algorithm for the
optimal performance of maintaining the passivity and stability of the Mindlin system. Considering different time slots , the graphic illustrations shown in Fig. 8 have given a clear picture of the plate deformation and plane-wave propagation including
TSENG AND LAWSON: INITIAL AND BOUNDARY CONDITIONS
energy reflection and transmission from and across the boundaries in terms of the to and fro motion of the constantly varying displacement. VI. CONCLUSION We have discussed predominantly the implementation of initial and boundary conditions in a MDWDF algorithm for Mindlin plate vibration, described by a fourth-order PDE system. The MDWDF algorithm, implemented by state output quantities, requires only local interconnections for which the behavior at any point in space is directly related to the points in its nearest neighborhood. This property is thus very suited to serve as a basis for the construction of different types of common boundary conditions applied to the area of mechanics. Based on the establishment of these edge conditions, several types of mixed edge boundary conditions encountered in practice can be easily established in a very simple and efficient manner without remodifying the whole algorithm. Two cases with different types of boundary conditions were investigated for demonstrating the capacities of the MDWDF approach: a thin isotropic plate and a thick isotropic plate. We were able to demonstrate, for the thin elastic plate, that the Chladni figures we have produced were similar to those produced experimentally by Waller. Results obtained from the simulations of the thick elastic plate bodies have clearly shown the plate deformation and plane-wave propagation including energy reflection and transmission from and across the boundaries. The fourth-order system we have investigated is of sufficient complexity to enable the effective consideration of a suitable parallel implementation. This is the subject of ongoing research. REFERENCES [1] A. Fettweis, “Digital filters related to classical filter networks,” AEÜ, vol. 25, pp. 79–89, 1971. , “Wave digital filters: Theory and practice,” Proc. IEEE, vol. 74, [2] pp. 270–327, 1986. [3] A. Fettweis and G. Nitsche, “Transformation approach to numerical integrating PDEs by means of WDF principles,” Multidim. Syst. Signal Processing, vol. 2, pp. 127–159, 1991. [4] , “Numerical integration of partial differential equations using principles of multidimensional wave digital filters,” J. VLSI Signal Processing, vol. 3, pp. 7–24, 1991. [5] G. Nitsche, “Numerische Lösung Partieller Differentialgleichungen mit Hilfe von Wellendigitalfiltern,” Ph.D. dissertation, Fakultät Electrotech. Info., Ruhr-Universität, Bochum, Germany, 1993. [6] X. Liu and A. Fettweis, “Multidimensional digital filtering by using parallel algorithms based on diagonal processing,” Multidim. Syst. Signal Processing, vol. 1, pp. 51–66, Mar 1990. [7] A. Fettweis and K. Meerkötter, “On parasitic oscillations in digital filters under looped conditions,” IEEE Trans. Circuits Syst., vol. CAS-24, pp. 475–481, Sept. 1977. [8] C. T. Mullis and R. A. Roberts, “Synthesis of minimum roundoff noise fixed point digital filters,” IEEE Trans. Circuits Syst., vol. CAS-23, pp. 551–562, Sept. 1976. [9] A. Fettweis, “Discrete passive modeling of physical systems described by PDEs,” in Proc. EUSIPCO-92, 6th Eur. Signal Processing Conf., vol. 1, Brussels, Belgium, 1992, pp. 55–62. [10] G. Hemetsberger, “Stability verification of multidimensional Kirchoff circuits by suitable energy functions,” in Proc. IEEE Int. Conf. Acoustics, Speech, Signal Processing, vol. 6, Adelaide, Australia, 1994, pp. 13–16. [11] C. Q. Xu, S. C. Bass, and X. Wang, “Accommodating linear and nonlinear boundary conditions in wave digital simulations of PDE systems,” J. Circuits Syst. Comput., vol. 7, pp. 563–597, 1997.
1663
[12] C. A. J. Fletcher, Computational Techniques for Fluid Dynamics.-1: Fundamental and General Techniques, 2nd ed. New York: SpringerVerlag, 1991. [13] A. Fettweis, H. Levin, and A. Sedlmeyer, “Wave digital lattice filters,” Int. J. Circuit Theory Applicat., vol. 2, pp. 203–211, 1974. [14] H. Krauss, R. Rabenstein, and M. Gerken, “Simulation of wave propagation by multidimensional digital filters,” J. Simulation Practice Theory, vol. 4, pp. 361–382, 1996. [15] C. H. Tseng and S. S. Lawson, “Incorporating initial and boundary conditions into the multidimensional wave digital filter,” in Proc. IEEE Int. Conf. Acoustics, Speech, Signal Processing, Orlando, FL, May 2002. [16] S. D. Bilbao, “Wave and scattering methods for the numerical integration of partial differential equations,” Ph.D dissertation, Dept. Elect. Eng., Stanford Univ., Stanford, CA, 2001. [17] Computer Solutions Europe AB, Partial Differential Equation Toolbox, The Math Works, Natick, MA, 1996. [18] K. F. Graff, Wave Motion in Elastic Solids. New York: Dover, 1975. [19] R. D. Mindlin, “Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates,” Trans. ASME J. Appl. Mech., vol. 18, pp. 31–38, 1951. [20] K. Liew, C. Wang, Y. Xiang, and S. Kitipornchai, Vibration of Mindlin Plates: Programming the p-Version Ritz Method, first ed. Amsterdam, The Netherlands: Elsevier, 1998. [21] E. Hinton and B. S. Al-Janabi, Numerical Methods and Software for Dynamic Analysis of Plates and Shells, E. Hinton, Ed. Swansea, U.K.: Pineridge Press, 1987. [22] H. Huang, Static and Dynamic Analyzes of Plates and Shells. London, U.K.: Springer-Verlag, 1989. [23] S. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells, 2nd ed. New York: McGraw-Hill, 1959. [24] M. D. Waller, Chladni Figures: A Study in Symmetry. London, U.K.: G. Bell and Sons, 1961. [25] R. Nouta, “The use of Jaumann structures in wave digital filters,” Int. J. Circuit Theory Applicat., vol. 2, pp. 163–174, 1974. [26] A. A. Samarskij, Theorie der Differenzenverfahren (in German). Leipzig: Akademicshe Verlagsgesellschaft, 1984.
Chien-Hsun Tseng was born in Taiwan in 1968. He received the B.Sc. degree in applied mathematics from National Cheng Kung University, Tainan, Taiwan, the M.Math. degree from University of New South Wales, Sydney, Australia, and the Ph.D. degree in electrical engineering from Curtin University of Technology, Perth, Australia, in 1995, 1997, and 2000, respectively. From 2000 to 2001, he was a Research Scientist at Weierstrass Institute for Applied Analysis and Stochastics (WIAS), Berlin, Germany. Since June 2001, he has been with the University of Warwick, Warwick, U.K., where he is a Research Fellow in the School of Engineering. His current research interests include optimal digital and statistical filters design through optimization techniques, modeling of physical systems using digital signal processing techniques, and stochastic and numerical optimizations.
Stuart Lawson (M’80–SM’91) received the B.Sc. degree in mathematics (with first class honors) from Woolwich Polytechnic and the D.I.C., U.K., and the Ph.D. degree in electrical engineering from Imperial College, University of London, London, U.K., in 1971 and 1975, respectively. From 1971 to 1978, he was a Research Engineer with Plessey Telecommunications, U.K. He taught and researched at City University, London, U.K., from 1978 to 1987. Since 1988, he has been at the University of Warwick where he is now a Reader in the School of Engineering. His current research interests include digital filters and filter banks, image compression, modeling physical systems using digital signal processing techniques and sonar signal processing. Dr. Lawson is a Fellow of Institution Electrical Engineers and a member of the Advisory Panel for the IEE Signal Processing Professional Network.
