Numerical simulations of a mixed finite element method for damped plate vibration problems

1.   Introduction

As a basic structural unit, plates are widely used in many places, such as spacecrafts and aircrafts, ships, buildings, containers, etc. The vibration of plates caused by external forces can lead to serious damage to the entire structure of the machinery or building. One way to reduce the damage caused by vibration is by applying the viscous damping strategy. The vibration of damped plates is described by fourth-order differential equations, whose analytical solutions are often excessively difficult to obtain. Thus, the theoretical analysis and numerical calculation of the vibration of damped plates are of great research interest.

So far, a great number of studies have been conducted on the vibration problems of damped plates. Leissa et al. studied the free vibration of rectangular plates ^[1] and the vibrations of cantilevered shallow cylindrical shells of rectangular platforms ^[2]. Nair et al. discussed the quasi-degeneracies in plate vibration problems ^[3]. Wang et al. studied the vibration problems of flexible circular plates with initial deflection^[4]. Hui Li et al. studied the vibration of foam core^[5,6], considered the nonlinear vibration analysis of fiber metal laminated plates with multiple viscoelastic layers^[7] and considered the vibration damping of multifunctional grille composite sandwich plates^[8,9].

The numerical methods studied for the plate vibration problems include the integration method, finite difference method, finite element method, mixed finite element method, etc. For example, Rock et al. used the finite element method in the study of the free vibration and dynamic response of thick and thin plates ^[10]. Bezine proposed a mixed boundary integral as a finite element approach to plate vibration problems ^[11]. Qian et al. studied the vibration characteristics of cracked plates ^[12]. Xu et al. analyzed the vibration problems of thin plates using the integral equation method ^[13]. Duran et al. conducted the finite element analysis of the vibration problem of a plate coupled with a fluid ^[14]. Xiong et al. conducted an analysis of free vibration problems for a thin plate by the local Petrov-Galerkin method ^[15]. Dawe discussed a finite element approach to plate vibration problems ^[16]. Wu et al. utilized the mesh-free least-squares-based finite difference method for large-amplitude free vibration analysis of arbitrarily shaped thin plates ^[17]. Mora et al. analyzed the buckling and the vibration problems of thin plates using a piecewise linear finite element method ^[18]. Werfalli et al. analyzed the vibration of rectangular plates using Galerkin-based finite element method ^[19]. Yang et al. discussed a differential quadrature hierarchical finite element method and its application to thin plate free vibration ^[20]. The mixed finite element method is effective in solving differential equations. The general theory of this method was established by Brezzi and Babuska in 1970s to solve second order elliptic problems ^[21,22].

Later, Brezzi et al. used the mixed finite element method to solve second order elliptic problems in three variables ^[23]. Diegel et al. discussed the stability and convergence of a second order mixed finite element method for the Cahn-Hilliard equation ^[24]. Singh et al. performed the compositional flow modeling using a multi-point flux mixed finite element method ^[25]. Burger et al. studied a mixed finite element method for nonlinear diffusion equations ^[26].

The mixed finite element method is also effective in simulating fourth-order differential equations, including both biharmonic equations and vibration equations. For biharmonic equations, Monk et al. utilized a stabilized mixed finite element method for the biharmonic equation based on biorthogonal systems ^[27]. Stein et al. proposed a mixed finite element method with piecewise linear elements for the biharmonic equation on surfaces ^[28]. Meng et al. studied the optimal order convergence for the lowest order mixed finite element method of the biharmonic eigenvalue problem ^[29]. For vibration equations, Meng et al. studied a mixed virtual element method for the vibration problem of clamped Kirchhoff plate ^[30].

As far as we know, the current literature lacks studies that utilize the mixed finite element method to solve vibration equations for viscously damped plates. Therefore, this work seeks to establish the mixed finite element scheme for the initial boundary conditions of damped plate vibration problems and to verify the existence and uniqueness of the approximate solution for the semi-discrete and backward Euler fully discrete schemes. An error analysis is conducted, and numerical case studies are conducted to validate the effectiveness and precision of the mixed finite element scheme, as well as to quantify the influence of the damping coefficient on plate vibrations.

According to the theory of elasticity, there is a vibration equation of thin plate in ^[31],

In this article, we add the damped term and consider the damped plate vibration problem:

Where $D$ is the flexural rigidity, $m = \rho h$ is the mass per unit area, $\rho$ is the mass density of the plate, and $h$ is the thickness of the plate. $\lambda$ is the damping factor, $f$ is the smooth function, $w(x, y, t)$ is the flexible surface function, $\Omega$ is the piecewise smooth bounded polygon region, $(0, T]$ is the time interval, $\Psi(x, y), \Phi(x, y)$ are known functions.

