Vibration characteristics of composite plates considering fluid-structure coupling

Abstract

In order to investigate the influence of fluid-structure coupling on the dynamic characteristics of the amphibious aircraft, the free vibration characteristics of composite plates in contact with a bounded fluid are studied based on the combination of Nurbs-Bézier and penalty function in which the fluid-structure coupling effect is equivalent to added mass that is uniformly distributed over the wet area. The numerical accuracy and convergence are also examined by comparison with the open literatures. Using parametric analysis, the effects of elastic restrained edge stiffness, aspect ratio, fluid depth ratio, fiber orientation, initial pre-stress, and initial imperfection on the free vibration characteristics of composite plates are examined and discussed in detail. The results indicate that the proposed analysis method has high accuracy and the fluid-structure coupling has a significant impact on the vibration characteristics of composite plates under different fluid depth ratios.

Highlights

• (1) The proposed analysis method is capable of accurately predicting the vibration characteristic of composite plates in contact with a bounded fluid under different boundary conditions.

• (2) The dimensionless vibration frequency of composite plates in contact with the bounded fluid gradually increases as the elastic restrained edge stiffness increases.

• (3) The higher the vibration frequency, the smaller the plate aspect ratio. the effect of the added mass on the vibration characteristic of composite plates seems to be weakened by decreasing the aspect ratio.

• (4) Compared to the effect of the initial imperfection, the vibration frequency of composite plates is more sensitive to the pre-stress.

• (5) The vibration frequency of composite plates will increase at first and then decrease with increasing material orthotropy orientation under simply supported end boundary.

Introduction

Nowadays composite plates are used widely in many engineering applications including surface aircraft, underwater vehicles and ships. When the structure accelerates in a flow field, it is subjected to added mass forces from the flow field, and the added mass has a significant impact on the vibration characteristics of the structure, which cannot be ignored. Therefore, in order to prevent structural vibration fatigue and resonance damage, it is of great significance to study the prediction of underwater structural vibration considering added mass and damping. Numerous studies have been performed to investigate free and forced vibrations of composite plates in contact with a fluid.

During last decades, many scholars have applied numerical and experimental methods to study the added mass and damping of structures. Zhang et al. ¹ used the weak coupling method to analyze the fluid-structure coupling effect and investigated the convergence properties and computational efficiency of the method. Bermúdez et al. ² used finite element method to study the coupled fluid-solid vibrations of a structure in which the effect of the fluid was taken into account by means of an added mass formulation. Rajasankar et al. ³ applied the added mass to represent fluid-structure interaction effects in vibration analysis of the hulls in which a new eight-noded element had been developed to model the interface between the solid and the fluid. Liao et al. ⁴ developed a time-varying added mass method and applied the finite element method to simulate the outflow process of underwater vehicles. Kerboua et al. ⁵ developed a new method to study the vibration analysis of rectangular plates coupled with fluid by using a combination of the finite element method and Sanders’ shell theory. Cho et al.^6–7 presented a simple and effective method to investigate the natural vibration characteristics of bottom plate structures in contact with different fluid domains, taking into account fluid kinetic energy. Chen et al.^8–10 numerically studied the lift generation, underlying vorticity dynamics and coexistence of dual wing-wake interaction mechanisms during the rapid rotation of a low aspect ratio flapping wing based on the computational fluid dynamics. The finite element method (FEM) can be successfully applied to analyze the natural vibration of structures in contact with fluid. However, it needs a lot of time for model preparation and numerical calculations, and the fluid structure coupling calculation process is usually unstable, which is not applicable at the preliminary design stage. To simplify the calculation, the effects of interaction of structures in contact with fluid be taken as the added mass and damping in the analysis. Zhang et al. ¹¹ derived the added mass and damping formulations and analyzed the change of the added mass and damping during the ejection process of the projectile. Du et al. ¹² proposed the simplified hydrodynamic pressure expression based on the radiation theory in which the added mass is used for low frequency vibration and added damping is used for high frequency vibration. Song et al. ¹³ investigated the free vibration of the truncated conical shell with arbitrary boundary conditions considering elastic and inertia force constraints. Matos et al. ¹⁴ derived a closed-form solutions for the critical buckling pressures and characteristic frequencies of underwater composite cylinders accounting for added fluid mass and hydrostatic pressures. Zhou et al. ¹⁵ established the added mass and damping analogy model and its updated mathematical model to analyze the effect of turbulent flow on vibration structure. Wang et al. ¹⁶ developed a two-step theoretical approach to investigate the effect of compressive in-plane load on the equilibrium paths and vibration characteristics of a vertical symmetric laminated plate in contact with fluid. Mao et al. ¹⁷ studied the effect of skew on the vibration characteristics of the marine propellers by changing the added mass and damping matrices. Li et al. ¹⁸ developed a numerical method for computing the added mass and damping coefficients of the marine propeller by establishing the strongly coupled fluid-structure interaction models of the propellers. Zou et al. ¹⁹ presented a new method to compute the added mass and damping of the propeller based on the panel method and solved the added mass and damping of the ellipsoid and propellers. Liang et al. ²⁰ developed a simple procedure to analyze the vibration frequencies and mode shapes of cantilever plates by an empirical added mass formulation. Generally, the analytical expressions of the added mass and damping are relatively complicated that they are difficult to be utilized in practical application, especially for some complex structures. Therefore, it is meaningful work to simply the complicated expressions.

