Identification of mechanical parameters of fiber-reinforced composites by frequency response function approximation method

Abstract

This article proposes frequency response function approximation method to identify mechanical parameters of fiber-reinforced composites. First, a fiber-reinforced composite thin plate is taken as a research object, and its natural characteristic and vibration response under pulse excitation are solved based on the Ritz method and mode superposition method, so that the theoretical calculation of frequency response function of such composite plates can be realized. Then, the identification principle based on frequency response function approximation method is illustrated and its correctness is validated by comparing with other published literature in the verification example, and the specific identification procedure is also proposed. Finally, frequency response function approximation method is applied in a study case, where the elastic moduli, Poisson’s ratios, and loss factors of the TC300 carbon/epoxy composite thin plate are identified, and the influences of boundary conditions, approximation points, total number of modes, and calculation step size on the identification accuracy and efficiency are discussed. It has been proved that the proposed method can identify mechanical parameters of fiber composite materials with high precision and efficiency.

Introduction

Fiber-reinforced composites have excellent mechanical properties, good thermal stability, and capability on weight reduction, which are widely used in aeronautics, astronautics, automotives, naval vessels, and weapons industry.^1,2 Currently, there are a large number of such composite thin walled structures, such as solar panels, aircraft engine fan blades, and large wind turbine blades. The vibration problems of such composite blades and plates, such as an excessive vibration, wear and tear, and fatigue failure, have become more and more serious due to the frequent work in harsh environments, attracting an increasing number of researchers to solve these problems.^3,4 The mechanical property parameters of fiber-reinforced composite materials usually include longitudinal elastic modulus, transverse elastic modulus, in-plane shear modulus, Poisson ratio, and loss factor. Accurately identifying these mechanical property parameters is the basis of analyzing the natural property and vibration response of composite structures, which has an important engineering and academic significance in the field of theoretical analysis, dynamic design, fault diagnosis of composite structures, and so on.

At present, great progress has been made in the identification of mechanical property parameters of fiber-reinforced composites. For example, Deobald and Gibson ⁵ presented a material parameter identification method of a fiber-reinforced composite thin plate by combining modal analysis and Rayleigh–Ritz method, and the material parameters of the composite plate, such as the elastic moduli and Poisson’s ratio values, were obtained based on the results of bending stiffness $D_{11}$ , $D_{12}$ , $D_{22}$ , and $D_{66}$ solved by the iterative algorithm. Soares et al. ⁶ proposed an indirect identification technique to predict the mechanical properties of composite plates. Based on the Mindlin plate theory and sensitivity analysis technique, the error function expressing the difference between measured frequencies of a composite plate and the corresponding numerical ones was constructed to find the desired mechanical parameters. Moussu and Nivoit ⁷ provided a non-destructive test technique to determine the material parameters of a composite plate by studying its free vibration under completely free boundary conditions, and based on the measured natural frequencies, the superposition method was employed to calculate the material parameters. Araújo et al. ⁸ presented a numerical/experimental method for the identification of material parameters of composite materials. By formulating an error function which expresses the deviation between numerical eigenvalues obtained by finite element method and experimental method, the elastic moduli of $E_1$ , $E_2$ , $G_{12}$ , $G_{13}$ , and $G_{23}$ and Poisson’s ratio values of ${\upsilon }_{12}$ and ${\upsilon }_{21}$ were obtained based on the optimization techniques. Lai and Ip ⁹ presented an iterative correction algorithm in which the modified analytical natural characteristic eventually matched those obtained by experiment to obtain the material parameters of fiber-reinforced composite plates. The Rayleigh–Ritz method was used to calculate the natural characteristic and the statistical Bayesian estimation was employed to minimize the error function of theoretical and measured eigenvalues and eigenvectors. In addition, the solution of the elastic modulus is very precise. Visscher et al. ¹⁰ presented a mixed numerical-experimental approach for the identification of the material damping properties of fiber-reinforced composite plates. The complex moduli were introduced in the theoretical modal, and the relation between the modal parameters and the material parameters was obtained using the modal strain energy method. Then, by comparing the numerical calculated modal parameters with experimental results, the elastic properties and damping properties of the composite plate could be determined. With this new technique, the problem of positing the loss tangent of Poisson’s ratio as zero is overcome. Visscher et al. ¹¹ presented an efficient method to identify the in-plane complex moduli of thin composite plates based on a mixed numerical-experimental approach, and by combining Ritz method and Bayesian estimation, the complex moduli of composite plate were derived. Hwang and Chang ¹² presented an identification method to obtain the elastic constants of composite plates by combining finite element analysis and optimum design. Based on the vibration testing, the error function of the natural frequencies obtained by experiment and the corresponding results obtained by finite element analysis were constructed; then, the iterative algorithm was used to calculate the optimum elastic constants which can satisfy the tolerance of the error function. This kind of finite element analysis method is easier to apply and very flexible. Rikards et al. ¹³ proposed a numerical-experimental method combining experiment with response surface approach to identify the elastic properties of composite plates. By minimizing the error function expressed as the difference between the natural frequencies obtained by response surface function, the expected elastic constants could be obtained. Lauwagie et al. ¹⁴ proposed a multi-model updating approach to identify the material parameters of composite plates based on the resonance frequencies of a number of freely suspended test plates. The model can ensure the uniqueness of the obtained properties. The finite element models of different test plates were simultaneously updated, and once the finite element models reproduced the measured frequencies, the updating procedure was halted, and the material properties of different layers can be retrieved from the finite element model’s database. Hwang et al. ¹⁵ proposed an elastic constants identification method of composite plates by combining a hybrid genetic algorithm and vibration testing. Meanwhile, a systematical two-step procedure was proposed to make up the problem of missing natural frequencies or natural frequencies with large errors due to the limitation in the measurement. Hence, the present method can identify the effective elastic constants of inhomogeneous composite plates Matter et al. ¹⁶ presented a numerical-experimental identification method to obtain the elastic and dissipative parameters of composite plates based on the natural frequencies, modal damping factors, and modal shapes of the specimens, and the error functions expressed as the differences of the experimental and numerical modal data were developed for the parameters estimation. And all the elastic properties and the major damping parameters can be estimated with a high accuracy.

Although the above researches have deeply investigated the identification technology of mechanical property parameters by combining theory and experiment, until now, most of them only utilize natural frequencies and modal damping ratio extracted from frequency response function (FRF) to identify elastic modulus and loss factor values of fiber-reinforced composites, and the whole information of FRF has not be fully employed to identify mechanical parameters. On one hand, FRF itself can objectively reflect the stiffness, damping, and other vibration information of the mechanical structure system, which is independent to the external excitation loads and its disturbances. On the other hand, FRF can also be easily obtained by the modal test technique, and its vertical axis represents the displacement admittance, that is, the inverse of the dynamic stiffness characteristics, which has the close relation with mechanical feature parameters of the structure system. Therefore, it is feasible that we can use the FRF data to identify the concerned elastic moduli and loss factors along the longitudinal, transverse, and shear direction of fiber-reinforced composites.

In this research, we try to find out a new identification solution for mechanical property parameters of fiber-reinforced composites, and frequency response function approximation method (FRFAM) is proposed to reach the target. In section “Theoretical calculation of FRF of fiber-reinforced composite thin plates,” first, a fiber-reinforced composite thin plate is taken as a research object, and its elastic moduli in different fiber directions are expressed as complex forms. Then, its natural characteristic and vibration response under pulse excitation are solved based on the Ritz method and mode superposition method, so that the mathematical calculation of FRF of such composite plates can be realized. In section “The identification principle based on FRFAM and its verification,” the identification principle based on FRFAM is illustrated and its correctness is validated by comparing with other published literature in the verification example. Finally, we go on to propose the specific identification procedure in section “Identification procedure,” and in section “A case study,” FRFAM is applied in a study case, where the elastic moduli, passion ratios, and loss factors of the TC300 carbon/epoxy composite thin plate are identified, and the influences of modal order, boundary conditions, and approximation points of such methods are also discussed. It has been proved that the proposed method can identify mechanical parameters of fiber composite materials with high precision and efficiency. On one hand, it can be used in the theoretical modeling and analysis process. On the other hand, it can be applied in the field of dynamic optimal design and fiber material failure prediction.

