Added-Mass Based Efficient Fluid–Structure Interaction Model for Dynamics of Axially Moving Panels with Thermal Expansion

Abstract

The paper considers the analysis of a traveling panel, submerged in axially flowing fluid. In order to accurately model the dynamics and stability of a lightweight moving material, the interaction between the material and the surrounding air must be taken into account. The lightweight material leads to the inertial contribution of the surrounding air to the acceleration of the panel becoming significant. This formulation is novel and the case complements our previous studies on the field. The approach described in this paper allows for an efficient semi-analytical solution, where the reaction pressure of the fluid flow is analytically represented by an added-mass model in terms of the panel displacement. Then, the panel displacement, accounting also for the fluid–structure interaction, is analyzed with the help of the weak form of the governing partial differential equation, using a Galerkin method. In the first part of this paper, we represent the traveling panel by a single partial differential equation in weak form, using an added-mass approximation of the exact fluid reaction. In the second part, we apply a Galerkin method for dynamic stability analysis of the panel, and present an analytical investigation of static stability loss (divergence, buckling) based on the added-mass model.

1. Introduction

In this study, we develop a mathematical model representing the behavior of an open draw in the paper making process, taking into account the effects of fluid–structure interaction and of thermal expansion of the material. Open draws are the spans, where the web travels without support and interacts with gas (fluid). Within these open draws, the paper web is subject to temperature forces. Stability of these parts is therefore an interesting research question for the applications. As for applied aspects of this study, axially moving materials have many applications in industry. Examples of systems with axially moving materials are band saws, power transmission belts, paper making processes, printing presses, manufacturing of plastic films and sheets. The model can also be applied to other open draws involving lightweight materials, such as the vibrations of a tape in a tape drive. In our paper, we however concentrate on paper making processes. Examples including transportation from the press nip to the drying section.

We have simplified the problem, while capturing relevant aspects of the physical phenomena. A key point is that the behavior of the moving material is strongly dependent on problem parameters such as panel and fluid velocities, densities of the fluid and the panel, bending rigidity, temperature, in-plane tension, and also external forces and geometric factors. For lightweight materials, the interaction between the material and the surrounding air must be taken into account because it significantly changes the dynamic response, as was first observed by [1].

In many previous investigations, axially moving materials have been modeled as traveling flexible strings, beams, panels and plates. Some classical articles have been made by Archibald and Emslie [2], Mote [3], Wickert and Mote [4], Wickert [5], Parker [6] and Wang et al. [7]. The majority of the performed studies have been devoted to isotropic moving materials. As for orthotropic moving materials, many important results, concerning modeling and the analysis of behavior, can be found in the monographs by Marynowski [8], Banichuk et al. [9], Tuovinen [10], and again Banichuk et al. [11].

In the article by Liu et al. [12], free vibrations of nonuniform rectangular plates varying in one or two directions are investigated. The plate was orthotropic and resting upon a Winkler-type spring foundation. In the article by Ghayesh and Amabili [13], the investigation of the transverse vibration response, stability, and bifurcations of an axially moving viscoelastic beam with time-dependent axial speed was considered. Research by Hozhabrossadati and Sani [14] deals with the determination of the static displacement function of a Euler–Bernoulli beam with two guided supports. An extensive literature review about mechanics of axially moving continua can also be found in Marynowski & Kapitaniak [15]. Note also that the papers by Banichuk et al. [16] and Banichuk and Ivanova are thematically adjointed [17].

Fundamental problems of fluid–structure interaction have been studied (Banichuk and Mironov [18,19], Ashley and Landahl [20], Anderson [21] and analyzed Bisplinghoff and Ashley [22], Ashley and McIntosh [23], Kornecki et al. [24], Vaughan and Raman [25], and Banichuk et al. [9]).

In a book by Banichuk et al. [9], following an approach similar to Ashley and Landahl [20], in the case of potential flow, the aerodynamic reaction was found as a functional fluid–structure interaction problem involving the solution of integrodifferential equation with a singular kernel. The problem can be further simplified by approximating the singular kernel. For a review on studies concerning circular cylindrical shells including fluid–structure interaction, see Amabili and Paidoussis [26]. For shells containing flowing fluid, see Amabili et al. [27,28], Amabili and Garziera [29,30]. Periodically supported plates subjected to axial flow have been considered in Tubaldi et al. [31].

In Liu and Chang [32], a fluid–structure interaction problem dealing with free vibration analysis of an axisymmetrical, nonuniform circular plate with quadratic change of thickness in contact with fluid was investigated. Moreover, taking into account thermoelastic behaviour, the forced vibration analysis of nonhomogeneous thermoelastic, isotropic, thin annular disk under periodic and exponential types of axisymmetric dynamic pressures applied on its inner boundary has been performed and an analytical benchmark solution has been obtained, e.g., Mishra et al. [33]. In the article by Mao et al. [34], the authors have studied the thermoelastic instability (TEI) of a functionally graded material (FGM) half-plane sliding against a homogeneous half-plane at the in-plane direction. The interaction of the frictional heat and thermal contact resistance is taken into account in the TEI analysis.

In this study, we formulate the dynamic aerothermomechanical problem for an axially moving panel, interacting with potential flow and subject to thermal expansion. This formulation is novel and the case complements our previous studies on the field. As for simplification of the model, we approximate the rigorous solution for the fluid reaction via an added-mass approach. As a result, the integrodifferential governing equation reduces to a partial differential equation, where we use added-mass approximations of the original inertial, Coriolis and centrifugal terms. In a steady state, the equation is further reduced into an ordinary differential equation. The resulting differential equations can be effectively solved for various initial and boundary conditions. As a result, we find solutions for some stationary and time-dependent problems. In the stationary case, we find the analytical solution of the stability problem and present expressions for the critical divergence velocity and for critical temperature of the loss of stability. We also study the dynamic stability of the system.