In traditional engineering practice, computer-aided design (CAD) and computer-aided engineering (CAE) are based on two different platforms. In order to achieve the integration of computer-aided design (CAD) and computer-aided engineering (CAE), Hughes et al. ²¹ introduced non-uniform rational B-spline (NURBS) basis functions into iso-parametric finite element analysis and established the method of geometric analysis (IGA). Unlike traditional finite element Lagrangian polynomial interpolation methods, iso-geometric analysis methods directly use spline basis functions in CAD as interpolation function representations in finite element analysis, which can accurately describe complex geometric configurations and eliminate geometric errors. Faroughi et al. ²² developed an isogeometric analysis method to investigate the vibration and stability analysis of in-plane laminated composite structures based on higher-order smooth NURBS basis functions. Alesadi et al. ²³ investigated the free vibration and linearized buckling analysis of laminated composite plates using the isogeometric approach and Carrera’s Unified Formulation. Lieu et al. ²⁴ investigated the Bending and free vibration analysis of in-plane bi-directional functionally graded plates with variable thickness based on isogeometric approach. Sederber et al. ²⁵ first proposed the concept of T-spline, which allows for the existence of T-nodes in the grid and achieves local refinement of the grid. Bazilevs et al. ²⁶ introduced T-splines into iso-geometric analysis, and research has shown that the T-spline basis function may have linear correlation, which will lead to severe numerical ill condition of the system stiffness matrix. Dokken et al. ²⁷ extended Nurbs and proposed LR splines, which achieve local refinement of the mesh by inserting node segments into tensor product form B-spline grids. However, LR spline basis functions may still have linear correlation and the construction of this spline basis function is relatively complex. Forsey and Bartels ²⁸ proposed a hierarchical B-spline (HB spline), which consists of a series of nested tensor product spline basis functions, enriching the local characteristics of the spline model. Vuong et al. ²⁹ proposed a method for constructing a set of linearly independent basis functions in the HB spline space and applied it to adaptive geometric analysis, achieving local refinement of the mesh for static heat conduction problems in L-shaped regions. Thai et al. ³⁰ applied the NURBS-based isogeometric approach for static, free vibration, and buckling analysis of laminated composite plates. Shojaee et al. ³¹ used the NURBS-based isogeometric finite element method for Free vibration and buckling analysis of laminated composite plates. Kapoor ³² studied the free vibration and dynamic characteristics of composite plates using Nurbs isogeometric finite element method. Zhao et al. ³³ developed a super-parametric shell element to investigate the dynamic characteristics of a rotating cylindrical shell considering parameter uncertainties. Compared to the classical Numerical methods, the proposed method utilizes NURBS basis functions which can be able to model accurately the geometries, instead of Lagrange polynomial and achieve the integration of computer-aided design (CAD) and computer-aided engineering (CAE). By employing the penalty function in the Nurbs-Bézier method, it is possible for the method to handle the complex working conditions. In addition, higher-order NURBS basis functions can provide the exact geometrical representation, high-order continuity and high accuracy.

According to the above literature survey, there is relatively little research on the effects of fluid depth ratio and boundary condition on the vibration characteristics of the composite plate. Therefore, a proposed method based on the combination of Nurbs-Bézier and penalty function is employed for the free vibration analysis of composite plates in contact with a bounded fluid in which the added mass and damping is calculated by using the LEWIS empirical formula. The numerical accuracy and applicability are compared and evaluated with those found in the open literature. Parametric studies are conducted to examine the effects of elastic restrained edge stiffness, aspect ratio, fiber orientation, initial pre-stress, and initial imperfection on the free vibration characteristics of composite plates.

A brief of NURBS function and Bézier extraction operator

Basis functions

A p-order Bézier curve can be defined as a linear combination of p + 1 Bernstein basis functions. The Bézier curve can be expressed as:

$$C\xi =\sum _{i=1}^{p+1}{\hat{P}}_i{B}_{i,p}(\xi )$$

(1)

$$B_{i,p}(\xi )=\frac{p!{(1-\xi )}^{p-(i-1)}{(1+\xi )}^{i-1}}{2^p(i-1)!(p+1-i)!}$$

(2)

Where $B_{i,p}(\xi )$ is the p order Bernstein basis function, ${\hat{P}}_i$ is the control point. The given p-order Bézier basis function can be determined by p + 1 bernstein basis functions. The Bézier basis functions on each element are the same, which have similar properties to the shape functions in finite element method.

A univariate p-order B-spline basis function is defined on an increasing node vector $\Xi =[{\xi }_1,{\xi }_2...{\xi }_{n+p+1}]$ . Based on Cox_deBoor recurrence formula, a p-order B-spline basis function $N_{i,p}(\xi )$ can be represented as a linear combination of two p − 1 order B-spline basis functions as

$$N_{i,p}(\xi )=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-1}(\xi )+\frac{{\xi }_{i+p+1}-\xi }{{\xi }_{i+p+1}-{\xi }_{i+1}}{N}_{\mathrm{i}+1,\mathrm{p}-1}(\xi )$$

(3)