Theoretical calculation of FRF of fiber-reinforced composite thin plates

The fiber-reinforced composite thin plate is made of fiber and matrix material with n layers, as seen in Figure 1. First, set up the coordinate system $\mathrm{xoy}$ at the middle surface, and suppose the length, width, and thickness of the composite thin plate can be expressed as a, b, and h, while the fiber direction within a layer is defined as $\theta$ from the x-axis of coordinate system $\mathrm{xoy}$ . Besides, each layer of the composite plate is located at $h_{k-1}$ and $h_k$ along the z-axis with the equal thickness. In this theoretical model, “1” represents the direction parallel to the fiber, “2” represents the direction perpendicular to the fiber, and “3” represents the direction perpendicular to the 1–2 surface. Suppose the pulse excitation $F(t)$ is applied at point Q $(x_1,y_1)$ on the composite plate, and the vibration response extracted at point L $(x_2,y_2)$ is $X(t)$ . The locations of pulse excitation $F(t)$ and vibration response $X(t)$ are shown in Figure 1.

Figure 1.

A theoretical model of fiber-reinforced composite thin plate.

Considering the influence of fiber direction, the elastic moduli of fiber-reinforced composites ¹⁷ can be expressed as

$$E_1^{*}=E_1^{\prime }(1+\mathrm{i}{\eta }_1),E_2^{*}=E_2^{\prime }(1+\mathrm{i}{\eta }_2),G_{12}^{*}=G_{12}^{\prime }(1+\mathrm{i}{\eta }_{12})$$

(1)

where $E_1^{*}$ and $E_2^{*}$ represent the complex moduli paralleled and perpendicular to the fiber, respectively; $G_{12}^{*}$ represents the complex shear modulus in the 1–2 surface; and $E_1^{\prime }$ , $E_2^{\prime }$ , and $G_{12}^{\prime }$ represent the real part of complex moduli $E_1^{*}$ , $E_2^{*}$ , and complex shear modulus $G_{12}^{*}$ , respectively. And the related Poisson’s ratios which are induced by the stress in “1” and “2” direction are ${\nu }_{12}$ and ${\nu }_{21}$ .

Based on the classical lamination theory, the displacement field can be expressed as

$$u(x,y,z,t)=u_0(x,y,t)-z\frac{\partial w_0(x,y,t)}{\partial x}$$

(2a)

$$v(x,y,z,t)=v_0(x,y,t)-z\frac{\partial w_0(x,y,t)}{\partial y}$$

(2b)

$$w(x,y,z,t)=w_0(x,y,t)$$

(2c)

where $u,v,\mathrm{and}w$ represent the displacement of any point of composite plates, and $u_0,v_0,w_0$ is the displacement in the midplane. Besides, t is the time period.

Because the concerned fiber-reinforced composite thin plate is symmetrical between the middle surface, and the displacement perpendicular to the xoy plane is decoupled from the displacement in the xoy plane, the displacements $u_0$ and $v_0$ of middle surface can be ignored. At the meantime, based on the assumed displacement field of the classical laminate theory, the normal strain ${\epsilon }_z$ and shear strain ${\gamma }_{\mathrm{yz}}$ and ${\gamma }_{\mathrm{xz}}$ of composite plates are equal to zero, that is, ${\epsilon }_z={\gamma }_{\mathrm{yz}}={\gamma }_{\mathrm{xz}}=0$ . In this way, by considering the relationship between strain and displacement, the strain of any point in the theoretical model can be obtained

$${\epsilon }_x=\frac{\partial u}{\partial x}=-z\frac{{\partial }^2{w}_0}{\partial x^2}$$

(3a)

$${\epsilon }_y=\frac{\partial v}{\partial y}=-z\frac{{\partial }^2{w}_0}{\partial y^2}$$

(3b)

$${\gamma }_{\mathrm{xy}}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}=-2z\frac{{\partial }^2{w}_0}{\partial x\partial y}$$

(3c)

The bending curvatures ${\kappa }_x$ and ${\kappa }_y$ and torsion curvature ${\kappa }_{\mathrm{xy}}$ of composite thin plate on the middle surface can be expressed as

$${\kappa }_x=-\frac{{\partial }^2{w}_0}{\partial x^2},{\kappa }_y=-\frac{{\partial }^2{w}_0}{\partial y^2},{\kappa }_{\mathrm{xy}}=-2\frac{{\partial }^2{w}_0}{\partial x\partial y}$$

(4)

Then, the strain of any point of composite plate can be simplified as

$${\epsilon }_x=z{\kappa }_x,{\epsilon }_y=z{\kappa }_y,{\gamma }_{\mathrm{xy}}=z{\kappa }_{\mathrm{xy}}$$

(5)

For the orthotropic material, the stress–strain relationship in the fiber coordinates can be defined as

$$\left\{\begin{array}{c}{\sigma }_1\hfill \\ \multicolumn{1}{c}{{\sigma }_2}\\ \multicolumn{1}{c}{{\sigma }_6}\end{array}\right\}=\left[\begin{array}{ccc}{Q}_{11}^{*}\hfill & Q_{12}^{*}& 0\\ \multicolumn{1}{c}{Q_{21}^{*}}& Q_{22}^{*}& 0\\ \multicolumn{1}{c}{0}& 0& Q_{66}^{*}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_1\hfill \\ \multicolumn{1}{c}{{\epsilon }_2}\\ \multicolumn{1}{c}{{\gamma }_6}\end{array}\right\}$$

(6)

where

$$Q_{11}^{\ast }=\frac{E_1^{\ast }}{1-{\nu }_{12}^{\ast }{\nu }_{21}^{\ast }},Q_{12}^{\ast }=Q_{21}^{\ast }=\frac{{\nu }_{12}^{\ast }{E}_2^{\ast }}{1-{\nu }_{12}^{\ast }{\nu }_{21}^{\ast }},Q_{22}^{\ast }=\frac{E_2^{\ast }}{1-{\nu }_{12}^{\ast }{\nu }_{21}^{\ast }},Q_{66}^{\ast }=G_{12}^{\ast },{\nu }_{21}^{\ast }={\nu }_{12}^{\ast }\frac{E_2^{\ast }}{E_1^{\ast }}.$$

When an angle of $\theta$ exists between fiber coordinates and global coordinates, the stress–strain relationship of the kth layer of composite plates in global coordinates can be calculated using the stress–strain transformation equation, which has the following form

$${\left\{\begin{array}{c}{\sigma }_x\hfill \\ \multicolumn{1}{c}{{\sigma }_y}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xy}}}\end{array}\right\}}^{(k)}={\left[\begin{array}{ccc}{\overline{Q}}_{11}^{*}\hfill & {\overline{Q}}_{12}^{*}& {\overline{Q}}_{16}^{*}\\ \multicolumn{1}{c}{{\overline{Q}}_{12}^{*}}& {\overline{Q}}_{22}^{*}& {\overline{Q}}_{26}^{*}\\ \multicolumn{1}{c}{{\overline{Q}}_{16}^{*}}& {\overline{Q}}_{26}^{*}& {\overline{Q}}_{66}^{*}\end{array}\right]}^{(k)}\left\{\begin{array}{c}{\epsilon }_x\hfill \\ \multicolumn{1}{c}{{\epsilon }_y}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xy}}}\end{array}\right\}$$