In this paper, the damping plate vibration equation is analyzed by the mixed finite element method. The advantage of the mixed finite element method lies in its ability to reduce the order of the high order differential equations by introducing intermediate variables, which often have physical meaning by themselves. Consequently, it can reduce the requirement for smoothness of the finite element space and hence simplify the interpolation space of the finite elements. Moreover, by using the mixed finite element method, both the unknown variables and the intermediate variables with realistic meaning can be obtained, hence increasing the precision of the discrete solutions. Compared to other methods, the mixed finite element method is easier to apply and more likely to yield meaningful solutions.

This article is divided into five sections. The first section introduces the research background of the plate vibration problems. The second section provides the variational formulation for the initial boundary conditions of damped plate vibration problems. The third and the fourth sections discuss the construction of the semi-discrete and fully discrete mixed finite element schemes for the initial boundary conditions of the damped plate vibration problems, respectively, followed by the verification of the existence and uniqueness of such schemes and the error analyses. Finally, the fifth section presents the numerical case studies aimed at validating the theoretical discussions in the previous sections.

2.   Variational format

Introducing auxiliary variables $\Delta w = -u, v = w_t$, where $\Delta = \frac{\partial}{\partial x^{2}}+\frac{\partial}{\partial y^{2}}$, we first rewrite Eq.$(1.1)$ into the following coupled system:

Multiplying both sides of $(2.1)(a)$ by $\varphi\in H_0^1(\Omega)$ and using Green's formula, we have

Multiplying both sides of $(2.1)(b)$ by $\psi\in H_0^1(\Omega)$ and using Green's formula, we obtain

Therefore, we have the following mixed weak formulation of $(2.1)$ : find $\{u, v\}:[0, T]\rightarrow H_0^1(\Omega)\times H_0^1(\Omega)$, such that

3.   Semi-discrete finite element scheme

First, we define the finite element space. Let $\Omega$ be a rectangular region whose boundaries are parallel to the two axes. The region $\Omega$ is divided into regular triangulation. $\jmath_h$ is a triangulation family whose region satisfies the regular hypothesis, $K$ represents the triangulation unit, and $h$ is the maximum diameter of the subdivision unit. $\Omega = \bigcup\limits_{K\epsilon\jmath_h}K$, $S_h = \{v_h\mid v_h\mid_K\in P_k(K), \forall K\epsilon\jmath_h\}\subset H^1(\Omega)$ is the finite element space composed of piecewise linear degree polynomials on $\jmath_h$.Then, the corresponding semi-discrete finite element scheme of $(2.2)$ is to find $\{u_h, v_h\}:[0, T]\rightarrow S_h^0\times S_h^0$, $S_h^0 = S_h\cap H_0^1(\Omega)$, such that

$R_h$ is an elliptic projection operator, which will be given below. The existence and uniqueness of semi-discrete finite element approximation scheme solutions and error analysis are given below.

Theorem 3.1. Existence and uniqueness of the solution of the semi-discrete finite element approximation scheme $(3.1)$.

Proof. $\{ \phi_i \} _{i = 1}^M$ be a set of bases of $S_h^0$. Then $u_h$ = $\sum\limits_{j = 1}^{M}u_j\phi_j$, $v_h$ = $\sum\limits_{j = 1}^{M}v_j\phi_j$. According to $(3.1)(a)$ and $(3.1)(b)$, we have the following equalities

Where $U(t) = (u_1(t), u_2(t), \cdots, u_N(t))^T,$

$V(t) = (v_1(t), v_2(t), \cdots, v_N(t))^T,$ $A = (\nabla\phi_j, \nabla\phi_i),$

$B = (\phi_j, \phi_i),$ $F = (f, \phi_i)$.

According to $(3.3)$, we deduce that

Substituting $(3.4)$ into $(3.2)$, we arrive at

$U(0)$ can be determined by $u_h(x, y, 0)$, and $(3.5)$ is an ordinary differential equation about vector ${U(t)}$. $A, \; BA^{-1}B$ are symmetric positive definite matrices. According to the theory of ordinary differential equations, it is easy to know that the solution of the semi-discrete finite element approximation scheme is existent and unique.

In the following discussion, we will derive the proof of the error estimates for semi-discrete schemes. For carrying out an analysis, we need to introduce a useful lemma. First, to give the error analysis, for $\forall 0\leq t\leq T$, we consider the elliptic projection operator $R_h:H_0^1\rightarrow S_h^0$ such that $(\nabla(u-R_hu), \nabla v_h) = 0, \; \forall v_h\in S_h^0$, which leads to the following estimate inequality.