The proposed formulation of problem and method of analysis are simplified to give clear physical motivation and realize effective solutions. First, simplification consists of representation of moving elastic plate (two-dimensional model) as moving in the axial direction elastic panel, describing by a one-dimensional deflecting function $w=w(x,t)$. The travelling panel is interacting with a two-dimensional model of fluid moving in axial and normal (transverse) directions. These assumptions give us the possibilities to apply complex variable techniques for finding analytical solutions of a two-dimensional hydrodynamical (plane) problem with arbitrary distribution of elastic deflections $w(x,t)$. Thus, the original coupled problem of fluid–structure interaction admits the reduction to a one-dimensional integro-differential equation. This equation can be used as a test for some numerical solutions with 3D hydrodynamical and 2D elastodynamical parts of the original problem.

The analytical approach, proposed in this paper, has some constraints. The model does not cover all possible aspects of the physical situation. The one-dimensional geometry of the panel limits the applicability to long and narrow draws. For short and wide draws, a plate model should be used, due to the localization of deflections near the free edges (Banichuk et al. [9]). In addition, paper as a material is slightly viscoelastic. It is known that the presence of viscosity, no matter how small, will qualitatively change the behavior of the stability exponents. However, the changes appear in the post-critical range, and the first critical velocity remains almost the same, so, for predicting the initial loss of stability, the elastic material model is sufficient. Moreover, the moving fluid is supposed to be ideal and the property of viscosity is excluded. If in reality liquid or gas are not ideal, but the viscosity is small, then the obtained ideal solution can be considered as a zero-order approximation in the small parameter method. Perturbation with respect to small viscosity parameter $\mu$ will be realized by the first order approximation.

Finally, we have not accounted for the temperature dependence of the Young’s modulus. This is expected not to have a significant effect on the results. The bending rigidity parameter is small, so it mainly has a singular perturbation effect on the partial differential equation analyzed. For small but finite bending rigidity, the behavior of the stability exponents describing free vibrations is qualitatively different from the case with no bending rigidity (Wang et al. [7] and Jeronen [35]), but the first critical velocity remains almost the same. Thus, even if the bending rigidity slightly varies due to the temperature dependence of the Young’s modulus, it will not have a major effect on the critical velocity predicted by the model.

2. Basic Relations and Aerothermoelastic Model

Consider a traveling thermoelastic panel subjected to a potential flow, where the free stream flows toward the right at velocity ${\nu }_{\infty }$ (with respect to the laboratory frame), see Figure 1. The governing equation describing small transverse vibrations of the panel is

$$m{\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}+2m{V}_0{\displaystyle \frac{{\partial }^{2}w}{\partial x\partial t}}+m{V}_0^2{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}-T{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+D{\displaystyle \frac{{\partial }^{4}w}{\partial x^4}}=q_{\mathrm{f}}+g,$$

(1)

where w is the transverse displacement, $q_{\mathrm{f}}\equiv q_{\mathrm{f}}\left(w\right)$ is the fluid reaction pressure, and $g\equiv g(x,t)$ represents external forces inside the domain. The mass per unit area is m. The panel travels axially at the velocity $V_0$ and is subjected to in-plane axial mechanical tension $T_{\mathrm{m}}$ and thermal compression $T_{\theta }$. The unit of $T_{\mathrm{m}}$ and $T_{\theta }$ is force per unit length. The bending rigidity of the panel is D. We have

$$T=T_{\mathrm{m}}-T_{\theta },T_{\mathrm{m}}=T_{0,}{T}_{\theta }={\displaystyle \frac{Eh}{1-\nu }}{\epsilon }_{\theta },D={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}}.$$

(2)

Here, h is the thickness of the panel, E is the Young’s modulus of the panel material, $\nu$ is its Poisson ratio, $T_0$ is a given constant tension, and ${\epsilon }_{\theta }$ is the generalized thermal strain defined by

$${\epsilon }_{\theta }={\displaystyle \frac{1}{h}}{\int }_{-h/2}^{h/2}{\alpha }_{\theta }\theta \left(z\right)\mathrm{d}z.$$

(3)

The coefficient of linear thermal expansion ${\alpha }_{\theta }$ is constant, and the function $\theta$ takes the form $\theta \left(z\right)={\theta }_{\mathrm{abs}}\left(z\right)-{\theta }_0$, where ${\theta }_{\mathrm{abs}}$ and the reference temperature ${\theta }_0$ are given in Kelvin. If the panel is at a uniform temperature of ${\theta }_0$, then $\theta \left(z\right)=0$, and consequently ${\epsilon }_{\theta }=0$. In the case where ${\alpha }_{\theta }$ and $\theta$ are constant, we have ${\epsilon }_{\theta }={\alpha }_0\theta$.

The reaction pressure $q_{\mathrm{f}}$ is unknown, to be solved from the flow model in terms of w. The reaction pressure describes how the surrounding fluid pushes back on the panel, in reaction to the panel vibration. The external forces g are considered known, and are allowed to vary dynamically. For finding the fluid reaction pressure, we apply techniques from the aerodynamics of thin aerofoils, constructing a Green’s function type solution, as was made in [18,19] and presented in the books [9,20,21,22,36]. In non-dimensional coordinates $x^{\prime }=x/\ell$, $t^{\prime }=t/\tau$ we have (primes omitted)