(7)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\overline{Q}}_{11}^{*}=Q_{11}^{*}{\cos }^4{\theta }_k+2(Q_{12}^{*}+2Q_{66}^{*}){\sin }^2{\theta }_k{\cos }^2{\theta }_k+Q_{66}^{*}{\sin }^4{\theta }_k}\\ \multicolumn{1}{c}{{\overline{Q}}_{12}^{*}=(Q_{11}^{*}+Q_{22}^{*}-4Q_{66}^{*}){\sin }^2{\theta }_k{\cos }^2{\theta }_k+Q_{12}^{*}({\sin }^4{\theta }_k+{\cos }^4{\theta }_k)}\\ \multicolumn{1}{c}{{\overline{Q}}_{22}^{*}=Q_{11}^{*}{\sin }^4{\theta }_k+2(Q_{12}^{*}+2Q_{66}^{*}){\sin }^2{\theta }_k{\cos }^2{\theta }_k+Q_{22}^{*}{\cos }^4{\theta }_k}\\ \multicolumn{1}{c}{{\overline{Q}}_{16}^{*}=(Q_{11}^{*}-Q_{12}^{*}-2Q_{66}^{*})\sin {\theta }_k{\cos }^3{\theta }_k+(Q_{12}^{*}-Q_{22}^{*}+2Q_{66}^{*}){\sin }^3{\theta }_k\cos {\theta }_k}\\ \multicolumn{1}{c}{{\overline{Q}}_{26}^{*}=(Q_{11}^{*}-Q_{12}^{*}-2Q_{66}^{*}){\sin }^3{\theta }_k\cos {\theta }_k+(Q_{12}^{*}-Q_{22}^{*}+2Q_{66}^{*})\sin {\theta }_k{\cos }^3{\theta }_k}\\ \multicolumn{1}{c}{{\overline{Q}}_{66}^{*}=(Q_{11}^{*}+Q_{22}^{*}-2Q_{12}^{*}-2Q_{66}^{*}){\sin }^2{\theta }_k{\cos }^2{\theta }_k+Q_{66}^{*}({\sin }^4{\theta }_k+{\cos }^4{\theta }_k)}\end{array}$$

(8)

where k represents the kth layer of composite plates and ${\theta }_k$ represents the angle between the fiber direction in fiber coordinates and the x direction in global coordinates.

The bending and twisting moment resultants ¹⁸ in composite plates can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\left[\begin{array}{c}{M}_x\hfill \\ \multicolumn{1}{c}{M_y}\\ \multicolumn{1}{c}{M_{\mathrm{xy}}}\end{array}\right]=\sum _{k=1}^n{\underset{z_{k-1}}{\overset{z_k}{\int }}\left[\begin{array}{c}{\sigma }_x\hfill \\ \multicolumn{1}{c}{{\sigma }_y}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xy}}}\end{array}\right]}^{\left(k\right)}z\mathrm{d}z=\sum _{k=1}^n\underset{z_{k-1}}{\overset{z_k}{\int }}{\left[\begin{array}{c}{\overline{Q}}_{11}^{*}{\overline{Q}}_{12}^{*}{\overline{Q}}_{16}^{*}\hfill \\ \multicolumn{1}{c}{{\overline{Q}}_{12}^{*}{\overline{Q}}_{22}^{*}{\overline{Q}}_{26}^{*}}\\ \multicolumn{1}{c}{{\overline{Q}}_{16}^{*}{\overline{Q}}_{26}^{*}{\overline{Q}}_{66}^{*}}\end{array}\right]}^{(k)}\left\{\begin{array}{c}z{\kappa }_x\hfill \\ \multicolumn{1}{c}{z{\kappa }_y}\\ \multicolumn{1}{c}{z{\kappa }_{\mathrm{xy}}}\end{array}\right\}z\mathrm{d}z}\\ \multicolumn{1}{c}{=\left[\begin{array}{c}{D}_{11}^{*}{D}_{12}^{*}{D}_{16}^{*}\hfill \\ \multicolumn{1}{c}{D_{12}^{*}{D}_{22}^{*}{D}_{26}^{*}}\\ \multicolumn{1}{c}{D_{16}^{*}{D}_{26}^{*}{D}_{66}^{*}}\end{array}\right]\left\{\begin{array}{c}{\kappa }_x\hfill \\ \multicolumn{1}{c}{{\kappa }_y}\\ \multicolumn{1}{c}{{\kappa }_{\mathrm{xy}}}\end{array}\right\}}\end{array}$$

(9)

where $D_{\mathrm{ij}}^{*}=\frac{1}{3}\sum _{k=1}^n{\left({\overline{Q}}_{\mathrm{ij}}^{*}\right)}^{\left(k\right)}(z_k^3-z_{k-1}^3),i,j=1,2,6$

The bending strain energy stored in the composite thin plate can be expressed as

$$U=\frac{1}{2}\int \underset{A}{\int }[M_x{\kappa }_x+M_y{\kappa }_y+M_{\mathrm{xy}}{\kappa }_{\mathrm{xy}}]\mathrm{d}A$$

(10)

where A represents the area of middle surface of composite plate.

The kinetic energy of composite thin plate can be expressed as

$$T=\frac{1}{2}\underset{A}{\int }\underset{-h/2}{\overset{h/2}{\int }}\rho {\left(\frac{\partial w}{\partial t}\right)}^2\mathrm{d}A\mathrm{d}\mathrm{z}$$

(11)

Assuming that the displacement of any vibration response point L $(x_2,y_2)$ in composite thin plates is $w_0$ , as seen in Figure 1, it can be expressed as

$$w_0=W(\xi ,\eta ){\mathrm{e}}^{\mathrm{i}\omega t}$$

(12)

where $\omega$ is the excitation frequency, and $W(\xi ,\eta )$ represents modal shape function^18,19 which can be defined as

$$W(\xi ,\eta )=\sum _{i=1}^M\sum _{j=1}^N{q}_{\mathrm{ij}}{P}_i\left(\xi \right)P_j\left(\eta \right)$$

(13)

where $q_{\mathrm{ij}}$ are the eigenvectors which need to be solved, and $P_i(\xi )$ $(i=1,\dots ,M)$ and $P_j(\eta )$ $(j=1,\dots ,N)$ are the orthogonal polynomials.

Then, the orthogonal polynomials can be obtained by implementing orthogonalization operation on polynomial function, which should satisfy the boundary condition of composite plate, and these polynomials have the following expressions

$$\begin{array}{c}\multicolumn{1}{c}{P_1\left(\xi \right)=\chi \left(\xi \right),P_1\left(\eta \right)=\kappa \left(\eta \right)}\\ \multicolumn{1}{c}{P_2\left(\varphi \right)=(\varphi -H_2)P_1\left(\varphi \right)}\\ \multicolumn{1}{c}{P_i\left(\varphi \right)=(\varphi -H_i)P_{i-1}\left(\varphi \right)-V_i{P}_{i-2}\left(\varphi \right)}\\ \multicolumn{1}{c}{\varphi =\xi ,\eta ,i>2}\end{array}$$

where $H_i$ and $V_i$ are the coefficient functions ¹⁸ and their expressions can be written as

$$H_i=\frac{\underset{0}{\overset{1}{\int }}W\left(\varphi \right){\left[P_{i-1}\left(\varphi \right)\right]}^2\varphi \mathrm{d}\varphi }{\underset{0}{\overset{1}{\int }}W\left(\varphi \right){\left[P_{i-1}\left(\varphi \right)\right]}^2\mathrm{d}\varphi }$$

(14a)

$$V_i=\frac{\underset{0}{\overset{1}{\int }}W\left(\varphi \right)P_{i-1}\left(\varphi \right)P_{i-2}\left(\varphi \right)\varphi \mathrm{d}\varphi }{\underset{0}{\overset{1}{\int }}W\left(\varphi \right){\left[P_{i-2}\left(\varphi \right)\right]}^2\mathrm{d}\varphi },\varphi =\xi ,\eta$$

(14b)

where $W(\varphi )$ is the weighting function and usually it can be set to $W(\varphi )=1$ . Besides, $\chi (\xi )$ and $\kappa (\eta )$ are polynomial functions which can satisfy the different boundary conditions, such as the clamped, simply support, and free boundary, and they can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\chi \left(\xi \right)={\xi }^{\alpha }{(1-\xi )}^{\beta },\kappa \left(\eta \right)={\eta }^{\gamma }{(1-\eta )}^{\tau }}\\ \multicolumn{1}{c}{\xi =\frac{x}{a},\eta =\frac{y}{b}}\end{array}$$

(15)

For completely free boundary condition, we can set $\alpha =0$ , $\beta =0$ , $\gamma =0$ , $\tau =0$ , and for cantilever boundary condition, we can set $\alpha =2$ , $\beta =0$ , $\gamma =0$ , $\tau =0$ . Meanwhile, because of its small influence on natural characteristic, the damping is ignored; therefore, only the real parts of elastic moduli in stress–strain relationship are used.