$$N_{i,0}(\xi )=\{\begin{array}{cc}1\hfill & {\xi }_i\le \xi <{\xi }_{i+1}\hfill \\ \multicolumn{1}{c}{0}& \mathrm{other}\hfill \end{array}$$

(4)

Where $\xi$ is the coordinate of the variable parameter, n is the number of basis functions, ${\xi }_i$ is the ith knot.

The Nurbs basis function in the form of bivariate tensor product is defined as:

$$R_{i,j}^{p,q}(\xi ,\eta )=\frac{N_{i,p}(\xi )M_{j,q}(\eta )w_{i,j}}{\sum _{k=1}^n\sum _{t=1}^m{N}_{k,p}(\xi )M_{t,q}(\eta )w_{k,t}}$$

(5)

Where $w_{i,j}$ is the weight corresponding the B-spline basis function, $N_{i,p}(\xi )$ and $M_{j,q}(\eta )$ are p-order and q-order B-spline basis functions, respectively. When the weight $w_{i,j}$ are all equal to 1 and $\sum _{j=1}^m{N}_{j,p}(\xi )w_j=1$ , the NURBS basis function will degenerate into a B-spline basis function.

Therefore, the coordinate vector of any point on the Nurbs surface can be expressed as:

$$z(\xi ,\eta )=\sum _{i=1}^n\sum _{j=1}^m{R}_{i,j}^{p,q}(\xi ,\eta )P_{i,j}$$

(6)

Where $P_{i,j}$ is the control point, $\Xi =[{\xi }_1,{\xi }_2...{\xi }_{n+p+1}]$ and $\mathrm{\Omega }=[{\eta }_1,{\eta }_2...{\eta }_{m+q+1}]$ are the node vectors in the $\xi$ and $\eta$ directions, respectively (Figure 1).

Figure 1.

Bivariate basis function.

Bézier extraction operator

In order to integrate geometric analysis programs with existing finite element programs, Bézier extraction approach is widely applied to various spline basis functions. The core idea of Bézier extraction is to use the Bernstein basis functions to represent globally B-spline basis functions based on node insertion algorithm. Figure 2 shows the Bézier extraction process of univariate B-splines.

Figure 2.

Bezier extraction of univariate B-spline: (a) Bézier curve, (b) Nurbs curve.

If a new node is inserted in the parameter field $[{\xi }_i,{\xi }_{i+1}]$ , the node vector $\xi \prime \in [{\xi }_i,{\xi }_{i+1}]$ will become $\Xi \prime =[{\xi }_1...{\xi }_i\xi \prime {\xi }_{i+1}...{\xi }_{n+p+1}]$ , p < i < n + 1, and a new basis function and control point will be added. The coordinate of the new control point $P\prime$ can be calculated by the initial control point P as:

$$\begin{array}{c}\multicolumn{1}{c}{P_j\prime =\{\begin{array}{c}{P}_{1}j=1\hfill \\ \multicolumn{1}{c}{{\alpha }_j{P}_j+(1-{\alpha }_j)P_{j-1}1<j<m}\\ \multicolumn{1}{c}{P_{n}j=m}\end{array}}\end{array}$$

(7)

$$\begin{array}{c}\multicolumn{1}{c}{{\alpha }_j=\{\begin{array}{c}11\le j\le i-p\hfill \\ \multicolumn{1}{c}{\frac{\hat{\xi }-{\xi }_j}{{\xi }_{j+p}-{\xi }_j}i-p+1\le j\le i}\\ \multicolumn{1}{c}{0j\ge i+1}\end{array}}\end{array}$$

(8)

Based on the node insertion algorithm in equations (8) and (9), a B-spline basis function can be transformed into a Bernstein basis function. The control point can be represented as a linear combination of control points P in the B-spline space.

$$\hat{P}=C^{T}P$$

(9)

In order to keepno geometric or parametric change during node insertion, it is necessary to

$$x(\xi )=P^{T}N(\xi )={\hat{P}}^{T}B(\xi )=(C^{T}P{)}^{T}B(\xi )=P^T\mathrm{CB}(\xi )$$

(10)

Therefore, the relationship between the B-spline basis function and the Bernstein basis function is

$$N(\xi )=\mathrm{CB}(\xi )$$

(11)

Where the linear operator C is called the global Bézier extraction operator. The element Bézier extraction operator $C^e$ can be obtained from the global Bézier extraction operator $C$ . Therefore, the element B-spline basis function can be represented as

$$N(\xi ){|}_{{\mathrm{\Omega }}^e}=N^e(\overline{\xi })=C^{e}B(\overline{\xi })$$

(12)

Where $\overline{\xi }=[-1,1]$ is the parameter space of the Bézier element. The multivariate Bézier extraction operator with dimension $d_p$ can be obtained by the Kronecker product of the univariate extraction operators,

$$C^e=C_1^e\otimes C_2^e...\otimes C_{d_p}^e$$

(13)

The Kronecker product of two matrices $A\otimes B$ is

$$A\otimes B=\left[\begin{array}{ccc}{A}_{11}B\hfill & A_{12}B\hfill & ...\hfill \\ \multicolumn{1}{c}{A_{21}B}& A_{22}B\hfill & ...\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill \end{array}\right]$$

(14)

The Nurbs basis function can be expressed as:

$$R^e(\xi )=W^e{C}^e\frac{B(\overline{\xi })}{{w^e}^T{C}^{e}B(\overline{\xi })}$$

(15)

Where $w^e$ is the weight associated with the B-spline basis function, $\hat{P}$ is the control point associated with the B-spline basis function, $W^e$ is the diagonal matrix of the element weight. Features of the element are indicated by a subscript e.

Theoretical formulation

Calculation of added mass of plates in contact with the bounded fluid

An analytical method for the vibration characteristics of plates with initial pre-stress is developed in this paper. The initial pre-stress distribution within the plates under non-uniform mechanical edge loading is developed by using buckling analysis of plates in which the amplitude of initial pre-stress is represented by the critical buckling load and load scale factor.

Consider a composite rectangular plate with length a, width b, total thickness h as shown in Figure 3 in which the plate under in-plane combined loads $N_x$ and $N_{\mathrm{xy}}$ . The longitudinal and transverse elastic restrained edge stiffnesses are defined as $k_1$ and $k_2$ , respectively.

Figure 3.

Geometry feature of composite plates with elastically restrained edges.

Figure 4 shows a composite rectangular plate in contact with the bounded fluid. The distance between the composite plate and the side of water tank are $L_x$ and $L_y$ , respectively. The depth of the fluid is $L_z$ , the immersion depth of the plate is $b_1$ , the water density is ${\rho }_{\mathrm{w}}$ .

Figure 4.

Geometry and coordinates of the composite plate in contact with the bounded fluid: (a) front view, (b) side view.

When a seaplane glides on the water surface, the inherent characteristics and vibration response of the structure are affected by the water load, so it is necessary to consider the effect of added mass and damping on the natural frequency of the seaplane. This paper takes a composite plate as the research object and uses the LEWIS empirical formula to calculate the added mass of the plate in contact with the bounded fluid in which the added mass is uniformly distributed and incorporated into the vibration equation. ³⁴ Assuming the composite plate is not completely immersed in water, and the fluid depth ratio has a linear relationship with the added mass, thus the added mass is defined as:

$$m_{\mathrm{ad}}=\frac{1}{2}{C}_v{\rho }_w\pi a^2$$

(16)

Where $C_v$ is correction factor related to fluid depth ratio, the water density is ${\rho }_{\mathrm{w}}$ , the length of the plate is $a$ , the width of the plate is $b$ , the thickness of the plate is $h$ . In this paper, the correction factor is defined as $C_v=1.0(b_1/b)$ .

Geometric mapping

In order to simply the calculation, the quadrilateral physical domain needs to be mapped to the square computational domain by coordinate transformation. To transform any irregular quadrilateral physical domain $x-y$ in Figure 5(a) into square computational domain $\xi -\eta$ in Figure 5(b) in natural coordinate system, the linear coordinate transformation formula can be defined as follows:

$$\begin{array}{c}x=\sum _{i=1}^8{N}_i(\xi ,\eta )x_i\hfill \\ y=\sum _{i=1}^8{N}_i(\xi ,\eta )y_i\hfill \end{array}$$

(17)

Where $x_i$ and $y_i$ are the coordinates of node i in the quadrilateral physical domain, $N_i(\xi ,\eta )$ is the shape function associated with node i. The Lagrangian interpolation function is defined as:

$$N_i(\xi ,\eta )=\frac{1}{4}(1+\xi {\xi }_i)(1+\eta {\eta }_i)(\xi {\xi }_i+\eta {\eta }_i-1)i=1,3,5,7$$

$$N_i(\xi ,\eta )=\frac{1}{2}(1-{\xi }^2)(1+\eta {\eta }_i)i=2,6$$

(18)

$$N_i(\xi ,\eta )=\frac{1}{2}(1+{\xi }_i\xi )(1-{\eta }^2)i=4,8$$

Figure 5.

The transformation of geometry and coordinate system: (a) quadrilateral physical domain, (b) square computational domain.

Based on the equation (17), the first-order and second-order partial derivative matrices of the functions in two coordinate systems are written as follows:

$$\left[\begin{array}{c}({)}_{,x}\hfill \\ \multicolumn{1}{c}{({)}_{,y}}\end{array}\right]=[J_0{]}^{-1}\left[\begin{array}{c}({)}_{,\xi }\hfill \\ \multicolumn{1}{c}{({)}_{,\eta }}\end{array}\right]$$

(19)

$$\left[\begin{array}{c}({)}_{,\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{({)}_{,\mathrm{yy}}}\\ \multicolumn{1}{c}{({)}_{,\mathrm{xy}}}\end{array}\right]=[J_1{]}^{-1}\left[\begin{array}{c}({)}_{,\xi \xi }\hfill \\ \multicolumn{1}{c}{({)}_{,\xi \eta }}\\ \multicolumn{1}{c}{({)}_{,\eta \eta }}\end{array}\right]-[J_1{]}^{-1}[J_2][J_0{]}^{-1}\left[\begin{array}{c}({)}_{,\xi }\hfill \\ \multicolumn{1}{c}{({)}_{,\eta }}\end{array}\right]$$

(20)

Where $J_0$ , $J_1$ , $J_2$ are Jacobian transformation matrices.

$$\begin{array}{c}[J_{,0}]=\left[\begin{array}{cc}{x}_{,\xi }\hfill & y_{,\xi }\hfill \\ \multicolumn{1}{c}{x_{,\eta }}& y_{,\eta }\hfill \end{array}\right],[J_2]=\left[\begin{array}{cc}{x}_{,\xi \xi }\hfill & y_{,\xi \xi }\hfill \\ \multicolumn{1}{c}{x_{,\eta \eta }}& y_{,\eta \eta }\hfill \\ \multicolumn{1}{c}{x_{,\xi \eta }}& y_{,\xi \eta }\hfill \end{array}\right],\hfill \\ [J_1]=\left[\begin{array}{ccc}{x}_{,\xi }^2\hfill & y_{,\xi }^2\hfill & 2x_{,\xi }{y}_{,\xi }\hfill \\ \multicolumn{1}{c}{x_{,\eta }^2}& y_{,\eta }^2\hfill & 2x_{,\eta }{y}_{,\eta }\hfill \\ \multicolumn{1}{c}{x_{,\xi }{x}_{,\eta }}& y_{,\xi }{y}_{,\eta }\hfill & x_{,\xi }{x}_{,\eta }+x_{,\eta }{y}_{,\xi }\hfill \end{array}\right]\hfill \end{array}$$

(21)

The governing equations of composite plates

In this paper, the first-order shear deformation theory (FSDT) is adopted to derive the theoretical model for the vibration analysis of the composite plate, taking into account geometric nonlinearity and initial imperfection. Based on the classical nonlinear plate theory, the corresponding strain components can be written as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}{\epsilon }_x\hfill \\ \multicolumn{1}{c}{{\epsilon }_y}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xy}}}\end{array}\right\}=\left\{\begin{array}{c}{\epsilon }_x^m\hfill \\ \multicolumn{1}{c}{{\epsilon }_y^m}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xy}}^m}\end{array}\right\}+z\left\{\begin{array}{c}{\chi }_x\hfill \\ \multicolumn{1}{c}{{\chi }_y}\\ \multicolumn{1}{c}{{\chi }_{\mathrm{xy}}}\end{array}\right\}}\end{array}$$