$$q_{\mathrm{f}}(x,t)=-{\displaystyle \frac{{\rho }_{\mathrm{f}}}{\ell }}\left({\displaystyle \frac{\ell }{\tau }}{\displaystyle \frac{\partial }{\partial t}}+v_{\infty }{\displaystyle \frac{\partial }{\partial x}}\right){\int }_{-1}^{1}N(\xi ,x)\left({\displaystyle \frac{\ell }{\tau }}{\displaystyle \frac{\partial }{\partial t}}+v_{\infty }{\displaystyle \frac{\partial }{\partial x}}\right)w(\xi ,x)\mathrm{d}\xi ,$$

(4)

$$N(\xi ,x)\equiv {\displaystyle \frac{1}{\pi }}ln\left|{\displaystyle \frac{1+\mathsf{\Lambda }(\xi ,x)}{1-\mathsf{\Lambda }(\xi ,x)}}\right|,\mathsf{\Lambda }(\xi ,x)=\sqrt{{\displaystyle \frac{(1-x)(1+\xi )}{(1-\xi )(1+x)}}},$$

(5)

$$\left(x,t\right)\in \Omega ,\Omega =\left\{\left(x,t\right):-1<x<1,0<t<\infty \right\},$$

(6)

where $v_{\infty }=0$ means no free-stream flow in laboratory coordinates and $v_{\infty }=V_0$ means that the whole air mass moves with the panel.

Briefly stated, this Green’s function solution describes two-dimensional potential flow in the plane excluding a linear cut. To make an analytical solution possible, since we consider small vibrations around the trivial equilibrium position only, we geometrically approximate the panel as a straight line segment. The region $z=0$, $-1<x<1$ is cut out from the plane, approximating the space occupied by the panel. The vibrations of the panel are taken into account in the no-penetration boundary condition on the velocity of the flow. This leads to Equations (4)–(6). A detailed treatment can be found in the book [9].

In (4), ${\rho }_{\mathrm{f}}$ is the density ($\mathrm{kg}/{\mathrm{m}}^3$) of the surrounding medium. The spatial scaling factor ℓ is the half-length of the span, and $\tau$ is an arbitrary scaling factor for time coordinate. For a physically meaningful time scaling, one can choose, e.g., $\tau =\ell /C$, where $C=\sqrt{T/m}$, the critical velocity of a traveling ideal string or membrane. The unit of $\tau$ is $\left[\tau \right]=s$.

Added-mass models are sometimes used for taking into account the inertial effects of fluid–structure interaction, expressed by an inertial operator in a simplified approximate setting. In this study, we consider an added-mass approximation for the above model of fluid–structure interaction of an axially moving panel. We will use the constructed approximation of the total inertial term $m{\mathrm{d}}^{2}w/\mathrm{d}{t}^2$ for the solution of the problems of vibration and stability. As was shown in the book [9], we may use the following approximation:

$${\int }_{-1}^{1}N(\xi ,x)f\left(\xi \right)\mathrm{d}\xi \approx \mu {\int }_{-1}^1\delta (\xi ,x)f\left(\xi \right)\mathrm{d}\xi =\mu f\left(x\right),x\in \left[-1,1\right],$$

(7)

where $\delta (\xi ,x)$ is the Dirac delta distribution, and $\mu$ is the constant

$$\mu \equiv \stackrel{\mathrm{mean}}{{}_{{}_{x\in (-1,1)}}}{\int }_{-1}^{1}N(\xi ,x)\mathrm{d}\xi ={\displaystyle \frac{1}{2}}{\int }_{-1}^1\left[{\int }_{-1}^{1}N(\xi ,x)\mathrm{d}\xi \right]\mathrm{d}x\approx {\displaystyle \frac{\pi }{4}}.$$

(8)

The symbol $\mathrm{mean}(·)$ denotes the average of the indicated quantity over the indicated interval. The idea behind the approximation (7) is that $N(\xi ,x)$ is a singular but integrable ($L^1$) kernel, which decays quickly as $\left|\xi -x\right|$ increases; it can be shown that $N(\xi ,x)<ln(4/\left|\xi -x\right|)$. Hence, for any sufficiently well-behaved f, most of the contribution to the integral on the left-hand side comes from the singularity.

In physical terms, we approximate the fluid reaction pressure as a pointwise local effect, whereas in reality a disturbance in the flow propagates along the fluid, causing the pressure field to change globally (limited by the speed of sound in the medium). Especially from a mathematical modeling perspective, recall that a potential flow will instantly reconfigure itself to satisfy the Laplace equation. Any disturbance anywhere in the fluid domain will instantly cause global (even if small) effects. These effects are ignored in the present approach; this makes the solution process significantly simpler at the cost of some accuracy. Another approach is to treat the singular integrals numerically, which has been performed earlier in [9,35].

Let us now apply (7) and (8) to Equation (1). For now, consider only the inertial terms and $q_{\mathrm{f}}$. In other words, in Equation (1), let $T=0$, $D=0$ and $g\equiv 0$, as we may easily add these terms back when finished. We approximate $q_{\mathrm{f}}$ by inserting Equation (7) and $\mu =\pi /4$ from Equation (8) into Equation (4). We have

$${\displaystyle \frac{m}{{\tau }^2}}{\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}+2{\displaystyle \frac{m{V}_0}{\ell \tau }}{\displaystyle \frac{{\partial }^{2}w}{\partial x\partial t}}+{\displaystyle \frac{m{V}_0^2}{{\ell }^2}}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}=-{\rho }_f{\displaystyle \frac{\pi }{4}}\left({\displaystyle \frac{\ell }{{\tau }^2}}{\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}+2{\displaystyle \frac{v_{\infty }}{\tau }}{\displaystyle \frac{{\partial }^{2}w}{\partial x\partial t}}+{\displaystyle \frac{v_{\infty }^2}{\ell }}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}\right).$$