Then, substituting equation (12) into equations (10) and (11), we can obtain the maximum bending strain energy $U_{max}$ stored in the plate and the maximum kinetic energy $T_{max}$ , which can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{U_{max}=\frac{1}{2}\int \underset{A}{\int }{D}_{11}^{\prime }{\left(\frac{{\partial }^{2}W}{\partial x^2}\right)}^2+2D_{12}^{\prime }\frac{{\partial }^{2}W}{\partial x^2}\frac{{\partial }^{2}W}{\partial y^2}+D_{22}^{\prime }{\left(\frac{{\partial }^{2}W}{\partial y^2}\right)}^2}\\ \multicolumn{1}{c}{+4(D_{16}^{\prime }\frac{{\partial }^{2}W}{\partial x^2}+D_{26}^{\prime }\frac{{\partial }^{2}W}{\partial y^2})\frac{{\partial }^{2}W}{\partial x\partial y}+4D_{66}^{\prime }{\left(\frac{{\partial }^{2}W}{\partial x\partial y}\right)}^2]\mathrm{d}A}\end{array}$$

(16)

$$T_{max}=\frac{\rho h{\omega }^2}{2}\int \underset{A}{\int }{\left(W\right)}^2\mathrm{d}A$$

(17)

where $D_{\mathrm{ij}}^{\prime }$ represents the real part of complex stiffness coefficient in equation (9).

Next, define the Lagrangian energy function L as

$$L=T_{max}-U_{max}$$

(18)

By minimizing partial derivative of the Lagrangian energy function L with respect to $q_{\mathrm{ij}}$ in the following equation

$$\frac{\partial L}{\partial q_{\mathrm{ij}}}=0,i=1,2,\dots ,M,j=1,2,\dots ,N$$

(19)

Substituting equation (18) into equation (19), the eigenvalue equation can be obtained as

$$(\mathit{K}-{\omega }^2\mathit{M})\mathit{q}=0$$

(20)

where K and M represent stiffness matrix and mass matrix of the structural system, respectively, and $\mathit{q}=(q_{11},q_{12},\dots ,q_{\mathrm{ij}}{)}^T$ is the eigenvector.

In order to make equation (20) have nonzero solution or nontrivial solution, the determinant of the coefficient matrix should be equal to zero

$$|\mathit{K}-{\omega }^2\mathit{M}|=0$$

(21)

By solving equation (21), the natural frequencies of the composite thin plate can be obtained. Then, substituting the eigenvector corresponding to a certain natural frequency into equation (13), the modal shape corresponding to this natural frequency can be easily obtained. Repeating the steps above, the all concerned modal shape will be obtained.

After the natural frequencies and modal shapes of the composite plate are solved, assume the pulse excitation $F(x,y,t)$ and the vibration response^19,20 $X(x,y,t)$ to have the following expression

$$\begin{array}{c}\multicolumn{1}{c}{F(x,y,t)=f\left(t\right)\delta (x-x_1)\delta (y-y_1)}\\ \multicolumn{1}{c}{f\left(t\right)=\{\begin{array}{c}{f}_0\sin (\omega t)\hfill \\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c},\hfill \\ \multicolumn{1}{c}{,}\end{array}\begin{array}{c}0\le t\le t_1\hfill \\ \multicolumn{1}{c}{t>t_1}\end{array}}\end{array}$$

(22)

$$X(x,y,t)=\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{W}_{\mathrm{ij}}(x,y)T_{\mathrm{ij}}\left(t\right)$$

(23)

where $\omega$ is the excitation circular frequency, $t_1$ is the excitation time, $W_{\mathrm{ij}}(x,y)$ is the modal shape, and $T_{\mathrm{ij}}$ is the mode component.

Considering the dynamic balance of composite plate and ignoring damping and inertia moment, the following equations can be obtained

$$\frac{\partial M_x}{\partial x}+\frac{\partial M_{\mathrm{xy}}}{\partial y}-Q_x=0$$

(24a)

$$\frac{\partial M_{\mathrm{xy}}}{\partial x}+\frac{\partial M_y}{\partial y}-Q_y=0$$

(24b)

$$\frac{\partial Q_x}{\partial x}+\frac{\partial Q_y}{\partial y}+F-\rho h\frac{{\partial }^2{w}_0}{\partial t^2}=0$$

(24c)

where $\rho$ represents the density of composite plate; $Q_x$ and $Q_y$ represent shear force in surface which is vertical to x axis and y axis; and $M_x$ , $M_y$ , and $M_{\mathrm{xy}}$ are the bending and twisting moment resultants which have been shown in equation (9).

Substituting equation (9) into equation (24), the forced vibration differential equation of composite thin plate without considering damping can be simplified as

$$\begin{array}{c}\multicolumn{1}{c}{D_{11}^{\prime }\frac{{\partial }^{4}X}{\partial x^4}+2(D_{12}^{\prime }+2D_{66}^{\prime })\frac{{\partial }^{4}X}{\partial x^2\partial y^2}+4D_{16}^{\prime }\frac{{\partial }^{4}X}{\partial x^3\partial y}+D_{22}^{\prime }\frac{{\partial }^{4}X}{\partial y^4}}\\ \multicolumn{1}{c}{+4D_{26}^{\prime }\frac{{\partial }^{4}X}{\partial x\partial y^3}+\rho h\frac{{\partial }^{2}X}{\partial t^2}=F}\end{array}$$

(25)

According to modal shape superposition method, the displacement response of a continuous system can be expressed as the series of modal shape function, and the partial differential equations under physical coordinates can be converted to the second-order ordinary differential equations under generalized coordinates. Thus, the continuous system would be simplified as the single degree of freedom system.

Substituting equation (23) into equation (25), the following equation is obtained

$$\begin{array}{c}\multicolumn{1}{c}{D_{11}^{\prime }\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{T}_{\mathrm{ij}}\frac{{\partial }^4{W}_{\mathrm{ij}}}{\partial x^4}+2(D_{12}^{\prime }+2D_{66}^{\prime })\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{T}_{\mathrm{ij}}\frac{{\partial }^4{W}_{\mathrm{ij}}}{\partial x^2\partial y^2}}\\ \multicolumn{1}{c}{+4D_{16}^{\prime }\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{T}_{\mathrm{ij}}\frac{{\partial }^4{W}_{\mathrm{ij}}}{\partial x^3\partial y}+D_{22}^{\prime }\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{T}_{\mathrm{ij}}\frac{{\partial }^4{W}_{\mathrm{ij}}}{\partial y^4}}\\ \multicolumn{1}{c}{+4D_{26}^{\prime }\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{T}_{\mathrm{ij}}\frac{{\partial }^4{W}_{\mathrm{ij}}}{\partial x\partial y^3}+\rho h\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{W}_{\mathrm{ij}}\frac{d^2{T}_{\mathrm{ij}}}{d{t}^2}=F}\end{array}$$

(26)

Besides, according to displacement variation principle, the modal shape function $W_{\mathrm{ij}}$ should satisfy the following equation

$$\begin{array}{c}\multicolumn{1}{c}{D_{11}^{\prime }\frac{{\partial }^4{W}_{\mathrm{ij}}}{\partial x^4}+2(D_{12}^{\prime }+2D_{66}^{\prime })\frac{{\partial }^4{W}_{\mathrm{ij}}}{\partial x^2\partial y^2}+4D_{16}^{\prime }\frac{{\partial }^4{W}_{\mathrm{ij}}}{\partial x^3\partial y}+D_{22}^{\prime }\frac{{\partial }^4{W}_{\mathrm{ij}}}{\partial y^4}}\\ \multicolumn{1}{c}{+4D_{26}^{\prime }\frac{{\partial }^4{W}_{\mathrm{ij}}}{\partial x\partial y^3}=\rho h({\omega }_{\mathrm{ij}}{)}^2{W}_{\mathrm{ij}}}\end{array}$$

(27)

Substituting equation (27) into equation (26) and through simplifying, the following equation can be derived

$$\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }(\rho h\frac{{\mathrm{d}}^2{T}_{\mathrm{ij}}}{\mathrm{d}{t}^2}+\rho h{\left({\omega }_{\mathrm{ij}}\right)}^2{T}_{\mathrm{ij}})W_{\mathrm{ij}}=F$$