Lemma 3.1. ^[32] $\forall u\in H_0^{k+1}$, such that

Corollary 3.1. $\forall u\in H_0^{k+1}$, such that

Lemma 3.2. ^[33] The family $S_h$ is based on a family of quasiuniform triangulations $\jmath_h$ and $S_h$ consists of piecewise polynomials of degree at most $k-1$, and then one may show the inverse inequality:

In the next analysis, we will discuss the proof of semi-discrete error estimates based on the elliptic projection in detail.

Theorem 3.2. Let $\{ u, v \}$ and $\{ u_h, v_h \}$ be the solutions of $(2.2)(a)\ and \ (2.2)(b)$ and $(3.1)(a) \ and\ (3.1)(b)$, respectively, we have $L^2$-mode and $H^1$-mode error estimations of variable $\{ u, v \}$:

Proof. To simplify, we now rewrite the errors as $u-u_h = (u-R_hu)+(R_hu-u_h) = \rho+\theta, \; v-v_h = (v-R_hv)+(R_hv-v_h) = \eta+\xi$.

$\forall \varphi_h, \psi_h\in S_h^0$, subtracting $(2.2)(a)$ from $(3.1)(a)$, subtracting $(2.2)(b)$ from $(3.1)(b)$, and applying the elliptic projection operator, we have the error equation:

Choosing $\varphi_h = \xi, \psi_h = \theta$, add $(3.13)$ and $D\times(3.14)$, we have

The Young inequality with $\varepsilon$ and corollary $3.1$ being applied to $(3.15)$, we easily obtain

Integrating from $0$ to $t$ on both sides of $(3.16)$, because $\xi(0) = \theta(0) = 0$, we have

We use Gronwall inequality to get

Thus, we have $L^2$-mode error estimation of variable $\{u, v\}$:

Using lemma 3.1 and triangle inequality, we finish the proof of $(3.9) and (3.10)$.

Theorem 3.3. Let $\{u, v\}$ and $\{u_h, v_h\}$ be the solutions of $(2.1)$ and $(2.2)$, respectively. When $\{u, v\}$ is smooth enough, we have the error estimation of variable $\{u_t, v_t\}$:

Similar to Theorem $3.2$, we give a simple proof.

Proof. First, taking the derivative of the variable $t$ of the error equation $(3.13)-(3.14)$, we obtain

Choosing $\varphi_h = \xi_t$ in $(3.20)$, $\psi_h = \theta_t$ in $(3.21)$, we have

With the same method as theorem $3.2$, we easily obtain

Using lemma 3.1 and inverse inequality, we have $H^1$-mode error estimation of variable ${u}$:

In the same way, we have $H^1$-mode error estimation of variable ${v}$:

Hence, we finish the proof of theorem $3.2$.

4.   Full discrete finite element scheme

Let $0 = t_0 < t_1 < \cdots < t_N = T$ be the subdivision of step $\tau = \frac{T}{N}$ in time interval $[0, T]$, $t_n = n\tau, n = 0, 1\cdots, N$, $U^n\in S_h^0$ stand for the approximation of $u(t_n)$, when $t = t_n = n\tau$. For any function $\phi$ on $[0, T]$, define:

Choosing $t = t_n$, we have a format equivalent to $(2.1)$:

Where $R_u^n = \partial_tu^n-u_t^n = \frac{1}{\tau}\int_{t_{n-1}}^{t_n}(t_{n-1}-s)u_{tt}(s)ds, R_v^n = \partial_tv^n-v_t^n = \frac{1}{\tau}\int_{t_{n-1}}^{t_n}(t_{n-1}-s)v_{tt}(s)ds$.

Then, the fully discrete finite element approximation scheme is described as: find $\{U^n, V^n\}$: $[0, T]\rightarrow S_h^0\times S_h^0$, $S_h^0 = S_h\cap H_0^1(\Omega)$, such that

Similarly, we give proof of the existence and uniqueness of the fully discrete finite element scheme solution and error analysis.

Theorem 4.1. Existence and uniqueness of the solution of the fully discrete finite element approximation scheme $(4.2)$.

Proof. Let $\{ \phi_i \} _{i = 1}^M$ be a set of bases of $S_h^0$. We have $U^n = \sum\limits_{i = 1}^{M}u_i^n\phi_i, V^n = \sum\limits_{i = 1}^{M}v_i^n\phi_i$. According to $(4.2)(a)$ and $(4.2)(b)$, we have

where

According to $(4.4)$, we easily arrive at

Substitute $(4.5)$ into $(4.3)$ to obtain