(9)

Collecting terms that have the same derivative of w, we obtain

$${\displaystyle \frac{1}{{\tau }^2}}\left(m+m_a\right){\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}+2{\displaystyle \frac{V_0}{\ell \tau }}\left(m+m_a{r}_v\right){\displaystyle \frac{{\partial }^{2}w}{\partial x\partial t}}+{\displaystyle \frac{V_0^2}{{\ell }^2}}\left(m+m_a{r}_v^2\right){\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}=0,$$

(10)

where we have defined

$$m_a={\displaystyle \frac{\pi \ell {\rho }_f}{4}},r_v={\displaystyle \frac{v_{\infty }}{V_0}}.$$

(11)

The expressions (10) and (11) reduce to a classical one- or three-term single-parameter added-mass model by choosing $r_v=0$ (i.e., $v_{\infty }=0$, no free-stream flow in laboratory coordinates) or $r_v=1$ (i.e., $v_{\infty }=V_0$, the whole air mass moves with the panel).

3. Non-Dimensional Initial Boundary Value Problem

Using (10) and (11), the governing equation for dynamical behavior can be written in non-dimensional form as

$${\displaystyle \frac{1}{{\tau }^2}}\left(m+m_a\right){\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}+2{\displaystyle \frac{V_0}{\ell \tau }}\left(m+m_a{r}_v\right){\displaystyle \frac{{\partial }^{2}w}{\partial x\partial t}}+{\displaystyle \frac{V_0^2}{{\ell }^2}}\left(m+m_a{r}_v^2\right){\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}$$