(28)

Then, multiply $W_{\mathrm{kl}}(x,y)(k,l=1,2,3,\dots )$ on both sides of equation (28) and perform integral operation along x–y plane; the following equation can be obtained

$$\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }\int \underset{A}{\int }(\rho h\frac{{\mathrm{d}}^2{T}_{\mathrm{ij}}}{\mathrm{d}{t}^2}+\rho h{\left({\omega }_{\mathrm{ij}}\right)}^2{T}_{\mathrm{ij}})W_{\mathrm{ij}}{W}_{\mathrm{kl}}\mathrm{d}A=\int \underset{A}{\int }F{W}_{\mathrm{kl}}\mathrm{d}A$$

(29)

Furthermore, by utilizing the orthogonality of modal shape of plates, which has the following form

$$\int \underset{A}{\int }\rho h{W}_{\mathrm{ij}}{W}_{\mathrm{kl}}\mathrm{d}A=0(i\ne j\mathrm{or}k\ne l)$$

(30)

The generalized vibration differential equation without damping can be expressed as

$$\frac{{\mathrm{d}}^2{T}_{\mathrm{ij}}\left(t\right)}{\mathrm{d}{t}^2}+{\left({\omega }_{\mathrm{ij}}\right)}^2{T}_{\mathrm{ij}}\left(t\right)=\frac{P_{\mathrm{ij}}\left(t\right)}{M_{\mathrm{ij}}}$$

(31)

where $P_{\mathrm{ij}}(t)$ and $M_{\mathrm{ij}}$ are generalized force and generalized mass, respectively, which have the following form

$$P_{\mathrm{ij}}\left(t\right)=\int \underset{A}{\int }F{W}_{\mathrm{ij}}(x,y)\mathrm{d}A$$

(32)

$$M_{\mathrm{ij}}=\int \underset{A}{\int }\rho h{\left(W_{\mathrm{ij}}(x,y)\right)}^2\mathrm{d}A$$

(33)

Meanwhile, considering the following property of the function of $\delta$

$$\{\begin{array}{c}\underset{-\infty }{\overset{+\infty }{\int }}\delta (x-x^{*})\mathrm{d}x=1\hfill \\ \multicolumn{1}{c}{\delta (x-x^{*})=0(x\ne x^{*})}\end{array}$$

(34)

The generalized force $P_{\mathrm{ij}}(t)$ can be simplified as

$$P_{\mathrm{ij}}\left(t\right)=f\left(t\right)W_{\mathrm{ij}}(x_1,y_1)$$

(35)

Similarly, under the assumption of small damping, the generalized vibration differential equation with damping can be expressed as

$$\frac{{\mathrm{d}}^2{T}_{\mathrm{ij}}\left(t\right)}{\mathrm{d}{t}^2}+2{\zeta }_i{\omega }_{\mathrm{ij}}\frac{\mathrm{d}{T}_{\mathrm{ij}}}{\mathrm{d}t}+{\left({\omega }_{\mathrm{ij}}\right)}^2{T}_{\mathrm{ij}}\left(t\right)=\frac{P_{\mathrm{ij}}\left(t\right)}{M_{\mathrm{ij}}}$$

(36)

where ${\zeta }_i$ is the ith modal damping ratio of the composite plate.

At the zero initial condition, the solution of equation (36) can be expressed as Duhamel integral form, which has the following expression

$$T_{\mathrm{ij}}(t)=\frac{W_{\mathrm{ij}}(x_1,y_1)}{{\omega }_{\mathrm{d}}{M}_{\mathrm{ij}}}\underset{0}{\overset{t}{\int }}f(\tau ){\mathrm{e}}^{-{\zeta }_i{\omega }_{\mathrm{ij}}(t-\tau )}\sin {\omega }_{\mathrm{d}}(t-\tau )\mathrm{d}\tau$$

(37)

where ${\omega }_{\mathrm{d}}$ is the circular frequency which can be expressed as ${\omega }_{\mathrm{d}}=\sqrt{1-{\zeta }_i^2}{\omega }_{\mathrm{mn}}$ .

Equation (37) can be solved using Simpson numerical integration, and substituting its results into equation (23), the vibration response $X(x,y,t)$ of composite plate in equation (38) under pulse excitation $F(t)$ can be obtained using modal shape superposition method

$$X(x,y,t)=\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }{W}_{\mathrm{ij}}(x,y)T_{\mathrm{ij}}\left(t\right)$$

(38)

Finally, the time-domain wave of excitation force and displacement response based on self-designed MATLAB program can be plotted, and the fast Fourier transformation (FFT) operation is applied to excitation and response data. Consequently, the Fourier function of $F(t)$ and $X(t)$ of composite plate can be transformed as $F(\omega )$ and $X\left(\omega \right)$ , and according to the definition of FRF, the FRF $H\left(\omega \right)$ between excitation point Q and response point L can be expressed as

$$H\left(\omega \right)=\frac{X\left(\omega \right)}{F\left(\omega \right)}$$

(39)

The identification principle based on FRFAM and its verification

The identification principle

In section “Theoretical calculation of FRF of fiber-reinforced composite thin plates,” the FRF of composite thin plates under pulse excitation is obtained through theoretical method. However, in most cases, the experimental modal test is mainly conducted to obtain the FRF expressed in equation (38). If the experiment FRF data can be first obtained, while the theoretical FRF data can also be calculated by theoretical method, then we can modify the theoretical FRF data with the experimental FRF data as a reference. Because the mechanical parameters of fiber-reinforced composites, such as $E_1$ , $E_2$ , $G_{12}$ , ${\nu }_{12}$ , ${\eta }_1$ , ${\eta }_2$ , and ${\eta }_{12}$ , are the main factors that affect the theoretical FRF, we can identify these parameters by enabling the theoretical FRF approximate to the experimental one, which is named as “frequency response function approximation method.” In the following section, its identification principle will be described in detail.

The schematic diagram to obtain the FRF of composite thin plates under pulse excitation by the above two methods is shown in Figure 2, which mainly contains three mapping forms: (1) the mapping between excitation force in the experiment produced by the hammer and pulse excitation force in the theoretical calculation, (2) the mapping between experimental vibration response and theoretical vibration response, and (3) the mapping between experimental test procedure and theoretical analysis procedure. By applying pulse excitation force at same position of the composite plate and obtaining the resulting vibration response at same location, we can modify the mechanical parameters using FRFAM. Finally, when the theoretical FRF can approximate the experimental FRF closely, the input mechanical parameters in the theoretical calculation can be equivalent to the real mechanical parameters of the composite plate. Figure 3 gives the schematic diagram of identification principle of this method.

Figure 2.

The schematic diagram to obtain FRF of composite thin plate.

Figure 3.

The schematic diagram of identification principle.

First, the experimental FRF and theoretical FRF with the two methods described in the previous paragraph should be obtained, which can be plotted on the same graph, as shown in Figure 3. Then, the whole concerned frequency range can be divided into multiple frequency ranges, and in each frequency range, two to three modes of composite thin plate should be included. Finally, we can use the approaching calculation techniques to construct the frequency relative error function $e_{\mathrm{fre}}$ between the experimental natural frequencies and the theoretical natural frequencies in each frequency range, which can have the following form

$$e_{\mathrm{fre}}=\sum _{i=1}^R{\left(\frac{\left|\Delta f_i\right|}{{\hat{f}}_i}\right)}^2$$

(40)

where R represents the number of modes in the whole frequency range, $\Delta f_i$ represents the difference between the ith natural frequency obtained by experiment and theoretical calculation, and ${\hat{f}}_i$ is the ith natural frequency obtained by experiment.

Then, taking the average material parameters as the center, such as $E_1^0,E_2^0,G_{12}^0,{\nu }_{12}^0$ , which are usually provided by the composite material manufacturer, and with the consideration of parameters error $R_{\mathrm{err}}$  = 10%–20%, the range of material parameters can be determined as follows