(22)

Where

$$\epsilon ={\epsilon }^m+z\chi ,{\epsilon }^m=\{{\epsilon }_x^m,{\epsilon }_y^m,{\gamma }_{\mathrm{xy}}^m\},\chi =\{{\chi }_x,{\chi }_y,{\chi }_{\mathrm{xy}}\}$$

$$\begin{array}{c}\multicolumn{1}{c}{{\epsilon }^m=\left\{\begin{array}{c}{\epsilon }_x^m\hfill \\ \multicolumn{1}{c}{{\epsilon }_y^m}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xy}}^m}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial w^{*}}{\partial x}\hfill \\ \frac{\partial v}{\partial y}+\frac{1}{2}{\left(\frac{\partial w}{\partial y}\right)}^2+\frac{\partial w}{\partial y}\frac{\partial w^{*}}{\partial y}\hfill \\ \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}+\frac{\partial w}{\partial x}\frac{\partial w}{\partial y}+\frac{\partial w^{*}}{\partial x}\frac{\partial w}{\partial y}+\frac{\partial w}{\partial x}\frac{\partial w^{*}}{\partial y}\hfill \end{array}\right\}}\\ ={\epsilon }_L+{\epsilon }_0+{\epsilon }_{\mathrm{NL}},\chi =\left\{\begin{array}{c}{\chi }_x\hfill \\ \multicolumn{1}{c}{{\chi }_y}\\ \multicolumn{1}{c}{{\chi }_{\mathrm{xy}}}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial {\phi }_x}{\partial x}\hfill \\ \frac{\partial {\phi }_y}{\partial y}\hfill \\ \frac{\partial {\phi }_x}{\partial y}+\frac{\partial {\phi }_y}{\partial x}\hfill \end{array}\right\}\hfill \end{array}$$

(23)

$${\gamma }_{\mathrm{xz}}=\frac{\partial w}{\partial x}+\frac{\partial w^{*}}{\partial x}+{\phi }_x,{\gamma }_{\mathrm{yz}}=\frac{\partial w}{\partial y}+\frac{\partial w^{*}}{\partial y}+{\phi }_y$$

Where $u$ , $v$ , and $w$ are the displacements at each point of the middle surface along the x, y, z coordinate axes, ${\phi }_x$ and ${\phi }_y$ are the rotation angles of the normal of the middle surface about the x, y coordinate axes, ${\epsilon }_L$ , ${\epsilon }_0$ , and ${\epsilon }_{\mathrm{NL}}$ are the in-plane linear strain, initial strain and nonlinear strain, respectively. $w^{*}$ is the initial deflection along the z coordinate axis. Features of the middle surface are indicated by a subscript 0.

According to linear elastic theory, the stress-strain relationship of the plate can be written as:

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}{\sigma }_x\hfill \\ \multicolumn{1}{c}{{\sigma }_y}\\ \multicolumn{1}{c}{{\tau }_{\mathrm{xy}}}\end{array}\right\}=\left[\begin{array}{ccc}{Q}_{11}\hfill & Q_{12}\hfill & Q_{16}\hfill \\ \multicolumn{1}{c}{Q_{11}}& Q_{12}\hfill & Q_{26}\hfill \\ \multicolumn{1}{c}{Q_{61}}& Q_{62}\hfill & Q_{66}\hfill \end{array}\right]\left\{\begin{array}{c}{\epsilon }_x\hfill \\ \multicolumn{1}{c}{{\epsilon }_y}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xy}}}\end{array}\right\},\left\{\begin{array}{c}{\tau }_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{{\tau }_{\mathrm{yz}}}\end{array}\right\}}\\ =k_s\left[\begin{array}{cc}{Q}_{44}\hfill & Q_{45}\hfill \\ \multicolumn{1}{c}{Q_{54}}& Q_{55}\hfill \end{array}\right]\left\{\begin{array}{c}{\gamma }_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{{\gamma }_{\mathrm{yz}}}\end{array}\right\}\hfill \end{array}$$

(24)