$$-T{\displaystyle \frac{1}{{\ell }^2}}{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+D{\displaystyle \frac{1}{{\ell }^4}}{\displaystyle \frac{{\partial }^{4}w}{\partial x^4}}={\displaystyle \frac{1}{\ell }}g(x,t).$$

(12)

Define

$$C=\sqrt{T/m,}$$

(13)

which is the critical velocity of a traveling string or membrane with no bending rigidity. Let us define also the non-dimensional quantities

$$\gamma ={\displaystyle \frac{\ell }{m}}{\rho }_f,\beta ={\displaystyle \frac{D}{m{C}^2{\ell }^2}}={\displaystyle \frac{D}{{\ell }^{2}T}},$$

$$c={\displaystyle \frac{V_0}{C}},\kappa ={\displaystyle \frac{v_{\infty }}{C}},\alpha ={\displaystyle \frac{\ell }{\tau C}},r_m={\displaystyle \frac{m_a}{m}}.$$

(14)

We will also use the non-dimensional external load

$$g^{\prime }\left(x^{\prime },t^{\prime }\right)={\displaystyle \frac{\ell }{m{C}^2}}g\left(\ell x^{\prime },\tau t^{\prime }\right)={\displaystyle \frac{\ell }{T}}g(\ell x^{\prime },\tau t^{\prime }),$$

(15)

where g is the original dimensional load function. Then, we multiply Equation (12) by the factor ${\ell }^2/m{C}^2$, and use the definitions (13)–(15). We obtain

$${\alpha }^2\left(1+r_m\right){\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}+2\alpha c\left(1+r_v{r}_m\right){\displaystyle \frac{{\partial }^{2}w}{\partial x\partial t}}+\left(c^2-1+c^2{r}_v^2{r}_m\right){\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+\beta {\displaystyle \frac{{\partial }^{4}w}{\partial x^4}}=g(x,t),$$

(16)

where all the variables and functions are in non-dimensional form, and the domain is

$$\Omega =\left\{-1<x<1,t>0\right\}.$$

We will use Equation (16) to determine the dynamical behavior and static loss of stability. As for the boundary conditions, we will use the simply supported conditions:

$${\left(w\right)}_{x=-1}={\left(w\right)}_{x=1}=0,$$

(17)

$${\left({\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}\right)}_{x=-1}={\left({\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}\right)}_{x=1}=0.$$

(18)

A direct simulation of dynamical behavior requires also two initial conditions

$${\left(w\right)}_{t=0}=g_1\left(x\right),$$

(19)

$${\left({\displaystyle \frac{\partial w}{\partial t}}\right)}_{t=0}=g_2\left(x\right),$$

(20)

where $g_1\left(x\right)$ and $g_2\left(x\right)$ are known. For determining free vibration modes, no initial conditions are needed.

4. Static Problem of Aerothermoelastic Stability

In the context of static analysis, we set

$${\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}(x,t)=0,{\displaystyle \frac{{\partial }^{2}w}{\partial x\partial t}}=0,g(x,t)=0$$

(21)

in Equation (16).

For the static stability problem, we use the Euler approach of determining a nontrivial steady-state solution, and the associated critical values of the problem parameters of interest. We concentrate on finding the critical panel velocity and critical temperature, which play a role analogous to Euler’s critical compression force for the axially compressed beam.

Using (21) in (16), the steady-state equation is

$$\left[c^2\left(1+r_m{r}_v^2\right)-1\right]{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+\beta {\displaystyle \frac{{\partial }^{4}w}{\partial x^4}}=0,$$

(22)

with boundary conditions (19) and (20). Equation (22) describes the buckling of a traveling panel submerged in potential flow.

It is convenient to transform the original problem (22) with boundary conditions (17) and (18) into a sequence of two problems. Let us introduce a new variable $\psi$:

$$\psi \left(x\right)={\displaystyle \frac{{\partial }^{2}w}{\partial x^2}},-1\le x\le 1.$$

(23)

This transforms Equation (22) into

$${\displaystyle \frac{{\partial }^2\psi }{\partial x^2}}+\lambda \psi =0,-1\le x\le 1,$$

(24)

$$\psi (-1)=0,\psi \left(1\right)=0,$$

(25)

where the introduced eigenvalue $\lambda$ (the load parameter of the stability problem) is

$$\lambda ={\displaystyle \frac{1}{\beta }}\left[c^2\left(1+r_m{r}_v^2\right)-1\right].$$

(26)

The solution of (24) is

$$\psi \left(x\right)=A_{1}sin\left(\sqrt{\lambda }\left({\displaystyle \frac{x+1}{2}}\right)\right)+A_2\left(cos\left(\sqrt{\lambda }\left({\displaystyle \frac{x+1}{2}}\right)\right)\right),$$

where $A_1$ and $A_2$ are arbitrary constants. With the help of boundary conditions (25), the solution is reduced to

$$\psi \left(x\right)=A_{1}sin\left(\sqrt{\lambda }\left({\displaystyle \frac{x+1}{2}}\right)\right),$$

(27)

$$\lambda ={\left({\displaystyle \frac{j\pi }{2}}\right)}^2,j=1,2,\cdots$$

(28)

The minimal value of $\lambda$ corresponds to $j=1$, and consequently ${\lambda }_{min}={\pi }^2/4$. Thus, we have the following basic relation for the critical parameters of static loss of stability:

$${\displaystyle \frac{{\pi }^2}{4}}={\displaystyle \frac{{\ell }^2}{D}}\left(m{V}_0^2+m_a{v}_{\infty }^2-T\right)={\displaystyle \frac{{\ell }^2}{D}}\left(m{V}_0^2+m_a{v}_{\infty }^2-T_0+{\displaystyle \frac{Eh}{1-\nu }}{\epsilon }_0\right).$$

(29)

As follows from (29) and the expressions for ${\epsilon }_{\theta }$ and $m_a$, the critical velocity of static loss of stability (divergence, buckling) and the critical temperature are expressed as

$${\left(V_0^2\right)}_{\mathrm{div}}={\displaystyle \frac{1}{m}}\left({\displaystyle \frac{{\pi }^{2}D}{4{\ell }^2}}+T_0-{\displaystyle \frac{Eh}{1-\nu }}{\alpha }_0\theta -{\displaystyle \frac{\pi \ell }{4}}{\rho }_f{v}_{\infty }^2\right),$$

(30)

$${\theta }_{\mathrm{div}}={\displaystyle \frac{1-\nu }{Eh{\alpha }_{\theta }}}\left({\displaystyle \frac{{\pi }^{2}D}{4{\ell }^2}}+T_0-m{V}_0^2-{\displaystyle \frac{\pi \ell }{4}}{\rho }_f{v}_{\infty }^2\right).$$

(31)

The influence of the added mass of the fluid on the critical velocity is represented by

$${\left({\displaystyle \frac{V_0}{v_{\infty }}}\right)}_{\mathrm{div}}=\sqrt{K-{\displaystyle \frac{m_a}{m}}},K={\displaystyle \frac{1}{v_{\infty }}}\left({\displaystyle \frac{{\pi }^{2}D}{4m{\ell }^2}}+{\displaystyle \frac{T}{m}}\right),$$

(32)

and is shown in Figure 2 by the solid lines corresponding to $K=1,2,3,4$.

5. Dynamic Problem of Free Vibrations

For numerical analysis, a discrete approximation will be used for the partial differential Equation (16) with boundary conditions (17) and (18) and initial conditions (19) and (20). In this context, Galerkin methods are especially convenient and well-known techniques to solve the problem. We will space-discretize using the finite element method with $C^2$ continuous Hermite elements.

To characterize the dynamics of the system, let us perform a numerical dynamic stability analysis. Following [37], we set the load as $g\equiv 0$, and formulate an eigenvalue problem, using the time-harmonic trial function

$$w(x,t)=exp\left(st\right)W\left(x\right),$$

(33)

where s is the stability exponent (a complex number) and $W\left(x\right)$ is the vibration mode. We will solve the problem for eigenvalue–eigenfunction pairs $(s,W)$. Loss of stability occurs at such values of the axial drive velocity $V_0$, where at least one eigenvalue s transitions to the positive half-plane (i.e., where the real part $\mathrm{Re}s$ becomes positive). The critical velocity is the smallest positive $V_0$ such that stability is lost.

In the original continuum problem, there will be a countably infinite number of solutions, corresponding to the spectrum of vibration modes. When the problem is discretized, if there are N degrees of freedom in the discretization, it will have $2N$ solutions because, upon inserting (33) into (16), the resulting ordinary differential equation will be a quadratic polynomial in s. This results in a quadratic eigenvalue problem, which can be converted into a twice larger generalized linear eigenvalue problem by the companion form technique, see [38]. However, not all of the numerical solutions will be solutions of the continuum problem. Attempting to represent an infinite spectrum using a finite basis leads to an aliasing phenomenon. Best accuracy will be obtained for the first few of the lowest modes (slowest vibration, smallest imaginary part $\left|\mathrm{Im}s\right|$). Below, we will investigate the four lowest modes.

Because the material is fully elastic, and hence the model ignores dissipation, then if s is an eigenvalue, $-s$ is also. This is in addition to the general symmetry (also valid for viscoelastic models) that, if s is an eigenvalue, then also $\mathrm{conj}\left(s\right)$ is. This latter property only reflects the physical symmetry with respect to the x axis ($z=0$) in the problem setup.

We then represent the displacement as a Galerkin series:

$$W\left(x\right)=\sum _{n=1}^{\infty }{c}_n{\mathrm{\Psi }}_n\left(x\right),$$

(34)

where ${\mathrm{\Psi }}_n$ are the global shape functions, and $c_n$ are the global degrees of freedom. In practice, we will truncate the Galerkin series at a finite value of n (determined by the mesh).

The basis functions ${\mathrm{\Psi }}_n\left(x\right)$ are defined piecewise, using an affine coordinate transformation $x(\xi ,n)=a_n\xi +b_n$, where n is the index of the global element, and $a_n$, $b_n$ are the coordinate transformation coefficients that map the reference element $\xi \in [0,1]$ onto the nth global element. The local shape functions $N_j\left(\xi \right)$ on the reference element $\xi \in [0,1]$ are defined as follows:

$$\begin{array}{cc}\hfill N_1\left(\xi \right)& =-6{\xi }^5+15{\xi }^4-10{\xi }^3+1,\hfill \\ \hfill N_2\left(\xi \right)& =-3{\xi }^5+8{\xi }^4-6{\xi }^3+\xi ,\hfill \\ \hfill N_3\left(\xi \right)& =\frac{1}{2}\left(-{\xi }^5+3{\xi }^4-3{\xi }^3+{\xi }^2\right),\hfill \\ \hfill N_4\left(\xi \right)& =6{\xi }^5-15{\xi }^4+10{\xi }^3,\hfill \\ \hfill N_5\left(\xi \right)& =-3{\xi }^5+7{\xi }^4-4{\xi }^3,\hfill \\ \hfill N_6\left(\xi \right)& =\frac{1}{2}\left({\xi }^5-2{\xi }^4+{\xi }^3\right).\hfill \end{array}$$

These functions interpolate, respectively, ${w|}_{\xi =0}$, ${\partial w/\partial \xi |}_{\xi =0}$, ${\partial }^{2}w/\partial x^2{|}_{\xi =0}$, ${w|}_{\xi =1}$, ${\partial w/\partial \xi |}_{\xi =1}$, and ${\partial }^{2}w/\partial x^2{|}_{\xi =1}$. These are then copied and scaled via the coordinate mapping to each global element.

Compared to standard $C^0$ continuous FEM, one must be careful when dealing with the values of the derivative degrees of freedom (accounting for the scaling by $a_n$ due to the coordinate transformation); the rest of the discretization is performed exactly as in standard FEM.

This basis is especially suitable for paper materials because it works also for the study of viscoelastic materials of the Kelvin–Voigt type, where the interaction between the axial motion and the viscoelasticity introduces the fifth space derivative in the strong form of the governing equation. With $C^2$ continuous elements, up to sixth-order equations can be treated, at the cost of being unable to approximate solutions having only $C^0$ or $C^1$ continuity. For a free-vibration solution, which in the case of constant coefficients is smooth, this is not an issue in practice, so we may use this basis for our problem, although strictly speaking the solution of the elastic case lives in $H^2(\Omega )$ (of which $C^{1,ϵ}\left(\overline{\Omega }\right)$ is a subset, in one space dimension; see the embedding theorems in a book by Adams [39]).

Let us now develop the discrete eigenvalue problem. Inserting the time-harmonic trial function (33) into (16) and setting $g\equiv 0$, we have

$$s^2{\alpha }^2\left(1+r_m\right)W+2s\alpha c\left(1+r_{\nu }{r}_m\right){\displaystyle \frac{\partial W}{\partial x}}+\left(c^2-1+c^2{r}_{\nu }^2{r}_m\right){\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}+\beta {\displaystyle \frac{{\partial }^{4}W}{\partial x^4}}=0.$$

(35)

To obtain the weak form, we multiply (35) by the test function ${\mathrm{\Psi }}_j\left(x\right)$, and integrate over the non-dimensional space domain $\Omega =\left\{x\in (-1,1)\right\}$:

$$s^2{\alpha }^2\left(1+r_m\right){\int }_{-1}^{1}W{\mathrm{\Psi }}_j\mathrm{d}x+2s\alpha c(1+r_{\nu }{r}_m){\int }_{-1}^1{\displaystyle \frac{\partial W}{\partial x}}{\mathrm{\Psi }}_j\mathrm{d}x$$

$$+\left(c^2-1+c^2{r}_{\nu }^2{r}_m\right){\int }_{-1}^1{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}{\mathrm{\Psi }}_j\mathrm{d}x+\beta {\int }_{-1}^1{\displaystyle \frac{{\partial }^{4}W}{\partial x^4}}{\mathrm{\Psi }}_j\mathrm{d}x=0.$$

Applying integration by parts (twice in the fourth-order term) and the simply supported boundary conditions (17) and (18), we have the weak form

$$\begin{array}{c}{s}^2{\alpha }^2\left(1+r_m\right){\int }_{-1}^{1}W{\mathrm{\Psi }}_j\mathrm{d}x+2s\alpha c(1+r_{\nu }{r}_m){\int }_{-1}^1{\displaystyle \frac{\partial W}{\partial x}}{\mathrm{\Psi }}_j\mathrm{d}x\hfill \\ \hfill -\left(c^2-1+c^2{r}_{\nu }^2{r}_m\right){\int }_{-1}^1{\displaystyle \frac{\partial W}{\partial x}}\frac{\partial {\mathrm{\Psi }}_j}{\partial x}\mathrm{d}x+\beta {\int }_{-1}^1{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}\frac{{\partial }^2{\mathrm{\Psi }}_j}{\partial x^2}\mathrm{d}x=0.\end{array}$$

(36)

Inserting the Galerkin series (34) into the weak form (36) gives us the following system:

$$\sum _{n=1}^{n_{max}}\left\{s^2{\alpha }^2(1+r_m)A_{jn}+\right.$$

$$\left.+2s\alpha c(1+r_v{r}_m)B_{jn}+\left[-\left(c^2-1+c^2{r}_v^2{r}_m\right)C_{jn}+\beta D_{jn}\right]\right\}{c}_n=0.$$

(37)

Here, $j,n=1,2,3,\cdots$, and the matrices $A_{jn}$, $B_{jn}$, $C_{jn}$, and $D_{jn}$ are defined by

$$A_{jn}={\int }_{-1}^1{\mathrm{\Psi }}_n\left(x\right){\mathrm{\Psi }}_j\left(x\right)\mathrm{d}x,$$

(38)

$$B_{jn}={\int }_{-1}^1{\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_n\left(x\right)}{\mathrm{d}x}}{\mathrm{\Psi }}_j\left(x\right)\mathrm{d}x,$$

(39)

$$C_{jn}={\int }_{-1}^1{\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_n\left(x\right)}{\mathrm{d}x}}{\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_j\left(x\right)}{\mathrm{d}x}}\mathrm{d}x,$$