$U^0$ can be determined by $R_hu(x, y, 0)$. $A, BA^{-1}B$ are symmetric positive definite matrices, so the solution of $(4.6)$ is existent and unique, and the solution of $(4.5)$ is existent and unique. The existence and uniqueness of the solution are equivalent to problem $(4.2)(a)$ and $(4.2)(b)$.

Theorem 4.2. Let $\{u^n, v^n\}$ and $\{U^n, V^n\}$ be the solutions of $(4.1)$ and $(4.2)$, respectively, we have $L^2$-mode error estimation of variable $\{u^n, v^n\}$:

Proof. To simplify, we now rewrite the errors as $u^i-U^i = (u^i-R_hu^i)+(R_hu^i-U^i) = \rho^i+\theta^i, v^i-V^i = (v^i-R_hv^i)+(R_hv^i-V^i) = \eta^i+\xi^i$.

$\forall \varphi_h, \psi_h\in S_h^0$, subtracting $(4.1)(a)$ from $(4.2)(a)$, subtracting $(4.1)(b)$ from $(4.2)(b)$, and applying elliptic projection operator, we have the error equation:

Let $\varphi_h = \xi^i, \psi_h = \theta^i$. Adding $(4.8)$ and $D\times(4.9)$, we have

Where

Let's estimate $\sum M_i$ in turn:

Using the Young inequality with $\varepsilon$, lemma $3.1$ and corollary $3.1$, we obtain

Using Cauchy-Schwarz inequality and Young inequality with $\varepsilon$, we have

Substituting them into $(4.10)$, we have

Sum the above formula about $i$ from $1$ to $n$. Noticing that $\xi(0) = \theta(0) = 0$, we have

Using Gronwall Lemma, we have $\tau$ sufficiently small

Thus, we have $L^2$-mode error estimation of variable $\{u^n, v^n\}$:

Using lemma 3.1 and the triangle inequality, we finish the proof of theorem $4.2$.

Next, we give the $H^1$-mode error estimate of $\{u^n, v^n\}$.

Theorem 4.3. Letting $\{u^n, v^n\}$ and $\{U^n, V^n\}$ be the solutions of $(4.1)$ and $(4.2)$, respectively, we have $H^1$-mode error estimation of variable $\{u^n, v^n\}$:

Proof. Choosing $\varphi_h = \theta^i$ in $(4.8)$, we have

Which leads to

The estimate of $M_j$ is as follows:

using Cauchy-Schwarz inequality, Young inequality with $\varepsilon$ and corollary $3.1$, we obtain

Using Cauchy-Schwarz inequality, Young inequality with $\varepsilon$ and Theorem $3.3$, we get

Using Young inequality with $\varepsilon$ and Lemma $3.1$, we have

Using Young inequality with $\varepsilon$, we deduce that

Using Cauchy-Schwarz inequality and Young inequality with $\varepsilon$, we get

Combining $(4.16)-(4.20)$ and using Theorem $4.2$, we have

Using lemma $3.1$ and $(4.21)$, we get

Choosing $\psi_h = \xi^i$ in $(4.9)$, we obtain

Which leads to

We estimate the terms on the right-hand side of $(4.23)$ one by one. Using Cauchy-Schwarz inequality, Young inequality with $\varepsilon$ and corollary $3.1$, we obtain

Using Cauchy-Schwarz inequality, Young inequality with $\varepsilon$ and Theorem $3.3$, we obtain

Combining $(4.24)-(4.26)$ and using Theorem $4.2$, it holds that

Using corollary $3.2$ and $(4.27)$, we have

5.   Numerical experiment

In this section, we provide numerical examples to validate the backward Euler full discretization mixed finite element scheme $(4.2)$ for the vibration problems of damped plates $(2.1)$. We not only validate the convergence order of the error estimate, but also simulate the vibration of damped plates to quantify the influence of damping coefficient on the frequency and amplitude of vibration.

Example $1$