The force and moment resultants of the plate can be obtained by integrating the stress along the thickness as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\left(\left(\begin{array}{c}{N}_x\hfill \\ \multicolumn{1}{c}{N_y}\\ \multicolumn{1}{c}{N_{\mathrm{xy}}}\end{array}\right)\left(\begin{array}{c}{M}_x\hfill \\ \multicolumn{1}{c}{M_y}\\ \multicolumn{1}{c}{M_{\mathrm{xy}}}\end{array}\right)\right)=\underset{-h/2}{\overset{h/2}{\int }}\left(\begin{array}{c}{\sigma }_x\hfill \\ \multicolumn{1}{c}{{\sigma }_y}\\ \multicolumn{1}{c}{{\tau }_{\mathrm{xy}}}\end{array}\right)(1,z)\mathrm{dz},\{Q_x,Q_y\}}\\ =\underset{-h/2}{\overset{h/2}{\int }}\{{\tau }_{\mathrm{xz}},{\tau }_{\mathrm{yz}}\}\mathrm{dz},\hfill \end{array}$$

(25)

The force and moment resultants can then be obtained from integration of equations (20), (22)–(24):

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}\left\{N\right\}\hfill \\ \multicolumn{1}{c}{\left\{M\right\}}\end{array}\right\}=\left[\begin{array}{cc}\left[A\right]\hfill & \left[B\right]\hfill \\ \multicolumn{1}{c}{\left[B\right]}& \left[D\right]\hfill \end{array}\right]\left\{\begin{array}{c}\left\{{\epsilon }^m\right\}\hfill \\ \multicolumn{1}{c}{\left\{\kappa \right\}}\end{array}\right\}=[\hat{D}]\left\{\begin{array}{c}\left\{{\epsilon }^m\right\}\hfill \\ \multicolumn{1}{c}{\left\{\kappa \right\}}\end{array}\right\},\left\{\begin{array}{c}{Q}_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{Q_{\mathrm{yz}}}\end{array}\right\}}\\ =\left[\begin{array}{cc}{A}_{44}\hfill & A_{45}\hfill \\ \multicolumn{1}{c}{A_{54}}& A_{55}\hfill \end{array}\right]\left\{\begin{array}{c}{\gamma }_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{{\gamma }_{\mathrm{yz}}}\end{array}\right\}=D_s\left\{\begin{array}{c}{\gamma }_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{{\gamma }_{\mathrm{yz}}}\end{array}\right\}\hfill \end{array}$$

(26)

Where $[A]$ is the extensional stiffness matrix, $[B]$ is the coupling stiffness matrix, and $[D]$ is the bending stiffness matrix, $[D_s]$ is the shear stiffness matrix, $\hat{D}$ is the material stiffness matrix.

$$\left\{N\right\}=[N_x,N_y,N_{\mathrm{xy}}]\left\{M\right\}=\{M_x,M_y,M_{\mathrm{xy}}\}\left\{Q\right\}=\{Q_{\mathrm{xz}},Q_{\mathrm{yz}}\}$$

(27)

$$\left[A\right]=\left[\begin{array}{ccc}{A}_{11}\hfill & A_{12}\hfill & A_{16}\hfill \\ \multicolumn{1}{c}{A_{21}}& A_{22}\hfill & A_{26}\hfill \\ \multicolumn{1}{c}{A_{16}}& A_{26}\hfill & A_{66}\hfill \end{array}\right],\left[B\right]=\left[\begin{array}{ccc}{B}_{11}\hfill & B_{12}\hfill & B_{16}\hfill \\ \multicolumn{1}{c}{B_{21}}& B_{22}\hfill & B_{26}\hfill \\ \multicolumn{1}{c}{B_{16}}& B_{26}\hfill & B_{66}\hfill \end{array}\right]$$

(28)

$$\left[D\right]=\left[\begin{array}{ccc}{D}_{11}\hfill & D_{12}\hfill & D_{16}\hfill \\ \multicolumn{1}{c}{D_{21}}& D_{22}\hfill & D_{26}\hfill \\ \multicolumn{1}{c}{D_{16}}& D_{26}\hfill & D_{66}\hfill \end{array}\right],D_s=\left[\begin{array}{cc}{A}_{44}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& A_{55}\hfill \end{array}\right]$$

Where $N_x$ , $N_y$ , and $N_{\mathrm{xy}}$ are the force resultants, $M_x$ , $M_y$ , and $M_{\mathrm{xy}}$ are the moment resultants, $Q_{\mathrm{xz}}$ , and $Q_{\mathrm{yz}}$ are the shear resultants, ${\epsilon }_x$ , ${\epsilon }_y$ , ${\gamma }_{\mathrm{xy}}$ , ${\gamma }_{\mathrm{xz}}$ , and ${\gamma }_{\mathrm{yz}}$ are the normal and shear strains, ${\chi }_x$ , ${\chi }_y$ , and ${\chi }_{\mathrm{xy}}$ are the curvature parameters, $h$ is the thickness of the plate (Appendix 1).

Equations of motion

For the vibration characteristics of composite plates, the equilibrium equation of the composite plate with elastically restrained edges is established by applying minimal potential energy method. The equilibrium equation for the composite plate can be written as follows:

$$\delta (T-(U+U_s-W))=0$$

(29)

Where $T$ $U$ , $U_s$ , and $W$ indicate the kinetic energy, the potential strain energy, the spring potential energy andthe external work. The kinetic energy is given as follows:

$$T=\underset{A}{\int }{\int }_{-h/2}^{h/2}\left({\left(\frac{\partial u}{\partial t}\right)}^2+{\left(\frac{\partial v}{\partial t}\right)}^2+{\left(\frac{\partial w}{\partial t}\right)}^2\right)dAdz$$

(30)

The potential strain energy stored in the plate is given as follows:

$$U=\underset{A}{\int }{\int }_{-h/2}^{h/2}({\sigma }_x{\epsilon }_x+{\sigma }_y{\epsilon }_y+{\tau }_{\mathrm{xy}}{\gamma }_{\mathrm{xy}}+{\tau }_{\mathrm{xz}}{\gamma }_{\mathrm{xz}}+{\tau }_{\mathrm{yz}}{\gamma }_{\mathrm{yz}})\mathrm{dAdz}$$

(31)

For the elastically restrained plate, the potential energy stored in the boundary springs is given by the following:

$$\begin{array}{c}\multicolumn{1}{c}{U_s=\frac{1}{2}{\int }_0^b[(k_2{u}^2+k_2{v}^2+k_2{w}^2+k_2{\phi }_x+k_2{\phi }_y){|}_{x=0}}\\ +(k_2{u}^2+k_2{v}^2+k_2{w}^2+k_2{\phi }_x+k_2{\phi }_y){|}_{x=a}]\mathrm{dx}\hfill \\ \frac{1}{2}{\int }_0^a[(k_1{u}^2+k_1{v}^2+k_1{w}^2+k_1{\phi }_x+k_1{\phi }_y){|}_{y=0}\hfill \\ +(k_1{u}^2+k_1{v}^2+k_1{w}^2+k_1{\phi }_x+k_1{\phi }_y){|}_{y=b}]\mathrm{dy}\hfill \end{array}$$