(40)

$$D_{jn}={\int }_{-1}^1{\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\Psi }}_n\left(x\right)}{\mathrm{d}{x}^2}}{\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\Psi }}_j\left(x\right)}{\mathrm{d}{x}^2}}\mathrm{d}x.$$

(41)

Further defining

$$\begin{array}{}\mathrm{(42)}& M_2\hfill & ={\alpha }^2(1+r_m)A_{jn},\hfill \mathrm{(43)}& M_1\hfill & =2\alpha c(1+r_v{r}_m)B_{jn},\hfill \mathrm{(44)}& M_0\hfill & =-\left(c^2-1+c^2{r}_v^2{r}_m\right)C_{jn}+\beta D_{jn},\hfill \end{array}$$

we obtain the (discrete) quadratic eigenvalue problem for eigenvalue–eigenvector pairs $(s,c)$:

$$\left\{M_2{s}^2+M_{1}s+M_0\right\}c=0.$$

(45)

Its companion form is the following twice larger generalized linear eigenvalue problem:

$$\left[\begin{array}{cc}-M_1& -M_0\\ I& 0\end{array}\right]\left[\begin{array}{c}\hfill sc\\ \hfill c\end{array}\right]=\left[\begin{array}{cc}{M}_2& 0\\ 0& I\end{array}\right]s\left[\begin{array}{c}\hfill sc\\ \hfill c\end{array}\right].$$