$$\begin{array}{c}\multicolumn{1}{c}{E_1^0(1-R_{\mathrm{err}})\le E_1\le E_1^0(1+R_{\mathrm{err}})}\\ \multicolumn{1}{c}{E_2^0(1-R_{\mathrm{err}})\le E_2\le E_2^0(1+R_{\mathrm{err}})}\\ \multicolumn{1}{c}{G_{12}^0(1-R_{\mathrm{err}})\le G_{12}\le G_{12}^0(1+R_{\mathrm{err}})}\\ \multicolumn{1}{c}{{\nu }_{12}^0(1-R_{\mathrm{err}})\le {\nu }_{12}\le {\nu }_{12}^0(1+R_{\mathrm{err}})}\end{array}$$

(41)

Furthermore, select an appropriate step size g (e.g. g = 1%) in the above range and construct iteration vectors of material parameters, such as ${\mathit{E}}_1,{\mathit{E}}_2,{\mathit{G}}_{12},{\mathbf{\nu }}_{12}$ , which can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{E}}_1=[E_1^1{E}_1^2\cdots E_1^n]}\\ \multicolumn{1}{c}{{\mathit{E}}_2=[E_2^1{E}_2^2\cdots E_2^n]}\\ \multicolumn{1}{c}{{\mathit{G}}_{12}=[G_{12}^1{G}_{12}^2\cdots G_{12}^n]}\\ \multicolumn{1}{c}{{\mathbf{\nu }}_{12}=[{\nu }_{12}^1{\nu }_{12}^2\cdots {\nu }_{12}^n]}\end{array}$$

(42)

where $Z^1=Z^0(1-R_{\mathrm{err}}),Z^2=Z^0(1-R_{\mathrm{err}})+2g{R}_{\mathrm{err}}{Z}^0,Z^n=Z^0(1-R_{\mathrm{err}})+2g(n-1)R_{\mathrm{err}}{Z}^0,Z=E_1,E_2,G_{12},{\nu }_{12}$ .

By iteratively calculating the material parameters in a permutation and combination manner, we can obtain the optimum estimation results of $E_1,E_2,G_{12},{\nu }_{12}$ when the frequency relative error function $e_{\mathrm{fre}}$ gets the minimum value.

Similarly, by extracting the multiple resonance peaks from the experimental FRF and theoretical FRF, we can construct the admittance relative error function $e_{\mathrm{rec}}$ , which has the following form

$$e_{\mathrm{rec}}=\sum _{r=1}^S{\left(\frac{\left|\Delta H_{\mathrm{ir}}\right|}{\left|{\hat{H}}_{\mathrm{ir}}\right|}\right)}^2$$

(43)

where i is the ith resonance peak, S is the number of approaching points on the ith resonance peak, $\left|\Delta H_{\mathrm{ir}}\right|$ represents the difference of amplitude between the ith experimental FRF and theoretical FRF at approaching point r, and $\left|{\hat{H}}_{\mathrm{ir}}\right|$ represents the amplitude of the ith experimental FRF at approaching point r.

Next, take the maximum loss factor $\eta$  = 0.03 and construct iteration vectors of loss factors, such as ${\mathbf{\eta }}_1$ , ${\mathbf{\eta }}_2$ , and ${\mathbf{\eta }}_{12}$ with an appropriate step size g (e.g. g = 1%) in a range of 0– $\eta$ , which can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\eta }}_1=[{\eta }_1^1{\eta }_1^2\cdots {\eta }_1^n]}\\ \multicolumn{1}{c}{{\mathbf{\eta }}_2=[{\eta }_2^1{\eta }_2^2\cdots {\eta }_2^n]}\\ \multicolumn{1}{c}{{\mathbf{\eta }}_{12}=[{\eta }_{12}^1{\eta }_{12}^2\cdots {\eta }_{12}^n]}\end{array}$$

(44)

where ${\eta }_{\mathrm{ij}}^1=0,{\eta }_{\mathrm{ij}}^2=g\eta ,\dots ,{\eta }_{\mathrm{ij}}^n=(n-1)g\eta$ .

According to modal strain energy method, the relations between modal loss factors and loss factors in different fiber direction can be expressed as

$${\eta }_i=\frac{{\eta }_1{U}_1+{\eta }_2{U}_2+{\eta }_{12}{U}_{12}}{U}$$

(45)

where U is the total strain energy of the composite plate, $U_1$ and $U_2$ represent the strain energy paralleled and perpendicular to the fiber, and $U_{12}$ represents the strain energy in the 1–2 surface, which have the following form

$$\begin{array}{c}\multicolumn{1}{c}{U_1=\sum _{k=1}^n\frac{1}{2}\underset{x=0}{\overset{a}{\int }}\underset{y=0}{\overset{b}{\int }}\underset{h_{k-1}}{\overset{h_k}{\int }}{\sigma }_1^k{\epsilon }_1^k\mathrm{d}x\mathrm{d}y\mathrm{d}z}\\ \multicolumn{1}{c}{U_2=\sum _{k=1}^n\frac{1}{2}\underset{x=0}{\overset{a}{\int }}\underset{y=0}{\overset{b}{\int }}\underset{h_{k-1}}{\overset{h_k}{\int }}{\sigma }_2^k{\epsilon }_2^k\mathrm{d}x\mathrm{d}y\mathrm{d}z}\end{array}$$

(46)

$$U_{12}=\sum _{k=1}^n\frac{1}{2}\underset{x=0}{\overset{a}{\int }}\underset{y=0}{\overset{b}{\int }}\underset{h_{k-1}}{\overset{h_k}{\int }}{\sigma }_6^k{\gamma }_6^k\mathrm{d}x\mathrm{d}y\mathrm{d}z$$

(46)

Then, according to the relation between modal damping ratios and modal loss factors, the modal damping ratio ${\zeta }_i$ can be expressed as

$${\zeta }_i=\frac{{\eta }_i}{2}$$

(47)

Finally, by iteratively calculating the loss factors in a permutation and combination manner, the optimum estimation results of ${\eta }_1$ , ${\eta }_2$ , and ${\eta }_{12}$ can be obtained when the admittance relative error function $e_{\mathrm{rec}}$ gets the minimum value.

Verification example

In order to verify the correctness of the method proposed in this article, we employ FRFAM to identify the mechanical parameters of graphite/epoxy composite thin plates in Han and Lee ²¹ with laminate configuration of $[0°/\pm 45°/90°{]}_s$ under cantilever boundary condition. The length, width, and thickness of such composite plates are 220, 180, and 0.82 mm, and the corresponding mechanical parameters used in Han and Lee ²¹ are listed in Table 1. Besides, the FRF curve is obtained using the impact hammer and the laser displacement meter, and the measured natural frequencies and modal damping ratios are listed in Table 2.

Table 1.

The mechanical parameters given in Han and Lee ²¹ and identified in this research by FRFAM.

Mechanical parameters	E1 (GPa)	E2 (GPa)	G12 (GPa)	ν 12	${\eta }_1$	${\eta }_2$	${\eta }_{12}$
Han and Lee 21 A	119	8.67	5.18	0.31	–	–	–
Present B	120	8.5	5.5	0.33	0.0025	0.0085	0.0495
Error (%) $|B-A|/A$	0.84	1.96	6.18	6.45	–	–	–

Table 2.

The natural frequencies and modal damping ratios measured in Han and Lee ²¹ and calculated in this research by FRFAM.

Modal order	Natural frequency			Modal damping ratio
	Han and Lee 21 A (Hz)	Present B (Hz)	Error (%) $|B-A|/A$	Han and Lee 21 C (%)	Present D (%)	Error (%) $|D-C|/C$
1	18.29	20.01	9.40	0.236	0.26	10.17
2	44.77	43.79	2.19	0.225	0.23	2.22
3	114.58	116.46	1.64	0.503	0.52	3.38
4	133.45	131.08	1.78	0.233	0.22	5.58
5	168.51	168.21	0.18	0.442	0.46	1.8

By utilizing the above data of dimensions, laminate configuration, and the experimental FRF of composite plate, the calculated FRF curve can be obtained by FRFAM. In addition, the identified mechanical parameters as well as the calculated errors are also listed in Table 1. Finally, the identified mechanical parameters are input into the theoretical model to recalculate the natural frequencies and modal damping ratios, as seen in Table 2, and the errors between the calculated frequency and damping results in this research and the ones given in Han and Lee ²¹ are also listed in Table 2.