(32)

The virtual work done by the edge loading is given as follows:

$$W={\int }_{\mathrm{\Gamma }}u·\mathrm{fd}\mathrm{\Gamma }$$

(33)

Where A is the integration area, t is the action time, $f$ is the edge load on the natural boundary, N is the number of B-spline basis functions, $k_1$ and $k_2$ are the longitudinal and transverse elastic restrained edge stiffnesses, $U=\{u,v,w,{\phi }_x,{\phi }_y{\}}^T$ represents the rotation and displacement vectors. $I_0$ , $I$ , and $I_2$ are indicated as follows

$$\{I_0,I_1,I_2\}={\int }_{-h/2}^{h/2}\rho (x,y)\{1,z,z^2\}\mathrm{dz}$$

(34)

Substituting the stress and strain equations into equations (29), the final discrete equilibrium equation by using the Nurbs-Bézier method is obtained as

$$(M+m_{\mathrm{ad}})U^{″}+C_{\mathrm{ad}}U\prime +\mathrm{KU}=f(t)$$

(35)

Where $K$ represents the total stiffness matrix of the composite plate, $K_g$ represents the geometric stiffness matrix of the composite plate, $M$ represents the mass matrix of the composite plate, $m_{\mathrm{ad}}$ represents the added mass matrixof the composite plate, $C_{\mathrm{ad}}$ represents the damping stiffness matrix of the composite plate, $f(t)$ represents the edge load of the composite plate, $U=\{u,v,w,{\phi }_x,{\phi }_y{\}}^T$ represents the rotationand displacement vectors of the control point. $U\prime$ represents the first order derivative of the rotation and displacement vectors. $U^{″}$ represents the second order derivative of the rotationand displacement vectors.

In order to transform the system equation into a general eigenvalue problem, the auxiliary equation is introduced:

$$(M-M)U\prime =0$$

(36)

Therefore, general eigenvalue equation can be written as follows:

$$\begin{array}{c}{\left[\begin{array}{cc}{C}_{\mathrm{ad}}\hfill & M+m_{\mathrm{ad}}\hfill \\ \multicolumn{1}{c}{M+m_{\mathrm{ad}}}& 0\hfill \end{array}\right]}_{2n\times 2n}{\left\{\begin{array}{c}U\prime \hfill \\ \multicolumn{1}{c}{U^{″}}\end{array}\right\}}_{2n\times 1}\hfill \\ +{\left[\begin{array}{cc}{K}_L+K_{\mathrm{NL}}+\lambda N_{\mathrm{cr}}{K}_g\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& M+m_{\mathrm{ad}}\hfill \end{array}\right]}_{2n\times 2n}{\left\{\begin{array}{c}U\hfill \\ \multicolumn{1}{c}{U\prime }\end{array}\right\}}_{2n\times 1}={\left\{\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\right\}}_{2n\times 1}\hfill \end{array}$$

(37)

Where $N_{\mathrm{cr}}$ represents the critical buckling loadof the composite plate, $\lambda$ represents the load scale factor. The total stiffness matrix is defined as $K=K_L+K_{\mathrm{NL}}$ . The damping matrix is defined as:

$$C_{\mathrm{ad}}=\alpha K+\beta M$$

(38)

Numerical results and discussion

Convergence and comparison studies

The convergence is first studied to prove the accuracy and accuracy of the analysis results by changing the number of B-spline basis functions in the trial function. A square plate with elastically restrained edges is considered. The geometric parameters and material properties of the aluminium plate are as follows: modulus of elasticity $E=70\mathrm{GPa}$ , Poisson’s ratio $v=0.33$ , aspect ratio $a/b=1$ and thickness-to-side ratio $h/a=0.01$ . The dimensionless transverse and longitudinal elastic restrained edge stiffness coefficients ( $K_{y0}^r=K_{\mathrm{yb}}^r=K_1$ , $K_{x0}^r=K_{\mathrm{xa}}^r=K_2$ ) are defined as follows:

$$\begin{array}{c}{K}_{x0}^r=\frac{k_{2}b}{\sqrt{D_{11}{D}_{22}}},K_{\mathrm{xa}}^r=\frac{k_{2}b}{\sqrt{D_{11}{D}_{22}}},K_{y0}^r=\frac{k_{1}a}{\sqrt{D_{11}{D}_{22}}},\hfill \\ K_{\mathrm{yb}}^r=\frac{k_{1}a}{\sqrt{D_{11}{D}_{22}}}\hfill \end{array}$$

(39)

Where, $D_{11}$ , $D_{22}$ are the bending stiffnesses, $\omega$ is the in-plane vibration frequency. It can be seen that the analysis results listed in Tables 1 and 2 satisfy the accuracy requirements when the number of basis function is $N_x=N_y=17$ . In addition, it is further known that the calculation results of the elastic restrained edge stiffness coefficients $K_1=K_2=0$ and $K_1=K_2={10}^8$ are in well agreement with results of simply supported and clamped end conditions by Liew and Ferreira.^35,36 A dimensionless frequency parameter of the plate is given as follows:

$$\varpi =\omega a^2\sqrt{\rho /G}$$

(40)

Where, $\omega$ is the frequency, $\rho$ is the mass density, G is the shear modulus.

Table 1.

Dimensionless natural frequencies ( $\omega a^2\sqrt{\rho /G}$ ) of a SSSS square Mindlin/Reissner plate.

Mode no	13 × 13	15 × 15	17 × 17	Ref. 35	Ref. 36
1	0.0967	0.0963	0.0963	0.0961	0.0963
2	0.2408	0.2404	0.2404	0.2419	0.2405
3	0.2409	0.2405	0.2404	0.2419	0.2407
4	0.3849	0.3844	0.3843	0.3860	0.3848

Table 2.

Dimensionless natural frequencies ( $\omega a^2\sqrt{\rho /G}$ ) of a CCCC square Mindlin/Reissner plate.

Mode no	13 × 13	15 × 15	17 × 17	Ref. 35	Ref. 36
1	0.1754	0.1752	0.1752	0.1743	0.1753
2	0.3571	0.3566	0.3566	0.3576	0.3572
3	0.3571	0.3566	0.3566	0.3576	0.3574
4	0.5253	0.5246	0.5246	0.5240	0.5273