(46)

Equation (46) can be solved using a standard solver.

Problem parameters values used in the numerical examples, corresponding to paper materials, are shown in

Table 1. Problem parameters used in the numerical examples.

${\rho }_{\mathrm{f}}$ [$\mathrm{kg}/{\mathrm{m}}^3$]	$1.225$
$v_{\infty }$ [$\mathrm{m}/\mathrm{s}$]	0
$T_0$ [$\mathrm{N}/\mathrm{m}$]	500
m [$\mathrm{kg}/{\mathrm{m}}^2$]	$0.08$
ℓ [$\mathrm{m}$]	1
h [$\mathrm{m}$]	${10}^{-4}$
E [$\mathrm{N}/{\mathrm{m}}^2$]	${10}^9$
$\nu$ [1]	$0.3$
$\tau$ [$\mathrm{s}$]	$\ell /\sqrt{T/m}$
$\alpha$ [$1/\mathrm{K}$]	$3\times {10}^{-6}$

. The value of the linear thermal expansion coefficient has been chosen from the range typical for paper materials according to Kouko [40]. Specifically, it is valid for the machine direction (i.e., x in our coordinate system).

Consider a setup where the temperature is kept constant, and the value of $\alpha$ is a constant ${\alpha }_0$ within the material. Equation (2) then gives us

$$\begin{array}{cc}\hfill T_{\theta }={\displaystyle \frac{Eh}{1-\nu }}{\epsilon }_{\theta }={\displaystyle \frac{Eh}{1-\nu }}{\alpha }_0\theta =\frac{{10}^9\mathrm{N}/{\mathrm{m}}^2\times {10}^{-4}\mathrm{m}}{1-0.3}3\times & {10}^{-6}\frac{1}{\mathrm{K}}\theta \hfill \\ & =\left(\frac{3}{7}\frac{\mathrm{N}/\mathrm{m}}{\mathrm{K}}\right)\theta \approx \left(0.42857\frac{\mathrm{N}/\mathrm{m}}{\mathrm{K}}\right)\theta .\hfill \end{array}$$

(47)

The effect is rather small, but, nevertheless, if the system is already operating near its stability limit, thermal expansion may reduce the critical velocity just enough to make the system lose stability. As an example, let us consider a rather extreme temperature difference, between $20{}^{\circ }\mathrm{C}$ and $90{}^{\circ }\mathrm{C}$, which may occur when a paper web initially at room temperature enters the dryer section. The temperature difference is $\theta =70\mathrm{K}$, which yields $T_{\theta }=(3/7)\times 70\mathrm{N}/\mathrm{m}=30\mathrm{N}/\mathrm{m}$. With a typical level of applied axial tension, $T_0=500\mathrm{N}/\mathrm{m}$, the decrease in tension due to thermal expansion, from $500\mathrm{N}/\mathrm{m}$ to $470\mathrm{N}/\mathrm{m}$, is $6\%$.

In the reference setup, where $\theta =0$ and $T=T_0$, the four pairs of s with the smallest magnitude, solved from Equation (46), are shown in Figure 3. Critical values of the non-dimensional axial drive velocity c are listed in