From Table 1, we can find that the maximum error between the identified material parameters and the parameters given in Han and Lee ²¹ is no more than 6.45%. In addition, from Table 2, we can compare the natural frequency and modal damping ratio results calculated by this research with the results given in Han and Lee, ²¹ then we can find out that the maximum errors of natural frequencies and modal damping ratios are less than 9.40% and 10.17%, respectively, which are also within an acceptable range. Therefore, the correctness and practicability of the proposed method can be verified.

Identification procedure

In this section, the identification procedure of mechanical parameters of fiber-reinforced composites is summarized, which is realized based on self-designed MATLAB program and can be divided into the following steps.

Solve natural characteristic and vibration response

First, a fiber-reinforced composite thin plate is taken as a research object, and its theoretical model is established based on the classical lamination theory. Meanwhile, its elastic moduli in different fiber directions are expressed as complex forms, and the natural characteristic and vibration response can be solved using Ritz method and mode superposition method, respectively.

Solve FRF of the composite thin plate under pulse excitation

Using equation (38) to calculate vibration response $X(t)$ in any point of composite thin plate under pulse excitation, the theoretical FRF $H(\omega )$ between excitation point Q and response point L can be solved as in equation (39).

Measure FRF of composite thin plate with the hammer

If the FRF measurement is done under the free boundary condition, it is necessary to use a rubber rope to hang the tested plate specimen, and if the FRF measurement is done under the constrained boundary condition, it is needed to ensure that the composite plate is effectively clamped by the fixture tools. In addition, in order to get good FRF results, lightweight accelerometer should be used to get response signal to reduce the effects of the additional mass and stiffness, and the FRF results and the coherence function results need to be tested simultaneously. It is suggested that their final values should be the average results, which are measured at least three times. Then, the coherence function can be used to evaluate the test accuracy of FRF data. If the coherence coefficients are above 0.9, the corresponding FRF results are accepted, otherwise the measurement should be re-conducted several times.

Identify the elastic coefficients $E_1,E_2,G_{12},\mathrm{and}{\nu }_{12}$ by the natural frequency approaching technique

In this step, we should separate the concerned frequency range into several sub-ranges and identify the theoretical and experimental natural frequencies in each sub-range, and then select an appropriate step size to construct iteration vectors of material parameters, so that the frequency relative error function $e_{\mathrm{fre}}$ between the theoretical and experimental natural frequency results can be established. Finally, by iteratively calculating the material parameters in a permutation and combination manner, the optimum estimation results of $E_1,E_2,G_{12},\mathrm{and}{\nu }_{12}$ can be obtained when the frequency relative error function $e_{\mathrm{fre}}$ gets the minimum value.

Identify the loss factors ${\eta }_1$ , ${\eta }_2$ , and ${\eta }_{12}$ by the displacement admittance approaching technique

This step is similar to the step (4), yet the approximating object is changed to the displacement admittance of theoretical and experimental FRF. First, set the number of approaching points in each sub-range of FRF data and construct iteration vectors of loss factors with an appropriate step size, so that the displacement admittance relative error function $e_{\mathrm{rec}}$ can be established. Then, by iteratively calculating the loss factors in a permutation and combination manner, the optimum estimation results of ${\eta }_1$ , ${\eta }_2$ , and ${\eta }_{12}$ can be obtained when the admittance relative error function $e_{\mathrm{rec}}$ gets the minimum value.

A case study

In this section, a TC300 carbon/epoxy composite plate is taken as a research object and its mechanical parameters are identified by FRFAM. The studied plate specimen is symmetrically laid, which is cut from the fiber composite material with total of 21 layers and laminate configuration of ((0/90)₅/0/(90/0)₅). The length, width, and thickness are 230, 130, and 2.36 mm, respectively, and each layer has the same thickness. Besides, its material parameter is provided by the manufacturer, with the longitudinal elastic modulus E₁ = 105.0 GPa, transverse elastic modulus E₂ = 7.9 GPa, shear modulus G₁₂ = 4.2 GPa, Poisson’s ratio ${\nu }_{12}=0.3$ , and density $\rho =1780\mathrm{kg}/{\mathrm{m}}^3$ .

Test system

In order to obtain the FRF of fiber-reinforced composite plates under pulse excitation, the FRF test system is set up. The instruments used in the test are as follows: (1) BK 4517-001 lightweight accelerometer, (2) LMS SCADAS 16-channel Mobile Front-End, (3) Dell M7700 computer and LMS TestLab 10B software, and (4) PCB 086C01 Hammer. After the repeated testing and comparison, the pulse excitation is applied at the point Q, which is 100 mm from the left edge and 30 mm from the upside edge of the composite plate. In addition, the accelerometer is glued firmly at response point L by super glue 502, which is 40 mm from the left edge and 30 mm from the downside edge of the composite plate. Then, the following setups and parameters are chosen: (1) sampling frequency of 3200 Hz, (2) frequency resolution of 0.125 Hz, and (3) force-exponential window for excitation signal and exponential window for response signals. Consequently, we can employ pulse excitation technique to carry out modal test, and by analyzing the response signal in frequency domain with LMS software, the natural frequency and the damping ratio of each mode can be obtained with the half power bandwidth method. The real picture of the FRF test of the composite thin plate under different boundary conditions can be seen in Figure 4(a) and (b).

Figure 4.

The real picture of the FRF test of the composite thin plate under different boundary conditions: (a) free boundary and (b) cantilever boundary.

Identification results of mechanical parameters under different boundary conditions

In this section, the FRFs of the composite plate under the free and cantilever boundary conditions are measured. First, pulse excitation is applied on the tested plate under the free boundary condition by the hammer, and the excitation signal and response signal are recorded at the same time by the data acquisition instrument. Since the lightweight accelerometer is used, we can only get FRF in the form of acceleration admittance. Then, after applying the integral operation, we can get FRF results in the concerned form of displacement admittance. Finally, we can use the approaching calculation techniques to draw the theoretical FRF of composite thin plate repeatedly, and when the frequency relative error function $e_{\mathrm{fre}}$ and admittance relative error function $e_{\mathrm{rec}}$ can get the minimum value (e.g. in the range of 5%–10%), then the optimum mechanical parameters are obtained. Figure 5 gives the theoretical and experimental FRFs in different frequency ranges which contain a total of seven modes of composite thin plate.

Figure 5.

Theoretical and experimental FRFs in different frequency ranges under the free boundary condition: (a) frequency range of 100–500 Hz, (b) frequency range of 600–900 Hz, and (c) frequency range of 900–1500 Hz.

Then, by applying the same test technique, the FRF data under the cantilever boundary condition can be obtained. Consequently, the theoretical and experimental FRF in different frequency ranges can be plotted using the approaching calculation techniques, as seen in Figure 6. If the constructed frequency relative error function $e_{\mathrm{fre}}$ and admittance relative error function $e_{\mathrm{rec}}$ can get the minimum value (e.g. in the range of 5%–10%), then the optimum mechanical parameters under cantilever boundary condition are obtained.

Figure 6.

Theoretical and experimental FRFs in different frequency ranges under the cantilever boundary condition: (a) frequency range of 0–200 Hz, (b) frequency range of 200–500 Hz, and (c) frequency range of 600–1200 Hz.

Finally, set five approximation points and calculation step size of 1% in the self-designed MATLAB program, and the mechanical parameters of fiber-reinforced composites, such as $E_1$ , $E_2$ , $G_{12}$ , ${\nu }_{12}$ , ${\eta }_1$ , ${\eta }_2$ , and ${\eta }_{12}$ can be identified by FRFAM, as listed in Table 3. The identified material parameter results are also compared with the ones provided by the manufacturer (which cannot provide loss factor values), and the corresponding identified errors are listed in Table 3. It can be seen that the maximum error under the two different boundary conditions is less than 9.8%, which is within an acceptable range. Besides, the identified loss factors along the longitudinal, transverse, and shear direction of the fiber-reinforced composite under the free boundary is smaller than the ones obtained under the cantilever boundary, which is due to the frictional effects in the clamping fixture. However, the identified errors of elastic moduli and Poisson’s ratio under the cantilever boundary condition are larger than the ones under the free boundary condition, which may be caused by the following reasons: (1) the FRF approaching performance by FRFAM under the cantilever boundary is not that good, especially in the higher modes, because the calculation errors of higher natural frequencies are greater than that under free boundary and (2) there are some measurement errors in the experimental test, because the FRF results are not obtained under the fully constrained cantilever boundary (the clamping fixture is not fully clamped), which will lead to the inexact FRF results in the theoretical calculation.