In this section, the accuracy and convergence of the analysis method for the free vibration analysis of the composite plate in contact with the bounded fluid is demonstrated by comparative study. A plate with dimensions and material properties a = 10 m, b = 1 m, h = 0.15 m, $\rho$  = 2400 kg/m³, E = 25 GPa and v = 0.15. It can be seen that the calculated results listed in Tables 3–5 satisfy the accuracy and convergence requirements for all fluid depth ratios (b₁/b = 0, 0.2, 0.4, 0.6, 0.8, and 1) when the number of polynomial terms is $N_x=N_y=17$ . The fundamental dimensionless vibration frequency of the plate in contact with fluid is given as follows:

$$\varpi =\omega a^2\sqrt{\rho h/D_{11}}$$

(40)

Table 3.

Dimensionless natural frequencies ( $\omega a^2\sqrt{\rho h/D_{11}}$ ) for a SSSS square plate in contact with the bounded fluid.

b 1 /b	13 × 13	15 × 15	17 × 17	Ref. 37	Ref. 38
0	3.139	3.136	3.136	3.169	3.142
0.2	2.989	2.985	2.985	3.064	3.013
0.4	1.938	1.935	1.935	2.196	2.075
0.6	1.327	1.322	1.322	1.496	1.356
0.8	1.031	1.027	1.027	1.173	1.017
1.0	0.914	0.909	0.909	1.036	0.857

Table 4.

Dimensionless natural frequencies ( $\omega a^2\sqrt{\rho h/D_{11}}$ ) for a CCCC square plate in contact with the bounded fluid.

b 1 /b	13 × 13	15 × 15	17 × 17	Ref. 37	Ref. 39
0	5.721	5.712	5.712	5.730	5.715
0.2	5.587	5.580	5.580	5.674	-
0.4	3.824	3.817	3.817	4.202	4.043
0.6	2.429	2.423	2.422	2.747	-
0.8	1.853	1.847	1.847	2.157	1.940
1.0	1.662	1.655	1.655	1.948	1.673

Table 5.

Dimensionless natural frequencies ( $\omega a^2\sqrt{\rho h/D_{11}}$ ) for a SSSS square plate in contact with fluid.

b 1 /b	Present	Ref. 37	Ref. 38	Ref. 40
0	3.136	3.169	3.142	3.139
0.2	2.985	3.064	3.013	3.020
0.4	1.935	2.196	2.075	2.191
0.6	1.322	1.496	1.356	1.414
0.8	1.027	1.173	1.017	1.031
1.0	0.909	1.036	0.857	0.860

To verify the accuracy and effectiveness of the proposed numerical method, A test result is cited in the paper. The plate is assumed to be uniform, homogenous, and isotropic. The material is mild steel, density $\rho =7850\mathrm{kg}/m^3$ , Elastic modulus $E=210\mathrm{GPa}$ , Poisson’s ratio $\nu =0.3$ , side length $L=0.6096m$ , thickness $h=0.01a$ . It is also observed that the mode sequence changes according to the fluid depth. Table 6 shows that there is a good agreement between the present results and experiment of Verma et al., ⁴¹ in which the maximum value of relative error compared with the experimental results is less than 10.0% for all the fluid depth ratios. The fundamental dimensionless vibration frequency of the plate in contact with fluid is given as follows:

$$\varpi =\omega a^2\sqrt{\rho /D_{11}}$$

(41)

Table 6.

Dimensionless natural frequencies ( $\omega a^2\sqrt{\rho /D_{11}}$ ) for a CFFF square plate in contact with the bounded fluid.

b 1 /b	Semi-analytical 41	ANSYS 41	Experiment 41	Present
0	3.48	3.47	3.51	3.497
0.25	2.562	2.587	2.622	2.715
0.5	2.064	2.042	1.987	2.139
0.75	1.877	1.886	1.715	1.7045

Parametric studies

In this section, numerical results for free vibrations of the composite plate in contact with the bounded fluid are obtained by the proposed analytical method. In the present study, two different boundary conditions are investigated: simply supported and clamped end conditions according to the definition of elastic restrained edge stiffness. The orthotropic material has the following dimensions and mechanical properties: a = 1 m, b = 1 m, h = 0.01 m, E₁ = 150 GPa, E₂ = E₃ = 9 GPa, G₁₂ = G₁₃ = 7.1 GPa, G₂₃ = 2.5 GPa, v₁₂ = v₂₃ = 0.3, $\rho$  = 1600 kg/m³. The fundamental dimensionless vibration frequency of the plate in contact with fluid is given as follows: $\varpi =\omega a^2\sqrt{\rho h/D_{11}}$ .

Effect of elastic restrained edge stiffness on vibration characteristics of composite plates

In Figure 6, the effect of elastic restrained edge stiffness on the vibration characteristics of composite plates in contact with the bounded fluid is examined. It can be observed from Figure 6 that the dimensionless vibration frequency of composite plates in contact with the bounded fluid gradually increases as the elastic restrained edge stiffness increases. On the contrary, the dimensionless vibration frequency of composite plates decreases as the fluid depth ratio increases. Significantly, the vibration frequency reduction rate of composite plates decreases as the fluid depth ratio increases and elastic restrained edge stiffness decreases. This illustrates that the vibration characteristics of composite plates is more sensitive to the fluid depth ratio. This is due to the improvement of bending stiffness of composite plates by raising the elastic restrained edge stiffness.

Figure 6.

Dimensionless vibration frequency $\varpi$ versus elastic restrained edge stiffness parameter (K) for composite plates: (a) effect of elastic restrained edge stiffness, (b) effect of fluid depth ratio.

Effect of the plate aspect ratio on vibration characteristics of composite plates

The effect of the aspect ratio $\gamma =a/b$ on the vibration characteristics of composite plates under different fluid depth ratios is shown in Figure 7. The dimensionless vibration frequency of composite plates is plotted against the aspect ratio for simply supported and clamped end conditions. The change of the dimensionless vibration frequency of composite plates is more pronounced while the aspect ratio $\gamma$ is in the range of 0.5∼1, and there is not obvious while the aspect ratio is over 1. The smaller the plate aspect ratio $\gamma$ , The higher the vibration frequency $\varpi$ . This can be explained by the effect of the added mass on the vibration characteristics of composite plates seems to be weakened by decreasing the aspect ratio. Together with the observations in Figure 7, it is concluded that the vibration frequency of composite plates under different boundary conditions can be reduced by dispersing added mass into the general mass matrix.

Figure 7.

Dimensionless vibration frequency $\varpi$ versus aspect ratio γ for composite plates under different fluid depth ratios (b₁/b): (a) simply supported end condition, (b) clamped supported end condition.

Effect of compressive pre-stress on vibration characteristics of composite plates