Table 2. Critical points, where the real part of at least one solution enters or exits the positive half-plane. No axial stream, $v_{\infty }=0$.

$\mathit{c}-1$ [${10}^{-4}$]	$\mathit{c}-1$ [${10}^{-4}$]	Mode	Type	Note
with $\mathit{\theta }=0$	with $\mathit{\theta }=70\mathbf{K}$
0.009	0.009	1	Loss of stability	Classic divergence gap.
0.036	0.038	1	Gain of stability
0.051	0.054	$1+2$	Loss of stability
0.081	0.086	$1+2$	Gain of stability	Other branch still unstable.
0.145	0.153	$1+2$	Gain of stability	End of unstable branch.
0.145	0.153	$2+3$	Loss of stability
0.226	0.240	$2+3$	Gain of stability	Other branch still unstable.
0.285	0.303	$3+4$	Loss of stability
0.325	0.346	$2+3$	Gain of stability	End of unstable branch.

. Based on the table, the lowest critical velocity is

$$\begin{array}{cc}\hfill V_0^{\mathrm{cr}}& =c^{\mathrm{cr}}C=(1+0.009\times {10}^{-4})\sqrt{T/m}\hfill \\ \hfill & =(1+0.009\times {10}^{-4})\sqrt{\frac{500\mathrm{N}/\mathrm{m}}{0.08\mathrm{kg}/{\mathrm{m}}^2}}\approx 79.05701\mathrm{m}/\mathrm{s}.\hfill \end{array}$$

In the thermally expanded case, when $\theta =70\mathrm{K}$, we have $T=470\mathrm{N}/\mathrm{m}$. Visually, the solution is not significantly different from Figure 3, so the graph is omitted. The critical points are listed in Table 2. The lowest critical velocity is

$$\begin{array}{cc}\hfill V_0^{\mathrm{cr}}& =c^{\mathrm{cr}}C=(1+0.009\times {10}^{-4})\sqrt{T/m}\hfill \\ \hfill & =(1+0.009\times {10}^{-4})\sqrt{\frac{470\mathrm{N}/\mathrm{m}}{0.08\mathrm{kg}/{\mathrm{m}}^2}}\approx 76.64862\mathrm{m}/\mathrm{s},\hfill \end{array}$$

showing a decrease of $3\%$ with regard to the reference setup. Note the different value of the normalization constant C because the tension is now lower.

In both cases, the presence of bending rigidity ($D>0$) pushes the first critical point above $c=1$, but only by a small amount because, for paper materials, D is small. From Equations (2) and (14), we find

$$\begin{array}{cc}\hfill D& ={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}}=\frac{{10}^9\mathrm{N}/{\mathrm{m}}^2\times {\left({10}^{-4}\mathrm{m}\right)}^3}{12(1-{0.3}^2)}=9.1575\times {10}^{-5}\mathrm{N}\mathrm{m},\hfill \\ \hfill \beta & =\frac{D}{{\ell }^{2}T}=\frac{9.1575\times {10}^{-5}\mathrm{N}\mathrm{m}}{{\left(1\mathrm{m}\right)}^2\times 500\mathrm{N}/\mathrm{m}}=1.8315\times {10}^{-7}.\hfill \end{array}$$

Nevertheless, because $\beta$ is the coefficient of the highest-order term in (16), the presence of finite bending rigidity, no matter how small, introduces a singular perturbation to the equation, changing its qualitative behavior (on singular perturbations, see, e.g., [41]).

Finally, for comparison, Figure 4 and

Table 3. Critical points, where the real part of at least one solution enters or exits the positive half-plane. Air mass moves with the panel, $v_{\infty }=V_0$. Note the scale of c.

$\mathit{c}-1$ [${10}^{-1}$]	$\mathit{c}-1$ [${10}^{-1}$]	Mode	Type	Note
with Type $\mathit{\theta }=0$	with $\mathit{\theta }=70\mathbf{K}$
$-0.391490$	$-0.391490$	1	Loss of stability	Classic divergence gap.
$-0.391464$	$-0.391462$	1	Gain of stability
$-0.391450$	$-0.391447$	$1+2$	Loss of stability
$-0.391421$	$-0.391416$	$1+2$	Gain of stability	Other branch still unstable.
$-0.391361$	$-0.391352$	$2+3$	Loss of stability
$-0.391360$	$-0.391351$	$1+2$	Gain of stability	End of unstable branch.
$-0.391282$	$-0.391268$	$2+3$	Gain of stability	Other branch still unstable.
$-0.391225$	$-0.391208$	$3+4$	Loss of stability
$-0.391187$	$-0.391167$	$2+3$	Gain of stability	End of unstable branch.

show the solution for $\theta =0$ and $T=T_0$, but with $v_{\infty }=V_0$, leading to $r_v=1$. This means that the air mass moves axially with the panel. It is seen that the shape of the solution remains the same, but the critical velocity is decreased to $61\%$ of its value in the previous case where $v_{\infty }=0$.

6. Conclusions

An effective partly-analytical and partly-numerical approach has been presented to model a thermoelastic panel moving at a constant velocity, subjected to potential flow. An aerothermoelastic model based on a homogeneous temperature distribution and an added-mass approximation for the fluid reaction pressure was developed and analyzed. This model was then applied to some particular cases of stationary and dynamic interaction. As a result, the stationary eigenvalue problem of static loss of stability (divergence, buckling) was solved analytically, by employing spectral analysis. The critical values of the problem parameters were found. In the case of nonstationary behavior of the panel, the analysis was performed using the Galerkin method. The presented approach has applications, for example, in elasticity, aeroelasticity, aerothermoelasticity and axially moving materials.