Table 3.

The identified mechanical parameters of fiber-reinforced composite by FRFAM under the two different boundary conditions and their errors.

Mechanical parameters	$E_1$ (GPa)	$E_2$ (GPa)	$G_{12}$ (GPa)	${\nu }_{12}$	${\eta }_1$	${\eta }_2$	${\eta }_6$
Manufacturer A	105.0	7.9	4.2	0.3	–	–	–
Free boundary B	112.8	8.2	4.3	0.3	0.0015	0.0055	0.0266
Cantilever boundary C	115.3	8.5	4.3	0.3	0.0017	0.0056	0.0268
Error (%) $|B-A|/A$	7.4	3.7	2.4	0.0	–	–	–
Error (%) $|C-A|/A$	9.8	7.5	2.4	0.0	–	–	–

Identification results under different total number of modes

Since the identification results of mechanical parameters of fiber-reinforced composite under free boundary condition are more accurate, in the following section, the influence of different total number of modes on the identification results under such boundaries are discussed. Still use the same setting in section “Identification results of mechanical parameters under different boundary conditions” with five approximation points and calculation step size of 1% in the MATLAB program, yet choose FRF data in the three different numbers of modes, such as 3 modes, 5 modes, and 7 modes, to identify mechanical parameters of the fiber-reinforced composite by FRFAM. Table 4 lists the identified results, the corresponding errors, and time expenditures.

Table 4.

The identified mechanical parameters of fiber-reinforced composite by FRFAM under different total number of modes, the corresponding errors, and time expenditures.

The total number of modes	$E_1$ (GPa)	$E_2$ (GPa)	$G_{12}$ (GPa)	${\nu }_{12}$	${\eta }_1$	${\eta }_2$	${\eta }_6$	Calculation time (s)
Manufacturer A	105.0	7.9	4.2	0.3	–	–	–	–
3 modes B	107.8	7.7	4.3	0.3	0.0015	0.0054	0.0264	285.4
5 modes C	110.6	7.7	4.3	0.3	0.0016	0.0055	0.0267	326.4
7 modes D	112.8	8.2	4.3	0.3	0.0015	0.0055	0.0266	436.1
Error (%) $|B-A|/A$	2.7	2.5	2.4	0.0	–	–	–	–
Error (%) $|C-A|/A$	5.3	2.5	2.4	0.0	–	–	–	–
Error (%) $|D-A|/A$	7.4	3.7	2.4	0.0	–	–	–	–

As can be seen from Table 4, the elasticity moduli and computation time will be rising gradually with the increase of the selected total number of modes of the fiber-reinforced composite, and the theoretical calculation errors also show an increasing trend. Therefore, it is not true that when the higher number of mode is considered, the more accurate identification results of mechanical parameters will be obtained, and how close the theoretical FRF can approach the experimental FRF is what counts most of all. Therefore, in this example, the experimental FRF in the first three modes is the better choice to construct frequency relative error function and admittance relative error function in FRFAM.

Identification results under different approximation points

In this section, in order to discuss the influence of the number of approximation points on identification accuracy and identification efficiency of mechanical parameters of fiber-reinforced composite, 3, 5, 7, and 9 approximation points are selected in the MATLAB program with the unchanged calculation step size of 1% and free boundary condition. Then, when the frequency relative error function and admittance relative error function in the first three modes get the minimum value, the optimum mechanical parameters can be identified, as seen in Table 5, and the corresponding time expenditures are also listed in Table 5.

Table 5.

The identified mechanical parameters of fiber-reinforced composite by FRFAM under different approximation points and time expenditures.

Approximation points	$E_1$ (GPa)	$E_2$ (GPa)	$G_{12}$ (GPa)	${\nu }_{12}$	${\eta }_1$	${\eta }_2$	${\eta }_{12}$	Calculation time (s)
3	107.8	7.7	4.3	0.3	0.0018	0.0057	0.0267	185.4
5	107.8	7.7	4.3	0.3	0.0015	0.0054	0.0264	285.4
7	107.8	7.7	4.3	0.3	0.0015	0.0054	0.0264	569.5
9	107.8	7.7	4.3	0.3	0.0015	0.0054	0.0264	886.9

As can be seen from Table 5, the identification results of elasticity moduli $E_1$ and $E_2$ , shear modulus $G_{12}$ , and Poisson ratio ${\nu }_{12}$ do not change under different approximation points, so the number of approximation points will not affect the identification accuracy of the elastic modulus and Poisson’s ratio. However, as the number of approximation points increases, the identified loss factors will be more accurate, yet the time expenditures in the calculation will also increase correspondingly. Therefore, it is necessary to choose the appropriate number of approximation points. In this example, five to seven approximation points are accurate enough to finish the FRF approximation work.

Identification results under different calculation step sizes

In this section, in order to evaluate FRFAM objectively, the influence of different calculation step sizes on identification accuracy and identification efficiency of mechanical parameters of the fiber-reinforced composite is discussed. First, choose the FRF data under free boundary condition in the first three modes and set five approximation points in the MATLAB program. Then, select the step size of 20%, 7.5%, 3.5%, 1%, 0.5%, and 0.1% to identify the concerned mechanical parameters along the longitudinal, transverse, and shear direction of the fiber composites. The identification results under different calculation step sizes and their time expenditures can be seen in Table 6, from which it can be discovered that calculation step sizes will have an important impact on the identification accuracy and efficiency. With the decrease of calculation step size, the theoretical calculation errors show a decreasing trend compared with the ones provided by the manufacturer, but the time expenditures in the calculation will increase sharply. Therefore, it is necessary to choose appropriate calculation step sizes. In this example, the step size of 1% is accurate enough to finish the FRF approximation work and its calculation time is less than 300 s, which is within an acceptable range. So, it can be used as the reference in the similar FRF approximation calculation program for the identification of mechanical parameters of other fiber composites.

Table 6.

The identified mechanical parameters of fiber-reinforced composite by FRFAM under different step sizes and their time expenditures.

Step size	$E_1$ (GPa)	$E_2$ (GPa)	$G_{12}$ (GPa)	${\nu }_{12}$	${\eta }_1$	${\eta }_2$	${\eta }_{12}$	Calculation time (s)
20%	109.2	7.7	4.3	0.3	0.0040	0.0080	0.030	50.4
7.5%	108.2	7.7	4.3	0.3	0.0023	0.0068	0.027	95.7
3.5%	107.9	7.7	4.3	0.3	0.0021	0.0063	0.027	165.9
1%	107.8	7.7	4.3	0.3	0.0015	0.0054	0.0264	285.4
0.5%	107.8	7.7	4.3	0.3	0.0015	0.0054	0.0264	412.9
0.1%	107.8	7.7	4.3	0.3	0.0015	0.0054	0.0264	2005.2

Conclusion

In this article, FRFAM is proposed to identify the mechanical parameters of fiber-reinforced composites, and its correctness and practicability have been validated in the verification example and the case study. It has been found that the boundary conditions, approximation points, total number of modes, and calculation step size have an important impact on the identification accuracy and efficiency. Because the identified errors of elastic moduli and Poisson’s ratio under the cantilever boundary condition are larger than the ones under the free boundary condition, it is suggested that the experimental FRF data should be obtained under free boundary condition. Meanwhile, after comparing with the identified results, the step size of 1% and five to seven approximation points would be accurate enough to finish the FRF approximation work. Finally, it is not true that when the higher number of mode is considered, the more accurate identification results of mechanical parameters will be obtained. The identification accuracy of FRFAM depends on how close the theoretical FRF can approach the experimental FRF, and usually the experimental FRF in the first three modes is the better choice to establish the frequency relative error function and admittance relative error function in FRFAM.