For the numerical calculation, let the space domain be $\Omega = [0, 4]\times[0, 4]$ and let the time domain be $[0, T] = [0, 1]$. Let $D = 1, m = 1, \lambda = 1$. The exact solution to the vibration problem of the damped plate $(2.1)$ is $w$ = cos$t$sin$(\frac{\pi}{4}x)$sin$(\frac{\pi}{4}y)$. The source term $f(x, y, t)$ can be obtained by inserting the given exact solution into the vibration equation $(2.1)$. The mixed finite element space is a double linear first-order polynomial. Keep the time step size $\tau = \frac{1}{100000}$ constant while varying the space step size $h_x = h_y = \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \frac{1}{16}$. Tables 1 and 2 show the space errors and convergence orders, respectively, of the $L^2-norm$ and $H^1-norm$ of the solutions to the backward Euler full discretization mixed finite element scheme $(4.2)$. Keep the space step size $h_x = h_y = \frac{1}{1024}$ constant while varying the time step size $\tau = \frac{1}{4}, \frac{1}{8}, \frac{1}{16}, \frac{1}{32}$. Tables 3 and 4 show the time errors and convergence orders, respectively, of the $L^2-norm$ and $H^1-norm$ of the solutions to the backward Euler full discretization mixed finite element scheme $(4.2)$. The second and third columns in Tables 1 and 2 show the space errors of the $L^2-norm$ and $H^1-norm$ for the solutions to the backward Euler full discretization mixed finite element scheme $(4.2)$, respectively. The fourth and fifth columns show their corresponding space convergence orders. The second and third columns in Tables 3 and 4 show the time errors of the $L^2-norm$ and $H^1-norm$ for the solutions to the backward Euler full discretization mixed finite element scheme $(4.2)$, respectively. The fourth and fifth columns show their corresponding time convergence orders.

The tables illustrate that the space convergence orders are 2 or 1, while the time convergence orders are uniformly 1, for the $L^2-norm$ and $H^1-norm$ of the solutions to the backward Euler full discretization mixed finite element scheme $(4.2)$ for the vibration problems of damped plates $(2.1)$. This is consistent with the theoretical results, and hence the conclusions of the theorem are validated.

When spatial step $h$ = $\frac{1}{1024}$, $w$ = cos$t$sin$(\frac{\pi}{4}x)$sin$(\frac{\pi}{4}y)$, we have

Example $2$

In this numerical example, we not only simulate the vibration of damped plates, but also validate the influence of damping coefficient on the frequency and amplitude of the vibration.

First, let $D = 100, m = 5, \lambda = 40,$ and the external force $f = 0$. Let the non-zero initial displacement of plate vibration be $w$ = sin$(\frac{\pi}{4}x)$sin$(\frac{\pi}{4}y)$. Vibrations at different moments are simulated. The vibration patterns at $t = 0.05, t = 0.2, t = 0.3, t = 1, t = 3 \; {\rm{and}}\; t = 5$ are shown in Figure 1, respectively.

By comparing Figure 1, it is noticed that the amplitude of vibration decreases over time. From $t = 3$, the amplitude changes increasingly slowly until it stabilizes at a fixed value.

Then, let $D = 10, m = 20, \lambda = 40,$ the initial vibration displacement $w = 0$, and the duration of external force = 0.1, i.e.

The change in vibration amplitude over time is studied. The vibration patterns at $t = 0.05, t = 0.1, t = 0.2, t = 0.5, t = 2 \; {\rm{and}}\; t = 2.5$ are shown in Figure 2, respectively.

By comparing Figure 2, it is observed that the amplitude increases from $t = 0.05$ to $t = 0.5$, and then starts to decrease and eventually stabilizes.

Finally, let $D = 10, m = 20$, the initial vibration displacement $w = 0$, and the duration of external force = 0.1, i.e.

The influence on vibration amplitude by changing the damping coefficient is studied. The damping coefficient is set at 10, 20,160 and 640. When $t = 0.15$, the vibration patterns when the damping coefficient is 10, 20,160 and 640 are shown in Figure 3, respectively. The influence of changing the damping coefficient on the vibration frequency is also studied. When the damping coefficient is 10, 20, 40 and 80, the changes of a certain point on the plate as a function of time are shown in Figure 4, respectively.

By comparing Figure 3, it is observed that when the external force is constant, a greater damping coefficient leads to a smaller vibration amplitude. The comparison between Figure 4 suggests that when the external force is constant, a greater damping coefficient leads to a lower frequency.

6.   Conclusions

In this article, we propose the semi-discrete and fully discrete finite element approximation schemes for the vibration equations of damped plates. The existence and the uniqueness of the solution are verified, and the order of convergence of errors is deduced. Moreover, the theoretical analysis is validated by numerical case studies, the pattern of plate vibration is simulated, and the influence of the damping coefficient on the frequency and amplitude of the plate vibration is elucidated. In the future, we attempt to discretize the time using the C-N scheme and approximate the space using elements of higher orders to obtain numerical solutions of higher precision while reducing the calculation load, in order to further improve the simulation of vibration problems of damped plates.

Acknowledgment

The research was supported by the NSFC of China (No. 12171287) and the NSFC of Shandong Province (No. ZR2021MA063).

Conflict of interest

The authors declare that they have no conflicts of interest to this work.