The effect of the compressive pre-stress on the vibration characteristics of composite plates under different fluid depth ratios is investigated in Figure 8. It is well known that the compressive pre-stress significantly affects the vibration frequency and mode shape of composite plates in contact with the bounded fluid. It can be observed that with the raise of the compressive load scale factor ${\lambda }_c$ varies from 0 to 1, the vibration frequencyof composite plates decreases first and then increases when the fluid depth ratio is relatively large. When the fluid depth ratio is relatively small, the vibration frequency of composite plates gradually decreases. This is due to the shape change of the first-order vibration mode when the pre-stress approaches the critical buckling load. This means that we can estimate the lowest vibration frequency of composite plates under different fluid depth ratios by vibration characteristic analysis. The trend of vibration frequency variation is basically consistent for different fluid depth ratios. We could also find that the greater the compressive load scale factor, the higher vibration frequency change rate. This is due to the fact that the mode shape of composite plates will change while the compressive load scale factor is over 0.9.

Figure 8.

Dimensionless vibration frequency $\varpi$ versus compressive load scale factor ${\lambda }_c$ for composite plates under different fluid depth ratios (b₁/b): (a) simply supported end condition, (b) clamped supported end condition.

Effect of tensile pre-stresson vibration characteristics of composite plates

In Figure 9, the numerical results are shown for composite plates partially in contact with fluid where the tensile load scale factor ${\lambda }_t$ varies from 0 to 1. Figure 9 shows that the vibration frequency of composite plates monotonically increases as the tensile load scale factor increases under different fluid depth ratios. It should be noted that the tensile load scale factor significantly affects the vibration frequency and mode shape of composite plates under different the fluid depth ratios (b₁/b). The results show that with the raise of the tensile load scale factor, the vibration frequency of composite plates continues to increase. The variation trend of vibration frequency variation is basically consistentfor different fluid depth ratios. We could also find that the greater the tensile load scale factor, the higher vibration frequency variation rate. It is worth emphasizing that the mode shape of composite plates for different pre-stress levels is almost unchanged.

Figure 9.

Dimensionless vibration frequency $\varpi$ versus tensile load scale factor ${\lambda }_t$ for composite plates under different fluid depth ratios (b₁/b): (a) simply supported end condition, (b) clamped supported end condition.

Effect of shear pre-stresson vibration characteristics of composite plates

The typical vibration frequency of composite plates in contact with fluid under different shear load scale factors is illustrated in Figure 10. It should be noted that the shear pre-stress significantly affects the wet vibration frequency and mode shape of composite plates under different fluid depth ratios (b₁/b). The results show that with the raise of the shear load scale factor ${\lambda }_s$ , the vibration frequency of composite plates gradually continues to decrease and the variation trend of vibration frequency is basically consistent for different fluid depth ratios. The results show that the greater the shear load scale factor, the higher vibration frequency variation rate, particularly the shear load scale factor ${\lambda }_s$ varies from 0.6 to 1.0. In addition, for different shear pre-stress levels, the mode shape of composite plates is almost unchanged.

Figure 10.

Dimensionless vibration frequency $\varpi$ versus shear load scale factor ${\lambda }_s$ for composite plates under different fluid depth ratios (b₁/b): (a) simply supported end condition, (b) clamped supported end condition.

Effect of initial imperfection on vibration characteristics of composite plates

The effect of initial imperfection amplitude on the vibration characteristics of composite plates under different fluid depth ratios is investigated in Figure 11. It can be observed that the vibration frequencies of composite plates in contact with a bounded fluid are always greater than its corresponding dry vibration frequencies owing to the fluid-structure coupling effect. As the initial imperfection amplitude increases, the normalized vibration frequencies ofcomposite plates gradually increase. It is due to increasing the initial imperfection amplitude of composite plates leads to increasing the geometric nonlinear stiffness.

Figure 11.

Dimensionless vibration frequency $\varpi$ versus initial imperfection amplitude $w_0$ for composite plates under different fluid depth ratios (b₁/b): (a) simply supported end condition, (b) clamped supported end condition.

Effect of the fiber orientation on vibration characteristics of composite plates

In Figure 12, the numerical analytical results are obtained for composite plates with different fiber orientation $\theta$ in which orientation angle $\theta$ varies from $0^{°}$ to ${90}^{°}$ . Figure 12(a) shows the vibration frequency-the fiber orientation curves under simply supported end boundary. The results indicate that the vibration frequency of composite plates increases at first and then decreases with increasing fiber orientation $\theta$ . When the fiber orientation is $\theta ={45}^{°}$ , the vibration frequency of composite plates will reach the extreme point. Figure 12(b) plots the vibration frequency-the fiber orientation curves under clamped supported end condition. The results indicate that the effect of the fiber orientation on vibration frequency of composite plates is greater when the fluid depth ratio is relatively small and the increase of the fiber orientation leads to the decrease of the vibration frequency. Compared to the clamped supported end boundary, the vibration frequency of composite plates is more sensitive to the simply supported end boundary.

Figure 12.

Dimensionless vibration frequency $\varpi$ versus the fiber orientation $\theta$ for composite plates under different fluid depth ratios (b₁/b): (a) simply supported end condition, (b) clamped supported end condition.

Conclusions

A study on vibration characteristic analysis of composite plates in contact with a bounded fluid is carried out considering different fluid depth ratios. The added mass of composite plates under different fluid depth ratios is calculated based on the Lewis’s empirical formula. The vibration equation of composite plates is established by using the combination of Nurbs-Bézier and penalty function. A parametric analysis is conducted to evaluate the effect of elastic restrained edge stiffness, fluid depth ratio, fiber orientation, pre-stress, and initial imperfection on vibration characteristic of composite plates. Several conclusions are drawn below.

• (1) The dimensionless vibration frequency of composite plates in contact with the bounded fluid gradually increases as the elastic restrained edge stiffness increases. Moreover, the dimensionless vibration frequency of composite plates decreases as the fluid depth ratio increases.

• (2) The smaller the plate aspect ratio $\gamma$ , The higher the vibration frequency $\varpi$ . This can be explained by the effect of the added mass on the vibration characteristics of composite plates seems to be weakened by decreasing the aspect ratio.

• (3) With the raise of the compressive load scale factor ${\lambda }_c$ varies from 0 to 1, the vibration frequency of composite plates decreases first and then increases.

• (4) With the raise of the tensile load scale factor, the vibration frequency of composite plates continues to increase.

• (5) With the raise of the shear load scale factor ${\lambda }_s$ , the vibration frequency of composite plates gradually continues to decrease and the variation trend of vibration frequency is basically consistent for different fluid depth ratios.

• (6) As the initial imperfection amplitude increases, the normalized vibration frequencies of composite plates gradually increase.

• (7) The vibration frequency of composite plates increases at first and then decreases with increasing fiber orientation $\theta$ .