Estimates for Divergence Velocities of Axially Moving Orthotropic Thin Plates

Some models for axially moving orthotropic thin plates are investigated analytically via methods of complex analysis to derive estimates for critical plate velocities. The linearized Kirchhoff plate theory is used, and the energy forms of steady-state models are considered with homogeneous and inhomogeneous tension profiles in the cross direction of the plate. With the help of the energy forms, some limits for the divergence velocity of the plate are found analytically. In numerical examples, the derived lower limits for the divergence velocity are analyzed for plates with small flexural rigidity.


INTRODUCTION
In industrial processes involving axially moving materials, such as the production of paper, steel, or textiles, high transport speed is desired, but it also may cause loss of stability. In the modeling of such systems, the dynamic behavior of strings, membranes, beams, and plates has been of general interest taking into account the transverse, Coriolis, and centripetal accelerations of the material motion. The first studies on them include Skutch (1897), Sack (1954), Archibald and Emslie (1958), Miranker (1960), Swope and Ames (1963) and Mote (1968Mote ( , 1972Mote ( , 1975. Sack (1954), Archibald and Emslie (1958), and Simpson (1973) studied the effects of axial motion on the frequency spectrum and eigenfunctions. In their research, it was shown that the natural frequency of each mode decreases as the transport speed is increased, and that the traveling string and beam both experience divergence instability at a sufficiently high speed. Wickert and Mote studied stability of axially moving strings and beams using modal analysis and a Green's function method (Wickert and Mote, 1990). They presented the expressions for the critical transport velocities analytically. However, recently, Wang et al. (2005) showed analytically that no static instability occurs for the transverse motion of an ideal string at the critical velocity. For axially moving beams with small flexural stiffness, Kong and Parker (2004) found closed-form expressions for the approximate frequency spectrum by a perturbation analysis.
The stability of travelling rectangular membranes and plates was first studied by Ulsoy and Mote (1982), and Lin andMote (1995, 1996). The stability of out-ofplane vibrations of axially moving rectangular membranes was studied by Shin et al. (2005). For the behavior of the membrane, it was found that the motion remains stable until a critical speed, at which static instability occurs. Lin (1997) studied stability of axially moving plates, and numerically showed that loss of stability of the plate occurs in a static (divergence) form at a sufficiently high speed. Banichuk et al. (2010) considered stability and studied the critical velocity and the corresponding critical shapes of an axially moving elastic plate. Later on, Banichuk et al. (2013b) applied the stability analysis to optimization of the magnitude of applied tension, taking also fatigue fracture into account.
In the recent studies concerning axially moving plates, material properties such as orthotropicity (Banichuk et al., 2011) or viscoelasticity (Marynowski, 2010) have been taken into consideration, and their effects on the plate behavior have been investigated. In Banichuk et al. (2011), divergence instability for travelling orthotropic rectangular plates, with two opposite edges simply supported and the other two edges free, was studied and an explicit expression for the limit velocity of stable axial motion was found. Hatami et al. (2009) studied free vibration of the moving orthotropic rectangular plate in sub-and supercritical speeds, and flutter and divergence instabilities at supercritical speeds. Their study was limited to simply supported boundary conditions at all edges.
Free vibrations of classical rectangular plates, which are not moving axially, have been discussed in the book by Gorman (1982). The case of orthotropic plates, specifically, has been studied by Biancolini et al. (2005) including all combinations of simply supported and clamped boundary conditions on the edges. Xing and Liu (2009) obtained exact solutions for free vibrations of stationary rectangular orthotropic plates considering three combinations of simply supported (S) and clamped (C) boundary conditions: SSCC, SCCC, and CCCC. Kshirsagar and Bhaskar (2008) studied vibrations and buckling of loaded stationary orthotropic plates. They found critical loads of buckling for all combinations of boundary conditions S, C, and F.
Tension inhomogeneities and their effects on the divergence instability of moving plates have been studied in Banichuk et al. (2013a). A linearly inhomogeneous tension profile, with the axial tension varying across the cross direction, was considered in the case of a moving isotropic plate. The inhomogeneities in tension were found to have a significant effect on the critical velocity of the plate and to change the buckling shapes dramatically compared to the case of homogeneous tension.
In this study, we investigate the energy forms corresponding to (the steady state of) a travelling orthotropic plate under homogeneous or inhomogeneous tension. It is shown that the critical velocity for an orthotropic plate subjected to any tension field (satisfying certain conditions) is greater than the critical velocity of an ideal membrane subjected to the same tension field. The differential equations for a travelling orthotropic plate subjected to an arbitrary tension field are derived from the corresponding energy form. For a linearly inhomogeneous tension profile, the stress field is solved with the help of the (Airy) stress function. For this type of inhomogeneity, we derive a guaranteed lower limit for the divergence velocity of the orthotropic thin plate and present a numerical example.

AXIALLY MOVING ORTHOTROPIC PLATES
Consider an axially moving orthotropic rectangular plate in a Cartesian coordinate system, supported at two edges x = 0 and x = . The plate is assumed to be subjected to tension T at the supported edges. The problem setup is shown in Fig. 1. The plate width is 2b, its thickness is h, and the length of the span is . Throughout this study, the plate is assumed to travel axially at a constant velocity V 0 . We denote the transverse displacement of the plate by the function w x y and the considered rectangular part of the plate by : The material parameters for the orthotropic plate are denoted by m (the mass per unit area), 12 and 21 (the Poisson ratios in the xy plane), E 1 and E 2 (the Young's moduli in the x and y directions, respectively), and G 12 (the shear modulus in the xy plane). We denote the components of in-plane tension by T xx , T yy and T xy = T yx .
First, we introduce bilinear forms (for complex-valued functions) that correspond to the energy of the described system. Bilinear energy forms, e.g., for stationary isotropic plates, have been discussed by Chen et al. (1998).
We consider bilinear forms corresponding to the strain energy due to in-plane tension m , the strain energy due to bending b , and the kinetic energy of the moving plate.
The bilinear form corresponding to the strain energy due to in-plane tension is given by (Timoshenko and Woinowsky-Krieger, 1959) wherev denotes the complex conjugate of v. The complex conjugate is needed, because, e.g., in linear stability analysis, in general the displacement function is allowed to be complex-valued. Using the linearity of the governing partial differential equation, the solution then splits

297
into two real-valued functions, being the real and imaginary parts of the complex displacement, respectively.
The bilinear form corresponding to the strain energy of bending is (Timoshenko and Woinowsky-Krieger, 1959) where the coefficients D i and A i for the flexural rigidities are defined as and We have assumed elastic compatibility, i.e., The bilinear form corresponding to the kinetic energy is given by is the operator of the total (or material, or Lagrange) derivative.
In this study, we concentrate on steady-state problems only. Note that when the axially moving material is viewed in the laboratory frame, the kinetic energy (7) has also a time-independent component, which is due to a centripetal contribution. For the rest of the discussion, we will denote by the time-independent part of the kinetic energy: Consider the form corresponding to the total energy of the considered axially moving orthotropic thin plate.
In the following, we derive the partial differential equation corresponding to w v = 0 for the energy defined in (10).

SAKSA AND JERONEN
We perform integration by parts for all the energy forms m , b and . For m in (2), integration by parts gives The in-plane tension components are assumed to satisfy the following equilibria: which describe the in-plane balance of linear momentum in the absence of external body forces and in-plane accelerations. On the free boundaries, the following freeof-traction boundary conditions are set: These boundary conditions require that the normal component of the stress tensor vanishes on the free boundaries. Boundary conditions for tension components on the supported edges depend on the case studied, and will be specified further below. Due to linearity,v satisfies the boundary conditions if and only if v does. We assume that for v, we have v =v = 0 at x = 0 . Thus, we obtain Integration (twice) by parts of b in (3) gives For w, we assume the (simply supported and free) boundary conditions Since w = 0 andv = 0 at x = 0 , also the derivatives of w andv in the y direction vanish at those edges, i.e., 2 w/ y 2 = 0 and v/ y = 0 at x = 0 . Thus, all the boundary terms in (15) zero out and what remains is Finally, for the energy form , the integration by parts gives Combining the results in (14), (20), and (21), the partial differential equation corresponding to subject to the boundary conditions (17)-(19) and tension equilibria (12). Equation (23) describes small transverse displacements of an axially travelling elastic plate when the material flow occurs in a steady state; the function w describes the displacement as it is appears in the laboratory coordinate system. In-plane tensions and stresses xx , xy , and yy are related by

SAKSA AND JERONEN
where h is the thickness of the plate (assumed constant). Referring to (12), the inplane stresses xx , xy and yy satisfy the equilibrium equations

AIRY STRESS EQUILIBRIUM FOR ORTHOTROPIC PLATES
In the following, we will present the Airy stress equilibrium for an orthotropic plate. With the assumption of small deflections, the strains xx , yy , and xy are defined with the help of the in-plane displacements u and v (in the x and y directions, respectively) as follows: for which the following equation holds: Stresses and strains are related to each other by Hooke's law (here, inverse relation): We introduce the stress function F : The stresses in (28) satisfy automatically (25). Inserting (28) and (27) into (26), we obtain (with the help of the elastic compatibility relation (6)): This is the Airy stress equilibrium for orthotropic plates, given also, e.g., in the book by Marynowski (2008). For isotropic plates, the equilibrium (29) is the biharmonic equation for the Airy stress function F . Note that (6) specifies E 1 21 = E 2 12 ; the coefficient of the mixed term only depends on quantities that relate symmetrically to both material axes. The orthotropic model can be reduced to an isotropic model by choosing G 12 = G H (Huber, 1923;Timoshenko and Woinowsky-Krieger, 1959), where In practice, the measured values for the shear modulus G 12 may significantly differ from this ideal value (Seo, 1999;Yokoyama and Nakai, 2007). In such cases, the full orthotropic model must be used.

ESTIMATES FOR DIVERGENCE VELOCITIES
Previously, we introduced the bilinear forms m , b , and corresponding to the strain energy due to in-plane tension, the strain energy due to bending, and the kinetic energy, respectively. Now, we derive some estimates based on them, and finally by combining these estimates, we obtain lower limits for the critical velocities of axially moving orthotropic plates and (ideal) membranes.
We begin with the bilinear form m defined above in (2). Let us assume that the in-plane tension components T xx , T yy , and T xy satisfy the following properties: T xy ≥ 0 in (33) where 2 is a positive constant that can be chosen freely. The reason for the requirements (31)-(33) will be justified below.
We will denote the norm of a complex-valued function (or scalar) a by a . Recall that, for complex-valued a, it holds that a 2 = aā. We have for m w w With the assumptions in (31)-(33), all the quantities on the last row of (34) are nonnegative. Thus, we obtain Similarly, we can estimate b w w corresponding to the strain energy due to bending. We have The constant D 2 − A 2 1 /D 1 ≥ 0, since and 0 ≤ √ 12 21 ≤ 1/2 and D 2 > 0. Thus, we see that For a nontrivial solution w (real part or imaginary part ≡ 0) with boundary conditions (17), it holds that We may now construct some estimates for the divergence velocities of axially moving plates and membranes. The divergence velocity V div 0 for an axially moving orthotropic plate is solved from the equation and is expressed as In the case of a constant tension T 0 at the supported edges, we have T xx = T 0 , T yy = 0, and T yy = 0. The divergence velocity V div 0 T xx =T 0 can be expressed in this special case as In the case of an ideal membrane, it holds that b w w = 0, and the divergence velocity V div,mem 0 (for a general tension profile) can be expressed as As a first result for divergence velocities, we see that the divergence velocity of an orthotropic plate is always greater than the divergence velocity of an ideal membrane. Since, by (38), we have b w w ≥ 0, from (41) and (43) we get the estimate Provided that (31)-(33) are satisfied by the in-plane tension components, we have by (35) and (43) that V div,mem 0 2 ≥ 0, and thus also Consider now a case in which the tension profile is linear at the supported edges x = 0 and x = . That is, T xx x=0 = T 0 + y where T 0 and are positive constants such that T xx is nonnegative along the edge. The parameter will be called the tension profile skew parameter. We introduce boundary conditions for the stresses: With the help of the boundary conditions in (46) and the relations in (28), we find the boundary conditions for F . Note that above the last two boundary conditions in (46) guarantee that (13) is satisfied. The solution to the boundary value problem (29) (with boundary conditions obtained from (46)) is where c 0 , c 1 and c 2 are arbitrary real constants. Now the tensions (inside ) are Assuming that tension in the x direction is positive in , i.e., assuming that we may derive more interesting estimates for the critical velocities. Inserting (48) into (34), we obtain where T 0 − b is a positive constant.

SAKSA AND JERONEN
Now, we can estimate the divergence velocity of an axially moving plate under linear tension by the divergence velocity of an axially moving plate under constant tension. By (50), (41), and (42), we have (for the absolute values of the compared velocities) The last inequality follows from (45).
In (51), we have two interesting results for the case of a linear tension profile. First, the critical velocity of an orthotropic plate under a linear tension profile can be estimated from below by the critical velocity of a plate under constant tension having the value that is the minimum of the considered linear tension profile. This result is practical, since the critical velocity of the plate under constant tension can be found analytically (e.g. Banichuk et al., 2011). Second, the divergence velocity of a plate under a linear tension profile is also estimated from below by the divergence velocity of an ideal membrane with constant tension having a simple analytical formula.
In the case of constant tension, we also have for the critical velocity that V div 0 ≥ T 0 /m. The result for constant tension and the first estimate in (51) in the case of an ideal membrane were also discussed in Banichuk et al. (2013a), here however in a slightly different manner using complex analysis.

STATIC STABILITY ANALYSIS FOR A PLATE UNDER LINEAR TENSION DISTRIBUTION
The partial differential equation for an axially moving orthotropic plate under a linear tension profile is where We present the solution of (52) and (17)-(19) in the following form: where f y/b is an unknown function. Introducing a new variable = y/b and inserting (54) into (52), we obtain the eigenvalue is defined as and the dimensionless bending rigidities are In (58), the normalization constant D 0 can be chosen freely, e.g., D 0 = D 1 . The boundary conditions (18)-(19) become The parameters 1 and 2 are given above in (16). Note that for an isotropic material H 1 = H 2 = H 3 = 1 with D 0 = D, and 1 = and 2 = 2 − . For comparison, see Banichuk et al. (2013a).
The strong form (55), (59)-(60) was discretized directly, with classical central differences of second-order asymptotic accuracy. Explicitly, the finite difference formulae used were where is the grid spacing, f is any sufficiently differentiable function, and the subscripts denote nodes in a uniformly spaced grid, f j ≡ f j and j = −1 + j − 1 for j = 1 2 N , with the last point placed at N = −1 + N − 1 = 1

SAKSA AND JERONEN
(which, given the number of nodes N , determines the grid spacing ). The relations (62) and (64) were used to discretize Eq. (55).
To account for the boundary conditions (59)-(60), the method of virtual points was used. The virtual points −1 and 0 were introduced at the left end of the domain, and the virtual points N +1 and N +2 at the right end of the domain. Columns were added to the discrete equation system to account for the corresponding virtual degrees of freedom f −1 , f 0 , f N +1 , and f N +2 . The boundary conditions (59)-(60)-which hold at the points 1 and N -were discretized using the formulae (61)-(63), allowing the virtual degrees of freedom to appear in the discretization. The resulting rows, representing the discretized boundary conditions, were then added to the discrete equation system.
Because the boundary conditions are homogeneous, it was possible to add them to the discrete system by rewriting the original discrete problem as a generalized linear eigenvalue problem where B is an identity matrix with the first two and last two rows zeroed out, A contains the differential operators (in (55) on the left-hand side), f is the discretized form of f . In (65), the first two and the last two rows of A contain the discretized boundary conditions. As for the numerical solution of the generalized eigenvalue problem (65), it is important to consider the symmetry properties of the square matrices A and B. The classical central finite difference formulae for the even (second and fourth) order derivatives, (62) and (64), are symmetric. Thus, the submatrix of A without the rows and columns corresponding to the boundary conditions is symmetric, since on the left-hand side of (55), there are only even-order derivatives of f . However, the boundary condition (60) contains odd-order derivatives, for which the central finite difference formulae, (61) and (63), are not symmetric. Altogether, the matrix A is not symmetric, and in the numerical eigenvalue solver, we apply the generalized Schur decomposition, also known as the QZ decomposition, applicable for nonsymmetric matrices.
In the following, we have studied the critical velocity of the axially moving orthotropic plate and the corresponding critical velocity in the case of thin material having a very small flexural rigidity. All the material and physical parameters have been chosen such that the model represents a paper material.
In Table 1, the results for the critical velocities are shown for three different values of the shear modulus G 12 . It can be seen that the smaller the value of the shear modulus is, the lower the critical velocity V cr 0 becomes. This effect was found also in the case of the isotropic plate (Banichuk et al., 2013a). The critical velocity for an orthotropic plate with G 12 = G H and homogeneous tension ( = 0) is V cr 0 = 83 4461 m/s (by keeping other parameter values the same). We see that for a large value of the shear modulus (1 3G H ), the critical velocity for a plate with an inhomogeneous tension profile can be larger than that of a plate with homogeneous tension but with a smaller value of the shear modulus.
In Fig. 2, the effect of the shear modulus variation on the buckling mode is visualized by presenting the displacement w at x = /2. As expected, the larger the value of the shear modulus is, the smaller the change in the buckling mode becomes when compared to the case with a homogeneous tension profile. In Fig. 3, both the effect of the value of the shear modulus G 12 and the effect of tension inhomogeneity are illustrated. The changes in the shapes of the critical eigenmodes correspond to the changes observed for the critical velocities (Table 1): the more symmetric (with respect to y) the critical eigenmode, the larger the critical velocity.
We also computed the relative differences (estimation errors) between the velocities compared in (51). We compared the values between the divergence velocities of the plate under linear tension V div 0 T xx =T 0 + y and the plate under constant tension V div 0 T xx =T 0 − b . For comparison with different estimates for V div 0 T xx =T 0 + y , we computed also the difference between it and V div, mem 0 T xx =T 0 − b . The relative differences were computed, i.e., how much (in terms of percentage) the estimate from below differs from the actual computed velocity.
In the computations for the relative differences, we let the values of E 1 and E 2 vary, since they are responsible for the magnitude of the flexural rigidity. Their relative magnitudes were kept the same as above, i.e., E 1 = 2E 2 . The in-plane shear modulus was chosen as G 12 = G H . For the study in which E 1 was varied, the relative tension profile skew parameter was kept constant at / max = 10 −5 . The relative differences between the divengence velocity and its estimates from below were also studied with respect to the relative tension profile skew parameter / max . In these computations, E 1 was chosen to be relatively small E 1 = 2 GPa (and thus E 2 = 1 GPa). Figure 4 compares the values of the critical velocity for a plate with linear tension distribution T 0 + y with the critical velocity for an ideal membrane with constant tension T 0 − b. Although T 0 − b /m gives a guaranteed estimate from below for the critical velocity V div 0 T xx =T 0 + y , we see that it does not provide a good approximation of the actual critical plate velocity.
In Fig. 5, the relative difference in the plate divergence velocity in the cases of linear and constant tension profiles is shown. From the left-hand side figure, we see that even for quite large values of / max (up to 10 %), the constant tension case gives a good approximation for the critical velocity of the plate with a linear tension distribution.

Figure 4
Relative difference between the critical velocities for a plate with a linear axial tension profile T xx y = T 0 + y and a membrane under constant tension T xx = T 0 − b applied at the edges x = 0 and x = . Left: The relative difference with respect to the value of / max . The values E 1 = 2 GPa, E 2 = 1 GPa were kept constant. Right: The relative difference with respect to E 1 . The tension profile skew parameter was kept constant at / max = 10 −5 .

Figure 5
Relative difference between the critical velocities for plates with a linear axial tension profile T xx y = T 0 + y and a constant tension T xx = T 0 − b applied at the edges x = 0 and x = . Left: The relative difference with respect to the value of / max . The values E 1 = 2 GPa, E 2 = 1 GPa were kept constant. Right: The relative difference with respect to E 1 . The tension profile skew parameter was kept constant at / max = 10 −5 .

CONCLUSION
Models for axially moving orthotropic plates under an inhomogeneous tension profile were studied. For the critical velocity of a moving orthotropic plate under any tension profile satisfying certain conditions, it was shown that the critical velocity is higher than that of an ideal membrane. The case of a plate under an inhomogeneous axial tension profile, with tension varying linearly in the cross direction, was studied in detail. It was shown analytically that the values for the divergence velocity of the plate can be estimated from below by analytically computable limits. A numerical example about the effects of the shear modulus on the critical velocities and buckling modes was given. It was seen that the greater the value of the shear modulus, the smaller the change in the buckling mode, when compared to the case with a homogeneous tension profile. The effect of the shear modulus on the critical velocity was minor, but for greater values of the shear modulus the plate was found to be more stable. The stabilizing effect of the shear modulus became more pronounced as the skew in the linear tension profile was increased.
The analytical results, especially the limits for the critical velocities, can be helpful for approximating the critical velocities in the cases of more complicated tension fields. Importantly, the limits for the critical velocities also guarantee that the results given by a numerical algorithm will be physically meaningful. For engineering applications, analysis for more general tension fields than those with a linear profile at the supported edges would be desirable. In those cases, the solution of the tension field itself becomes a more challenging problem (Gorman and Singhal, 1993). Developing good velocity estimates for more general cases remains a target for future studies.

FUNDING
This research was supported by the Academy of Finland (grant no. 140221), the Jenny and Antti Wihuri Foundation, and the Finnish Cultural Foundation.