Techniques for Identification of Bending and Extensional Elastic Stiffness Matrices on Thin Composite Material Plates Based on Virtual Field Method (VFM): Theoretical and Numerical Aspects

1. Introduction

Nowadays, the demand and the necessity for the use of materials with specificcharacteristics have increased in many engineering fields. Due to this necessity of making new materials, composite materials have been an alternative, or maybe the unique option, to attempt a large number of design requirements suchas high strength-to-weight ratio,high resistance to mechanical shocks, chemical attacks, corrosion, and fatigue, that cannot be obtained only from the commonly used structural materials (metals, ceramics, polymers and wood). Because of this, their applications are present in the main industries such as aerospace, automotive, marine, and sportive

Thanks to their flexibility characteristicsthere are many combinations and arrangements and, consequently, constitutive properties, that are possible to be achieved. This particularity represents one of the main advantages of these materials. Nevertheless, some factors related to arrangements as the number of layers and the orientation of the fibers can introduce abehavior called anisotropy that, in the most of the cases, is not required. The anisotropy makes the structural analysis more complex due to increasing the number of independent variables, as for example, the number of bending and extensional elastic stiffness constants.

Recently, wide part of works presented by scientific literature whose goal is to identify constitutive parameters of materials (these being composites or not) is based in the called “inverse problems”. Experimental data such as geometry, resultant forces and strain (or stress) fields are used as input data, and the unknown variables are the required constitutive parameters. In general, the solution is basically associated to two methods: iterative (or also called indirect methods) and non-iterative (or also called of direct methods). The first one is related to optimization problems where the design variables are constitutive parameters and the objective function represents, in general, a residue (or error) between experimental and numerical (generally obtained by finite elements) data. For the numerical simulations, it is considered structures that have the same geometrical characteristics and boundary conditions of the real ones. For each step, required parameters are checked out, and, the optimum represents the iteration whose residue (objective function) has its lowest value. Unlike indirect methods, the direct methods are ones where the required parameters are computed from the solution of constitutive equation(s) that are functions of these parameters.

In general, according to the type of experimental test, it is possible to separate the methods of elastic property identification in two categories: static (destructives and non-destructives) and dynamic methods (non-destructive), as shown in Fig.1. A large number of identification techniques that use data from these categories of tests have been proposed, especially ones dedicated to composite materials. It is possible to say that these techniques identify effective properties of the entire material. The way as each formulation is built, and, the adopted procedures and devices are the main differences among the many proposed methodologies.

Static tests with monotonic load are experimental tests that were more commonly used in the last years, and maybe the simplestones, for this material property identification. Despite the simplicity of these tests, some aspects render them less attractive than dynamic tests, such as the fact of requiring a number of samples with fiber orientations according to standard norms, for example, American Society for Testing and Materials (ASTM), which, in the most of cases, are not in accordance to the real characteristics of the required material. Furthermore, some variables difficult of controlling during the tests can contribute to worsen the experimental results, e.g., the presence of non-uniform stress fields near the ends of the sample from the clamped boundary conditions. For these reasons, dynamic tests have been considered an interesting alternative. In general, they are tests that combine experimental data with numerical methods, and allow the identification of elastic constants from only one unique sample or even from composite material part. Sample is usually thin plate (that reflect Kirchhoff’s hypotheses), cylindrical shell, or beam. In many cases, input data of the numerical methods are natural frequencies and/or mode shapes.

Figure 1.

Methodologies more used to identify elastic properties of materials.

Many authors have proposed to identifythe elastic constants by iterative procedures adopting Rayleigh-Ritz and finite element as numerical methods. The difference among the identification techniques based on these iterative procedures is basically in the way as the optimization problem is formulated, for example, the type of adopted search method to find the minimal, the boundary conditions, the geometric characteristics of the samples, the type of anisotropy of the test material, the type of experimental devices, and the numerical method used to compute the mode shapes (or operational modes) with their respective frequencies (Deobald& Gibson, 1988; Pedersen &Frederiksen, 1992; Lai & Lau, 1993; Ayorinde& Gibson, 1995; Rikards&Chate, 1998; Ayorinde& Yu, 1999, 2005; Rikards et al., 1999; Bledzki et al., 1999; Hwang & Chang, 2000; Araujo et al., 2000; Chakraborty&Mukhopadhyay, 2000; Rikards et al., 2001; Lauwagie et al., 2003; Lauwagie et al., 2004; Lee &Kam, 2006; Cugnoni et al., 2007; Bruno et al., 2008; Pagnotta&Stigliano, 2008; Diveyev&Butiter, 2008a, 2008b).

In works that do not use iterative process, natural frequencies and mode shapes, or operational frequencies and modes, are input data of an algorithm based on the differential equation that governs the transversal vibration of sample in a specific direction and under specific boundary conditions (Gibson, 2000; Alfano&Pagnotta, 2007). In this methodology, it can be included the use of Virtual Fields Method, VFM (Grédia, 1996, 2004; Grédiac& Paris, 1996; Grédiac&Pierron, 1998, 2006; Pierron et al., 2000, 2007; Pierron&Grédiac, 2000; Grédiac et al., 1999a, 1999b, 2001, 2006; Giraudeau&Pierron, 2003, 2005; Chalal et al., 2006; Toussaint et al., 2006; Avril&Pierron, 2007; Pierron et al., 2007; Avril et al., 2008; Giraudeau et al., 2006). For VFM, weighting functions are the called virtual fields. Due to sensitivity to experimental errors and the presence of noise during the dynamic testing, it was proposed the use of specific virtual fields named “special virtual fields” (Grédiac et al., 2002a, 2002b, 2003). In order to decrease the noise contribution, the use of more accurate experimental modal analysis techniques or the application of some signal smoothing (or filtering) technique is mandatory.

The majority of works identifies only the bending stiffness matrix or directly the engineering elastic constants. However, the extensional elastic stiffness matrix is also needed to model composite materials under multi-axial loads. In general, these stiffness matrices are independent. The extensional stiffness matrix relates the in-plane resultant forces to the midplane strains, and, the bending stiffness matrix relates the resultant moments to the plate curvatures. In a laminate composite, if only the stacking sequence of layers is changed, the bending matrix is changed but the extensional matrix remains the same. In other words, different laminates can have different bending stiffness matrices and the same extensional stiffness matrix. It is not possible to obtain the extensional matrix from the bending matrix. It will be possible only if the stacking sequence of layers and their thickness are known and, also, if the material is the same for all lamina.

Sometimes, it is more convenient to use effective laminate engineering constants rather than the laminate stiffness. These effective laminated engineering constants may be easily obtained from the extensional elastic constants. However, due to difficulties on experimental in-plane modal analysis, such as the necessity of using specific devices to measure in-plane displacements and to excite high frequencies, the identification of extensional elastic stiffness constants using modal testing is less attractive. The main challenge to perform in-plane vibration testing is the excitation and measurement ofonly in-plane and not out-of-plane vibration modes. Today there are some new techniques that are suitable for this kind of problems, for example, the excitation by piezoelectric (PZT) and measurements by digital image correlation.

In this chapter, a review about the VFM applied to compute bending elastic stiffness constants proposed by Grédiac& Paris, 1996 is presented. Furthermore, a formulation based on the VFM is proposed in order to identify the extensional elastic stiffness matrix ofKirchhoff’s thin plates. The linear system of equations that provides the required elastic constants is obtained from differential equations that govern the forced vibration of anisotropic, symmetric and non-damped plates under in-plane loads. The common procedures to find the weak form (or integral form) of these equations are applied here. The correct choice of weighting functions (which are the virtual fields) and mode shapes representsa key characteristic to the accuracy of the results. Numerical simulations using anisotropic, orthotropic, quasi-isotropic plates are carried out to demonstratethe accuracy of the methodology.

2. Identification of elastic constants using VFM

2.1. Review of the Virtual Fields Method - VFM

The VFM has been developed for extracting constitutive parameters from full-field measurements and it is associated to problems of identification of parameters from constitutive equations. Two cases are clearly distinguished: constitutive equations depending linearly on the constitutive parameters and non-linear constitutive equations. The type of constitutive equations is chosen a priori for its relevancy and objective of it is to determine the parameters which govern the constitutive equations. The main difficulty comes from the fact that the measured displacement or strain components are generally not directly related to the unknown parameters (Grédiac et al., 2006), and no closed-form solution for the displacement, strain and stress fields is available.

Mathematically, the VFM is based on the principle of virtual work and can be written as:

$$-{\displaystyle {\int }_V\sigma :{\epsilon }^{*}dV}+{\displaystyle {\int }_{S_f}T.u^{*}dS}+{\displaystyle {\int }_{V}f.u^{*}dV}={\displaystyle {\int }_V\rho \gamma udV}$$ E1

where V is volume of the solid, σ is the actual stress tensor, ε ^* is the virtual strain tensor, T is the distribution vector of loading tractions acting on the boundary, S _f is the part of the solid boundary where the tractions are applied, u^* is the virtual displacement vector, f is the distribution of volume forces acting on V, ρ is the density and γ the acceleration. Eq. (1) is verified for any kinematically admissible virtual field (u^*, ε ^*). Kinematically admissible means that u^* must be continuous across the whole volume and it must be equal to the prescribed displacement on the boundary S _f where displacements are prescribed. Let’s introduce the constitutive equations in the general case as:

$$\sigma =g\left(\epsilon \right)$$ E2

whereg is a function of the actual strain. Thus, when constitutive equations are introduced and volume forces are disregarded, Eq. (1) can be rewritten as:

$$-{\displaystyle {\int }_{V}g\left(\epsilon \right):{\epsilon }^{*}dV}+{\displaystyle {\int }_{S_f}T.u^{*}dS}={\displaystyle {\int }_V\rho \gamma u^{*}dV}$$ E3

It is possible to see in Eq. (3) that each virtual field originates a new equation involving the constitutive parameters. The VFM relies on this important property. It is a method based on setting virtual fields that provide a set of equations. This set of equations is used to extract the required unknown constitutive parameters. The correct choice of the virtual fields that combine to actual fields in Eq. (3) is the key issue of the method. Their number and their type depend on the nature of g in Eq. (3).

2.2. Review of the identification method of bending elastic stiffness matrix

The method proposed by Grédiac&Paris,1996, consists of obtaining elastic constants based on the partial differential equation that governs the transversal vibration of an anisotropic thin plate (Kirchhoff’s plate). This equation is given by:

$$D_{11}\frac{{\partial }^{4}w}{\partial x^4}+4D_{16}\frac{{\partial }^{4}w}{\partial x^3\partial y}+2\left(D_{12}+2D_{66}\right)\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+4D_{26}\frac{{\partial }^{4}w}{\partial x\partial y^3}+D_{22}\frac{{\partial }^{4}w}{\partial y^4}=-\rho h\frac{{\partial }^{2}w}{\partial t^2}$$ E4

where D _ij are thin plate bending stiffness constants; ρis the mass density of the material; h is the plate thickness; x and y are coordinates of the plate; t is time; and w(x,y,t) is the deflection function that represents the transversal displacement of a point of the plate at an instant t. Eq. (4) doesn’t state the global equilibrium of the plate since the excitation force and damping are not considered. However, for many composite materials, as for example, aeronautic carbon epoxy tested in this work, the damping is low enough to disregard its contribution in the formulation. Besides, if the input data refer to resonant response of the plate, the work provided by the excitation is balanced by internal dissipation of the plate. A detailed discussion about when excitation and damping should be considered in Eq. (4) can be found in Giraudeau&Pierron, 2006.

After some mathematical manipulations in Eq. (4), Grédiac& Paris,1996 obtained a linear system in which the unknown variables are the elastic constants. Briefly, the sequence of operations is as follows: (a) multiply both sides of Eq. (4) by an arbitrary weighting function; (b) integrate twice by parts along the plate domain; (c) eliminate the boundary integrals by applying the free-edge boundary conditions; (d) decompose the displacement function w(x, y, t) as a product of the deflection amplitude Φ and sin(ѡt), where w is the out-of-plane natural frequency of a particular mode shape of the plate; and (e) choose appropriate weighting functions and mode shapes to build the matrix of the linear system. At this point, as Grédiac& Paris, 1996 explain, the choice of mode shapes associated with the weighting function is extremely important for the accuracy of this method. Three particular modes are strongly dependent on the required coefficients D _ij : a twisting mode that strongly depends on terms D ₆₆,D ₁₆, and D ₂₆; a bending mode along direction 1 that strongly depends on terms D ₁₁, D ₁₂, and D ₁₆, and a bending mode along direction 2 that strongly depends on D ₂₂, D ₁₂, and D ₂₆. These modes present smooth curvatures and are generally among the first modes, with lower frequencies. If these modes are not found, it is recommended to use modes that have similar shapes to them. Furthermore, they are modes that can be approximated by quadratic functions with constant curvatures: x ², y ², and xy. For this reason, these quadratic functions were the weighting functions chosen by Grédiac& Paris, 1996. Thus, using these previous quadratic-weighting functions, the following simplified system of equations can be obtained:

$$\left[\begin{array}{cccccc}\cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ K_{xx}^{\left(j\right)}& 0& K_{yy}^{\left(j\right)}& 0& 2K_{xy}^{\left(j\right)}& 0\\ 0& K_{yy}^{\left(j\right)}& K_{xx}^{\left(j\right)}& 0& 0& 2K_{xy}^{\left(j\right)}\\ 0& 0& 0& 2K_{xy}^{\left(j\right)}& K_{xx}^{\left(j\right)}& K_{yy}^{\left(j\right)}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ K_{xx}^{\left(k\right)}& 0& K_{yy}^{\left(k\right)}& 0& 2K_{xy}^{\left(k\right)}& 0\\ 0& K_{yy}^{\left(k\right)}& K_{xx}^{\left(k\right)}& 0& 0& 2K_{xy}^{\left(k\right)}\\ 0& 0& 0& 2K_{xy}^{\left(k\right)}& K_{xx}^{\left(k\right)}& K_{yy}^{\left(k\right)}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \end{array}\right]\left\{\begin{array}{c}{D}_{11}\\ D_{22}\\ D_{12}\\ D_{66}\\ D_{16}\\ D_{26}\end{array}\right\}=-\frac{\rho h}{2}\left\{\begin{array}{c}\cdots \\ {\omega }_j^2{\displaystyle {\int }_S{\Phi }^{\left(j\right)}{x}^{2}dS}\\ {\omega }_j^2{\displaystyle {\int }_S{\Phi }^{\left(j\right)}{y}^{2}dS}\\ {\omega }_j^2{\displaystyle {\int }_S{\Phi }^{\left(j\right)}xydS}\\ \cdots \\ {\omega }_k^2{\displaystyle {\int }_S{\Phi }^{\left(k\right)}{x}^{2}dS}\\ {\omega }_k^2{\displaystyle {\int }_S{\Phi }^{\left(k\right)}{y}^{2}dS}\\ {\omega }_k^2{\displaystyle {\int }_S{\Phi }^{\left(k\right)}xydS}\\ \cdots \end{array}\right\}$$ E5

where indices j and k represent a specific mode shape; and S is the plate domain. Elements of the matrix of Eq. (5) are given by:

$$K_{xx}={\displaystyle \underset{S}{\int }\frac{{\partial }^2\Phi \left(x,y\right)}{\partial x^2}}dS,K_{yy}={\displaystyle \underset{S}{\int }\frac{{\partial }^2\Phi \left(x,y\right)}{\partial y^2}}dS,K_{xy}={\displaystyle \underset{S}{\int }\frac{{\partial }^2\Phi \left(x,y\right)}{\partial x\partial y}}dS$$ E6

Eq. (5) can be represented in matrix form as:

$$KD=C$$ E7

where, considering L as the number of modes used in the linear system of equations, K is a 3L x 6 matrix, D is a 6 x 1 matrix, and C is a 3L x 1 matrix. As can be seen, Eq. (7) is an overdetermined system of equations. The solution can be found by least squares:

$$D={\left(K^{T}K\right)}^{-\text{1}}\left(K^{T}C\right)$$ E8

2.3. Identification method of the extensional elastic stiffness matrix

In the general case of composite laminates, each lamina is assumed to have orthotropic material properties. After the assembly, the behavior can be anisotropic due to the interaction of different laminas. Considering a plate under plane state of stress and using Hooke’s generalized law, stresses can be integrated over its thickness yielding the following force-deformation equations:

$$\left\{\begin{array}{c}N\\ \cdots \\ M\end{array}\right\}=\left[\begin{array}{ccc}A& ⋮& B\\ \cdots & ⋮& \cdots \\ B& ⋮& D\end{array}\right]\left\{\begin{array}{c}\epsilon \\ \cdots \\ \kappa \end{array}\right\}$$ E9

where N and M are vectors that contain normal forces and resultant moments, respectively, A is the extensional elastic stiffness matrix, B is the coupling elastic stiffness matrix (B is a null matrix in the case of a symmetric laminate), D is the bending elastic stiffness matrix, ε and κ are vectors that contain middle plane linear strains and rotations, respectively. Considering a symmetrical (B = [0]) and fully anisotropic laminate under free-edge in-plane vibration (the plate is not under bending) and using the equilibrium relations, the following equations can be written:

E10

E11

where A _ij are the elements of matrix A (i,j = 1,2,6); ρ is the mass density; h is the plate thickness; x and y are the coordinates in the plate plane; t is the time, and u(x,y,t) and v(x,y,t) are functions that represent the displacements along x and y direction, respectively, of a point with coordinates (x,y) of the plate at an instant t. Multiplying Eq. (10) by a weighting function W(x, y) and integrating along the domain of the plate, we can obtain:

E12

whereΩ is the plate domain. Using chain rule, i.e.,

$$\frac{d}{dx}\left(f.g\right)=f\frac{dg}{dx}+g\frac{df}{dx}$$ E13

wheref and g are any two continuous functions, and Green’s theorem, i.e.,

$${\displaystyle {\int }_{\Omega }\frac{\partial f}{\partial x}}d\Omega ={\displaystyle {\int }_{\Gamma }f{n}_{x}dx}$$ E14

,

where Γ is the boundary domain, and n _x is the component of the normal unity vector in directions x, the left hand side terms of Eq. (12) can be written as:

E15

wheren _y is the component of the normal unity vector in direction y. From the constitutive equation, one has:

$$\left\{\begin{array}{c}{N}_x\\ N_y\\ N_{xy}\end{array}\right\}=\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{16}\\ A_{12}& A_{22}& A_{26}\\ A_{16}& A_{26}& A_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right\}$$ E16

whereN _x and N _y are axial forces per unit length along directions x and y, respectively, N _xy is the shear force per unit length along plane xy, ε _x and ε _y are the middle surface axial strain along directions x and y, respectively, and γ _xy is the shear angular strain along plane xy. Rewriting Eq. (16), one obtains:

$$N_x=A_{11}{\epsilon }_x+A_{12}{\epsilon }_y+A_{16}{\gamma }_{xy}=A_{11}\frac{\partial u}{\partial x}+A_{12}\frac{\partial v}{\partial y}+A_{16}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)$$ E17

,

$$N_y=A_{12}{\epsilon }_x+A_{22}{\epsilon }_y+A_{26}{\gamma }_{xy}=A_{12}\frac{\partial u}{\partial x}+A_{22}\frac{\partial v}{\partial y}+A_{26}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)$$ E18

,

$$N_{xy}=A_{16}{\epsilon }_x+A_{26}{\epsilon }_y+A_{66}{\gamma }_{xy}=A_{16}\frac{\partial u}{\partial x}+A_{26}\frac{\partial v}{\partial y}+A_{66}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)$$ E19

.

Multiplying Eq. (17) by a weighting function W and by n _x , which is the x component of unit normal vector n, and integrating on the boundary Γ, one obtains:

$$A_{11}{\displaystyle {\int }_{\Gamma }\frac{\partial u}{\partial x}W{n}_{x}d\Gamma }+A_{12}{\displaystyle {\int }_{\Gamma }\frac{\partial v}{\partial y}W{n}_{x}d\Gamma }+A_{16}\left({\displaystyle {\int }_{\Gamma }\frac{\partial u}{\partial y}}W{n}_{x}d\Gamma +{\displaystyle {\int }_{\Gamma }\frac{\partial v}{\partial x}W{n}_{x}d\Gamma }\right)={\displaystyle {\int }_{\Gamma }{N}_{x}W{n}_{x}d\Gamma }$$ E20

.

Multiplying Eq. (19) by a weighting function W and by n _y , which is the y component of the unit normal vector n, and integrating on the boundary Γ, yields:

$$A_{16}{\displaystyle {\int }_{\Gamma }\frac{\partial u}{\partial x}W{n}_{y}d\Gamma }+A_{26}{\displaystyle {\int }_{\Gamma }\frac{\partial v}{\partial y}W{n}_{y}d\Gamma }+A_{66}\left({\displaystyle {\int }_{\Gamma }\frac{\partial u}{\partial y}}W{n}_{y}d\Gamma +{\displaystyle {\int }_{\Gamma }\frac{\partial v}{\partial x}W{n}_{y}d\Gamma }\right)={\displaystyle {\int }_{\Gamma }{N}_{xy}W{n}_{y}d\Gamma }$$ E21

.

Substituting Eqs. (20) and (21) into Eq. (15), and reorganizing the terms, one can write:

$${\displaystyle \underset{T}{\overset{}{\int }}{N}_{x}W{n}_{x}dT+}{\displaystyle \underset{T}{\overset{}{\int }}{N}_{xy}W{n}_{y}dT-{\displaystyle \underset{\Omega }{\overset{}{\int }}\left[\begin{array}{l}{A}_{11}\left(\frac{\partial u}{\partial x}\frac{\partial W}{\partial x}\right)+A_{16}\left(\frac{\partial u}{\partial y}\frac{\partial W}{\partial x}+\frac{\partial u}{\partial x}\frac{\partial W}{\partial y}+\frac{\partial v}{\partial x}\frac{\partial W}{\partial x}\right)+\\ A_{12}\left(\frac{\partial u}{\partial y}\frac{\partial W}{\partial x}\right)+A_{66}\left(\frac{\partial u}{\partial y}\frac{\partial W}{\partial y}+\frac{\partial v}{\partial x}\frac{\partial W}{\partial y}\right)+A_{26}\left(\frac{\partial v}{\partial y}\frac{\partial W}{\partial y}\right)\end{array}\right]}d\Omega }\text{=}\rho h{\displaystyle \underset{\Omega }{\overset{}{\int }}\frac{{\partial }^{2}u}{\partial t^2}}Wd\Omega$$ E22

Now, if free-edge boundary conditions are considered, boundary integrals of Eq. (22) vanish. Considering that the plate is vibrating, functions u and v can be written as:

$$u\left(x,y,t\right)=U\left(x,y\right)\sin \left(\varpi t\right)$$ E23

,

$$v\left(x,y,t\right)=V\left(x,y\right)\sin \left(\varpi t\right)$$ E24

,

where $$\varpi$$ is the in-plane natural frequency associated to any in-plane mode, and U(x,y) and V(x,y) are the displacement amplitudes along directions x and y, respectively, of a point with coordinate (x,y). In this sense, the amplitude is only a function of coordinates x and y. Eliminating the boundary integrals of Eq. (22), substituting Eqs. (23) and (24) into Eq. (22), and considering any mode j, one obtains:

$$\begin{array}{l}{\displaystyle {\int }_{\Omega }[A_{11}\left(\frac{\partial U^{\left(j\right)}}{\partial x}\frac{\partial W}{\partial x}\right)+A_{16}\left(\frac{\partial U^{\left(j\right)}}{\partial y}\frac{\partial W}{\partial x}+\frac{\partial U^{\left(j\right)}}{\partial x}\frac{\partial W}{\partial y}+\frac{\partial V^{\left(j\right)}}{\partial x}\frac{\partial W}{\partial x}\right)}+A_{12}\left(\frac{\partial V^{\left(j\right)}}{\partial y}\frac{\partial W}{\partial x}\right)+\\ +A_{66}\left(\frac{\partial U^{\left(j\right)}}{\partial y}\frac{\partial W}{\partial y}+\frac{\partial V^{\left(j\right)}}{\partial x}\frac{\partial W}{\partial y}\right)+A_{26}\left(\frac{\partial V^{\left(j\right)}}{\partial y}\frac{\partial W}{\partial y}\right)]d\Omega =\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{U}^{\left(j\right)}Wd\Omega }\end{array}$$ E25

Now, if the same previous mathematical procedures used in Eq. (10) are used in Eq. (11), one obtains:

$$\begin{array}{l}{\displaystyle {\int }_{\Omega }[A_{16}\left(\frac{\partial U^{\left(j\right)}}{\partial x}\frac{\partial W}{\partial x}\right)+A_{12}\left(\frac{\partial U^{\left(j\right)}}{\partial x}\frac{\partial W}{\partial y}\right)}+A_{66}\left(\frac{\partial U^{\left(j\right)}}{\partial y}\frac{\partial W}{\partial x}+\frac{\partial V^{\left(j\right)}}{\partial x}\frac{\partial W}{\partial x}\right)+A_{26}\left(\frac{\partial U^{\left(j\right)}}{\partial y}\frac{\partial W}{\partial y}\right)+\\ +A_{26}\left(\frac{\partial V^{\left(j\right)}}{\partial y}\frac{\partial W}{\partial x}+\frac{\partial V^{\left(j\right)}}{\partial x}\frac{\partial W}{\partial y}\right)+A_{22}\left(\frac{\partial V^{\left(j\right)}}{\partial y}\frac{\partial W}{\partial y}\right)]d\Omega =\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{V}^{\left(j\right)}Wd\Omega }\end{array}$$ E26

2.3.1. Choice of the weighting functions

Eqs. (25) and (26) are theoretically valid for isotropic, orthotropic, or anisotropic plates, provided that the laminate is symmetrical. As it can be seen, the function W is arbitrary since itself and its first order derivative is continuous in the domain Ω. Amplitudes U(x,y) and V(x,y), and frequencies $${\varpi }_j$$ are obtained from dynamic tests. Dimensions of the plate and parameters ρ and h can also be easily measured on the sample plate. Thus, the next steps are the choice of a suitable numerical method to compute the derivatives and integrals in Eqs. (25) and (26). Furthermore, suitable mode shapes and weighting functions should be chosen. In this work, finite differences and Gauss-Legendre numerical integration scheme are used to compute these derivatives and integrals, respectively. For numerical reasons, modes with several sign changes in the mode shape are avoided because their numerical derivatives and integrals are more sensitive to errors (Grédiac& Paris, 1996). Generally, first modes present more smooth curvatures and are, at the same time, easier to be obtained experimentally. It is worth noting that in-plane modal analysis presents much higher frequencies than transverse modal analysis (bending modes). This is because stiffness along direction x and y is much higher than stiffness along the transversal direction of the plate. For numerical reasons, smooth mode shapes associated with smooth weighting functions are preferred. For all these reasons, and in order to simplify Eqs. (25) and (26), the following group of weighting functions are proposed:

W(x, y) = x², which applied to Eqs. (25) and (26) provides the following integral equations:

$$2{\displaystyle {\int }_{\Omega }\left[A_{11}\left(\frac{\partial U^{\left(j\right)}}{\partial x}x\right)+A_{16}\left(\frac{\partial U^{\left(j\right)}}{\partial y}x+\frac{\partial V^{\left(j\right)}}{\partial x}x\right)+A_{12}\left(\frac{\partial V^{\left(j\right)}}{\partial y}x\right)\right]}d\Omega =\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{U}^{\left(j\right)}{x}^{2}d\Omega }$$ E27

,

$$2{\displaystyle {\int }_{\Omega }\left[A_{16}\left(\frac{\partial U^{\left(j\right)}}{\partial x}x\right)+A_{66}\left(\frac{\partial U^{\left(j\right)}}{\partial y}x+\frac{\partial V^{\left(j\right)}}{\partial x}x\right)+A_{26}\left(\frac{\partial V^{\left(j\right)}}{\partial y}x\right)\right]d\Omega =\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{V}^{\left(j\right)}{x}^{2}d\Omega }}$$ E28

.

W(x, y) = y², which applied to Eqs. (25) and (26) provides the following integral equations:

$$2{\displaystyle {\int }_{\Omega }\left[A_{16}\left(\frac{\partial U^{\left(j\right)}}{\partial x}y\right)+A_{66}\left(\frac{\partial U^{\left(j\right)}}{\partial y}y+\frac{\partial V^{\left(j\right)}}{\partial x}y\right)+A_{26}\left(\frac{\partial V^{\left(j\right)}}{\partial y}y\right)\right]d\Omega =\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{U}^{\left(j\right)}{y}^{2}d\Omega }}$$ E29

,

$$2{\displaystyle {\int }_{\Omega }\left[A_{12}\left(\frac{\partial U^{\left(j\right)}}{\partial x}y\right)+A_{26}\left(\frac{\partial U^{\left(j\right)}}{\partial y}y+\frac{\partial V^{\left(j\right)}}{\partial y}y\right)+A_{22}\left(\frac{\partial V^{\left(j\right)}}{\partial y}y\right)\right]d\Omega =\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{V}^{\left(j\right)}{y}^{2}d\Omega }}$$ E30

.

W(x, y) = 2xy, which applied to Eqs. (25) and (26) provides the following integral equations:

$$\begin{array}{l}2{\displaystyle {\int }_{\Omega }[A_{11}\left(\frac{\partial U^{\left(j\right)}}{\partial x}y\right)+A_{16}\left(\frac{\partial U^{\left(j\right)}}{\partial y}y+\frac{\partial U^{\left(j\right)}}{\partial x}x+\frac{\partial V^{\left(j\right)}}{\partial x}y\right)+A_{12}\left(\frac{\partial V^{\left(j\right)}}{\partial y}y\right)+}\\ +A_{66}\left(\frac{\partial U^{\left(j\right)}}{\partial y}x+\frac{\partial V^{\left(j\right)}}{\partial x}x\right)+A_{26}\left(\frac{\partial V^{\left(j\right)}}{\partial y}x\right)]d\Omega =\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{U}^{\left(j\right)}\left(2xy\right)d\Omega }\end{array}$$ E31

,

$$\begin{array}{l}2{\displaystyle {\int }_{\Omega }[A_{16}\left(\frac{\partial U^{\left(j\right)}}{\partial x}y\right)+A_{66}\left(\frac{\partial U^{\left(j\right)}}{\partial y}y+\frac{\partial V^{\left(j\right)}}{\partial x}y\right)+A_{12}\left(\frac{\partial U^{\left(j\right)}}{\partial x}x\right)+}\\ +A_{26}\left(\frac{\partial U^{\left(j\right)}}{\partial y}x+\frac{\partial V^{\left(j\right)}}{\partial y}y+\frac{\partial V^{\left(j\right)}}{\partial x}x\right)+A_{22}\left(\frac{\partial V^{\left(j\right)}}{\partial y}x\right)]d\Omega =\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{V}^{\left(j\right)}\left(2xy\right)d\Omega }\end{array}$$ E32

Defining:

$${\displaystyle {\int }_{\Omega }\left(\frac{\partial U}{\partial x}x\right)d\Omega }=K_{uxx}$$ , $${\displaystyle {\int }_{\Omega }\left(\frac{\partial U}{\partial y}x\right)}d\Omega =K_{uyx}$$ , $${\displaystyle {\int }_{\Omega }\left(\frac{\partial U}{\partial x}y\right)}d\Omega =K_{uxy}$$ , $${\displaystyle {\int }_{\Omega }\left(\frac{\partial U}{\partial y}y\right)}d\Omega =K_{uyy}$$ , $${\displaystyle {\int }_{\Omega }\left(\frac{\partial V}{\partial x}x\right)d\Omega }=K_{vxx}$$ , $${\displaystyle {\int }_{\Omega }\left(\frac{\partial V}{\partial y}x\right)}d\Omega =K_{vyx}$$ , $${\displaystyle {\int }_{\Omega }\left(\frac{\partial V}{\partial x}y\right)}d\Omega =K_{vxy}$$ , $${\displaystyle {\int }_{\Omega }\left(\frac{\partial V}{\partial y}y\right)}d\Omega =K_{vyy}$$ .

Eqs. (34)-(38) can be written as:

$$A_{11}{K}_{uxx}^{\left(j\right)}+A_{16}\left(K_{uyx}^{\left(j\right)}+K_{vxx}^{\left(j\right)}\right)+A_{12}{K}_{vyx}^{\left(j\right)}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{U}^{\left(j\right)}{x}^{2}d\Omega }$$ E33

,

$$A_{16}{K}_{uxx}^{\left(j\right)}+A_{66}\left(K_{uyx}^{\left(j\right)}+K_{vxx}^{\left(j\right)}\right)+A_{26}{K}_{vyx}^{\left(j\right)}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{V}^{\left(j\right)}{x}^{2}d\Omega }$$ E34

,

$$A_{16}{K}_{uxy}^{\left(j\right)}+A_{66}\left(K_{uyy}^{\left(j\right)}+K_{vxy}^{\left(j\right)}\right)+A_{26}{K}_{vyy}^{\left(j\right)}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{U}^{\left(j\right)}{y}^{2}d\Omega }$$ E35

,

$$A_{12}{K}_{uxy}^{\left(j\right)}+A_{26}\left(K_{uyy}^{\left(j\right)}+K_{vxy}^{\left(j\right)}\right)+A_{22}{K}_{vyy}^{\left(j\right)}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{V}^{\left(j\right)}{y}^{2}d\Omega }$$ E36

,

$$\begin{array}{l}{A}_{11}{K}_{uxy}^{\left(j\right)}+A_{16}\left(K_{uyy}^{\left(j\right)}+K_{uxx}^{\left(j\right)}+K_{vxy}^{\left(j\right)}\right)+A_{12}{K}_{vyy}^{\left(j\right)}+\\ +A_{66}\left(K_{uyx}^{\left(j\right)}+K_{vxx}^{\left(j\right)}\right)+A_{26}{K}_{vyx}^{\left(j\right)}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{U}^{\left(j\right)}\left(2xy\right)d\Omega }\end{array}$$ E37

,

$$\begin{array}{l}{A}_{16}{K}_{uxy}^{\left(j\right)}+A_{66}\left(K_{uyy}^{\left(j\right)}+K_{vxy}^{\left(j\right)}\right)+A_{12}{K}_{uxx}^{\left(j\right)}+A_{26}\left(K_{uyx}^{\left(j\right)}+K_{vyy}^{\left(j\right)}+K_{vxx}^{\left(j\right)}\right)+\\ +A_{22}{K}_{vyx}^{\left(j\right)}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\rho h{\varpi }_j^2{\displaystyle {\int }_{\Omega }{V}^{\left(j\right)}\left(2xy\right)d\Omega }\end{array}$$ E38

,

or, in matrix form:

$$\begin{array}{l}\left[\begin{array}{cccccc}\cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ K_{uxx}^{\left(j\right)}& K_{vyx}^{\left(j\right)}& \left(K_{uyx}^{\left(j\right)}+K_{vxx}^{\left(j\right)}\right)& 0& 0& 0\\ 0& 0& K_{uxx}^{\left(j\right)}& 0& K_{vyx}^{\left(j\right)}& \left(K_{uyx}^{\left(j\right)}+K_{vxx}^{\left(j\right)}\right)\\ 0& 0& K_{uxy}^{\left(j\right)}& 0& K_{vyy}^{\left(j\right)}& \left(K_{uyy}^{\left(j\right)}+K_{vxy}^{\left(j\right)}\right)\\ 0& K_{uxy}^{\left(j\right)}& 0& K_{vyy}^{\left(j\right)}& \left(K_{uyy}^{\left(j\right)}+K_{vxy}^{\left(j\right)}\right)& 0\\ K_{uxy}^{\left(j\right)}& K_{vyy}^{\left(j\right)}& \left(K_{uyy}^{\left(j\right)}+K_{uxx}^{\left(j\right)}+K_{vxy}^{\left(j\right)}\right)& 0& K_{vyx}^{\left(j\right)}& \left(K_{uyx}^{\left(j\right)}+K_{vxx}^{\left(j\right)}\right)\\ 0& K_{uxx}^{\left(j\right)}& K_{uxy}^{\left(j\right)}& K_{vyx}^{\left(j\right)}& \left(K_{uyx}^{\left(j\right)}+K_{vyy}^{\left(j\right)}+K_{vxx}^{\left(j\right)}\right)& \left(K_{uyy}^{\left(j\right)}+K_{vxy}^{\left(j\right)}\right)\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \end{array}\right]\left\{\begin{array}{c}{A}_{11}\\ A_{12}\\ A_{16}\\ A_{22}\\ A_{26}\\ A_{66}\end{array}\right\}=\\ =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\rho h{\varpi }_j^2\left\{\begin{array}{c}\cdots \\ {\displaystyle {\int }_{\Omega }{U}^{\left(j\right)}{x}^{2}d\Omega }\\ {\displaystyle {\int }_{\Omega }{V}^{\left(j\right)}{x}^{2}d\Omega }\\ {\displaystyle {\int }_{\Omega }{U}^{\left(j\right)}{y}^{2}d\Omega }\\ {\displaystyle {\int }_{\Omega }{V}^{\left(j\right)}{y}^{2}d\Omega }\\ {\displaystyle {\int }_{\Omega }{U}^{\left(j\right)}\left(2xy\right)d\Omega }\\ {\displaystyle {\int }_{\Omega }{V}^{\left(j\right)}\left(2xy\right)d\Omega }\\ \cdots \end{array}\right\}\end{array}$$ E39

Eq. (39) can also be rewritten in a compact form as:

$$KA=C$$ E40

,

where, considering L modes, K is a 6L x 6 matrix in Eq. (40), A is a 6 x 1 matrix, and C is a 6L x 1 matrix in Eq. (40). Eq. (40) is an over determined system. Thus, the solution can be found by least squares:

$$A={\left(K^{\text{T}}K\right)}^{-\text{1}}\left(K^{\text{T}}C\right)$$ E41

,

from where the extensional elastic constants A _ij are computed.

3. Results and comments

A commercial finite element code (ANSYS 11.0) was used to give particular mode shapes and their corresponding natural frequencies from both in-plane and out-of-plane numerical modal analysis. Element SHELL99 was used and plates under free-edge boundary conditions were considered.

To exemplify the method proposed by Grédiac and Paris (1996), it was used an anisotropic plate with dimensions 0.450 x 0.350 x 0.0021 m and density 1500 kg/m³. It was used a laminate with 8 plies, [0 45 90 135]_S, and the following engineering elastic constants by ply: E ₁=120 GPa (Young’s module along the principal direction 1), E ₂=10 GPa (Young’s module along the principal direction 2), G ₁₂=4.9 GPa (shear module along the plane 1-2), and ν ₁₂=0.3 (Poisson’s ratio along the plane 1-2). A mesh of 651 nodes was used. Fig. 2 shows the three modes used to identify the required properties. As can be seen, depending on type of material anisotropy, it is not possible to find all three modes necessary to apply the method. In this case, it is necessary to find the more approximated ones.

Figure 2.

Numerical modes obtained from Ansys: (a) Mode shape 1, (b) Mode shape 2, (c) Mode shape 3.

Table 1 shows the bending elastic stiffness constants computed using the engineering constants and the classical theory of laminates, and it also shows the errors computed after applying the identification method. As can be observed, the technique is able to find very satisfactory results when it is used the correct modes. The problem of this technique is the high sensitivity to noise presence because of second-order derivatives. More results and comments about this method can be found in Grédiac& Paris, 1996.

In order to verify the accuracy of the extensional elastic stiffness identification method, it was used six graphite/polymer symmetric laminated plates (Table 2): a fully anisotropic with all A _ij different from zero (i,j = 1, 2, and 6); a cross-ply orthotropic with A ₁₁ = A ₂₂, A ₁₆ = A ₂₆ = 0, and A ₁₁ – A ₁₂ ≠ 2A ₆₆; a 0° unidirectional orthotropic; a 30° unidirectional orthotropic (generallyorthotropic); a +30°/-30 angle-ply orthotropic; and a quasi-isotropic with A ₁₁ = A ₂₂, A ₁₆ = A ₂₆ = 0, and A ₁₁ – A ₁₂ = 2A ₆₆. These laminates have 8 plies with the following engineering elastic constants by ply: E ₁=155 GPa, E ₂ =12.10 GPa, G ₁₂=4.4 GPa, and ν ₁₂=0.248. Dimensions considered were 0.450 x 0.350 x 0.003 m (x, y, and z plate coordinate axis, respectively), for the rectangular plate, 0.350 x 0.350 x 0.003 m, for the square plate, and, the density material was 1500 kg/m³. The plates were modeled using a mesh with 651 nodes, for the rectangular plates, and with 441 nodes, for the square plate. The extensional elastic stiffness constants are shown in Table 2. The terms “Aniso”, “Ortho” and “Quasi-iso” are simplifications of “Anisotropic”, “Orthotropic” and “Quasi-isotropic”, respectively. Only the first fifteen mode shapes were analyzed.

As can be seen in Table 2, the constants A ₁₂ and A ₆₆ for the orthotropic laminates I and IIare equals. They also are equals to laminates with the same characteristics but with 90° unidirectional fibers. 0° and 90° unidirectional laminated plates are only different in relation to A ₁₁ and A ₂₂ terms. These are inverted: A ₁₁ term of the 0° laminate is equal to A ₂₂ term of the 90° laminate, and vice-verse. For the generally orthotropic laminate III and orthotropic laminate IV the difference are only the A ₁₆ and A ₂₆ constants: they are nulls for the laminate IV and non-nulls for the laminate III. For laminate III all extensional elastic constants are non-nulls, similar to fully anisotropic laminates, what, consequently, originates to full extensional elastic stiffness matrix.

Bending elastic constants	N x mm	Errors (%)
D 11	64363.9	0.02
D 22	24155.8	0.04
D 12	8875.1	0.02
D 66	10032.7	1.22
D 16	6019.6	0.63
D 26	6019.6	0.64

Table 1.

Bending elastic stiffness constants of the tested anisotropic plate and the computed errors.

The key point of this technique of identification is related with the correct choice of mode shapes together to the weighting functions (virtual fields). The correct mode shapes are called here by “suitable modes” and the correct combination between these modes and the weighting functions are called by “suitable combinations”. The identification of the suitable modes is not difficult, as it will be shown in the next topics. But, the suitable combinations are more difficult because they depend on the type and the geometry of the material. Fortunately, there are some aspects that help finding the best choice. Unlike the bending stiffness identification method originally proposed, for this method there are a lot of modes and suitable combinations that give satisfactory results.

Aij; 10^8 [N/m]	Aniso; [90 0 0 45]S	Ortho I; [90 0 90 0]S	Ortho II; [0]4S	Ortho III; [30]4S	Ortho IV; [30 -30]2S	Quasi-iso; [90 45 0 -45]S
A 11	2.7865	2.5186	4.6724	2.7840	2.7840	1.9776
A 12	0.3610	0.0905	0.0905	0.9020	0.9020	0.6315
A 16	0.2692	0	0	1.4012	0	0
A 22	1.7096	2.5186	0.3648	0.6301	0.6301	1.9776
A 26	0.2692	0	0	0.4641	0	0
A 66	0.4025	0.1320	0.1320	0.9436	0.9436	0.6730

Table 2.

Laminated plates used for the verification of the method

[90 0 0 45]S anisotropic plate

Table 3 shows some errors computed for the anisotropic plate, rectangular and square. It was considered only the first fifteen in-plane modes shapes. Anisotropic plates, in general, give very satisfactory results using the combinations among suitable modes. This factor can be justified by the fact of these combinations be hardly involved with all required extensional elastic constants A _ij’s.

The numerical contribution of each mode to the computation of a specific constant cannot be jeopardized by numerical contribution of another mode during the solution of the system given by Eq. (40).The suitable modes are those that when associated with weighting functions do not null or give very low values for integrals of the right (K matrix) and/or left (C matrix) hand sides of Eq. (40). The suitable combinations are one composed by suitable modes and that give more accurate results. In the majority of the cases, combinations using a higher number of suitable modes can be suitable combinations. According to Table 3 is possible to see that using combinations with only two suitable modes very satisfactory results can be obtained. Satisfactory results would also be obtained using combinations with any modes since the number of suitable modes among all used modes is higher than non-suitable modes. But the accuracy of these results cannot be guaranteed for all combinations.

[0 90 0 90] _S orthotropic plate (ortho I)

In general, for the orthotropic and isotropic materials is more difficult to find the suitable combinations when it is compared to fully anisotropic materials. It is necessary to take care to correctly identifying the combinations that give the best results. In these types of materials not all combinations are among suitable modes that can be considered as being suitable combinations. According to values found to terms of the K and C matrices, Eq. (40), and using combinations among suitable modes, it is possible to see the following types of systems:

For the rectangular plate:

Type 1:

$$\left[\begin{array}{cccccc}{I}_1& 0& 0& I_4& 0& 0\\ 0& 0& I_3& 0& 0& 0\\ 0& I_2& 0& 0& I_4& 0\\ 0& I_1& 0& 0& I_4& 0\\ 0& 0& 0& 0& I_3& 0\\ 0& 0& I_3& I_1& 0& I_4\end{array}\right]\left\{\begin{array}{c}{A}_{11}\\ A_{16}\\ A_{66}\\ A_{12}\\ A_{26}\\ A_{22}\end{array}\right\}=\left\{\begin{array}{c}{I}_5\\ I_6\\ 0\\ 0\\ 0\\ I_7\end{array}\right\}$$ E42

Type 2:

$$\left[\begin{array}{cccccc}0& I_2& 0& 0& 0& 0\\ 0& I_1& 0& 0& I_3& 0\\ I_1& 0& I_2& I_3& 0& 0\\ 0& 0& I_2& 0& 0& 0\\ 0& 0& 0& I_1& 0& I_3\\ 0& I_1& 0& 0& I_4& 0\end{array}\right]\left\{\begin{array}{c}{A}_{11}\\ A_{16}\\ A_{66}\\ A_{12}\\ A_{26}\\ A_{22}\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ I_5\\ I_6\\ I_7\\ 0\end{array}\right\}$$ E43

Anisotropic rectangular plate – suitable modes: 2, 3, 6, 9, 10, 13, and 14
Suitable; combinations	Errors (%)
A 11	A 12	A 22	A 16	A 26	A 66
2-3-6-9	2.20	1.95	1.56	1.89	1.46	0.78
2-3-6	0.04	0.20	0.21	0.93	0.79	0.97
2-3-9	0.37	2.63	0.27	0.78	0.93	0.28
2-6-9	0.30	2.46	0.17	1.45	1.03	0.02
3-6-9	0.17	1.77	0.11	0.97	1.09	0.21
2-3	0.00	0.32	0.24	0.78	0.53	1.08
6-9	0.01	2.06	0.11	1.62	1.50	1.51
Anisotropic square plate – suitable modes: 2, 3, 6, 9, 10, 13, and 14
Suitable; combinations	Errors (%)
A 11	A 12	A 22	A 16	A 26	A 66
2-3-6-9	3.73	3.41	2.23	3.21	1.88	1.06
2-3-6	0.20	0.62	0.23	0.86	0.09	1.08
2-3-9	0.41	3.86	0.54	1.36	0.77	0.50
2-6-9	0.35	3.43	0.44	2.00	1.61	0.13
3-6-9	0.24	2.59	0.14	1.66	0.84	0.30
2-3	0.04	0.06	0.17	1.02	0.97	1.49
6-9	0.08	2.66	0.25	2.70	1.08	1.93

Table 3.

Errors computed for the anisotropic plate using some suitable combinations.

whereI ₁, I ₂, I ₃, I ₄, I ₅, I ₆, and I ₇ are the integral values of Eq. (40). These integrals can be negatives or positives depending on strain direction and reference coordinate axis. As can be observed, for these two types of systems, Eq. (42) and Eq. (43), all elastic constants are involved, and, however, combinations associated to only one type of system can be suitable combinations and sufficient to give correct results. For this plate were found the following suitable modes: 2, 3, 6, 8, 11, and 14.

For the square plate: for this plate, the suitable modes are: 2, 3, 6, 7, 11, and 12. The systems of equations are full, even though of the additional terms, that are null in rectangular plate, to be low in this square plate. It is observed in this plate that for each suitable mode there is another identical but out-of-phase at 90°: modes 2 and 3, 6 and 7, and, 11 and 12.

Table 4 shows errors computed to some suitable combinations for these orthotropic plates. As can be seen, very satisfactory results can be obtained using correct combinations of modes. For this type of orthotropy, it can be more difficult to compute an accurate value for constant A ₁₂. The identification of the suitable combinations is not so clear. As this technique of identification is associated to the solution of an equation system having different modal contributions, it is difficult to identify which modes compose a correct combination.

Ortho I rectangular plate – suitable modes: 2, 3, 6, 8, 11, and 14
Suitable; combinations	Errors (%)				Differences (N/m)
A 11	A 12	A 22	A 66	A 16	A 26
2-3-6-8-11-14	0.14	3.80	0.12	1.22	-2.7 x 10^4	0.8 x 10^4
2-3-6-8	0.11	2.90	0.24	0.74	-2.5 x 10^4	0.6 x 10^4
3-6-8-11	0.17	1.58	0.31	1.67	-1.4 x 10^4	0.5 x 10^4
2-8-11	0.17	0.29	0.23	1.46	-1.5 x 10^4	0.6 x 10^4
6-8-11	0.18	0.89	0.34	1.91	-4.8 x 10^4	1.0 x 10^4
3-8	0.02	5.75	0.91	0.13	-1.4 x 10^4	8.4 x 10^3
8-11	0.25	0.42	0.40	2.17	-2.6 x 10^4	1.1 x 10^4
Ortho I square plate – suitable modes: 2, 3, 6, 7, 11, and 12
Suitable; combinations	Errors (%)				Differences (N/m)
A 11	A 12	A 22	A 66	A 16	A 26
2-3-6-7-11-12	0.09	5.81	0.05	1.89	2.7 x 10^4	-1.6 x 10^4
2-3-6-11	0.10	5.09	0.06	1.92	-4.4 x 10^4	-8.7 x 10^4
3-6-7-12	0.38	2.33	0.08	2.92	6.4 x 10^4	-1.7 x 10^4
2-6-11	0.13	5.47	0.11	2.51	-1.0 x 10^5	-1.8 x 10^4
2-7-12	0.22	5.37	0.04	1.95	1.1 x 10^5	0.5 x 10^5
2-11	0.02	9.04	0.18	0.56	-1.3 x 10^4	-1.9 x 10^4
6-11	0.18	5.63	0.30	3.46	-1.1 x 10^5	-2.7 x 10^4

Table 4.

Errors computed for the ortho I plate using some suitable combinations.

[0] _4S orthotropic plate (ortho II)

For the rectangular plate: for this plate, modes 3, 10, 11, and 14 originate systems of type 1, Eq. (42), and modes 2, 5, 9, and 15 originate systems of type 2, Eq. (43).

For the square plate: for this plate, modes 4, 10, 12, and 14 originate systems of type 1, Eq. (51), and modes 2, 5, 9, and 13 originate systems of the type 2, Eq. (43).

Table 5 shows the errors computed to some suitable combinations for these orthotropic plates. Using correct combinations very satisfactory results can be obtained. Similar to ortho I plate, for this type of orthotropy, constant A ₁₂is more difficult of being accurately computed.

Ortho II rectangular plate – suitable modes: 2, 3, 5, 9, 10, 11, 14, and 15
Suitable; combinations	Errors (%)				Differences (N/m)
A 11	A 12	A 22	A 66	A 16	A 26
2-3-5	2.77	6.06	0.08	0.38	2.4 x 10^4	0.5 x 10^4
11-14-15	1.14	2.15	0.92	1.15	2.3 x 10^5	0.1 x 10^5
2-3	2.36	0.57	0.19	0.72	1.7 x 10^4	0.1 x 10^4
2-10	2.10	0.82	0.60	2.42	-2.0 x 10^5	-0.4 x 10^5
2-14	2.04	5.22	0.37	1.12	-3.9 x 10^4	-1.2 x 10^4
11-14	1.64	2.14	0.87	0.98	-1.8 x 10^5	-0.1 x 10^4
Ortho II square plate – suitable modes: 2, 4, 5, 9, 10, 12, 13, and 14
Suitable; combinations	Errors (%)				Differences (N/m)
A 11	A 12	A 22	A 66	A 16	A 26
2-4-5	4.13	8.63	0.39	0.45	-3.5 x 10^4	-0.1 x 10^4
2-4-14	2.43	0.78	0.72	2.89	4.2 x 10^4	1.1 x 10^4
2-4	3.35	0.10	0.21	1.25	-3.2 x 10^4	-0.1 x 10^4
2-14	1.63	1.59	1.30	4.72	1.1 x 10^5	0.3 x 10^5
5-14	2.67	13.54	2.21	2.82	5.6 x 10^4	-1.0 x 10^4
9-14	2.12	62.05	0.85	3.89	1.1 x 10^5	0.2 x 10^5

Table 5.

Errors computed for the ortho II plate using some suitable combinations.

[30] _4S orthotropic plate (ortho III)

For these plates, rectangular and square the computed matrices K and C, Eq.(40), are full matrices, similar ones of the anisotropic plates. Thus, the majority of combinations among suitable modes are suitable combinations. Combinations that are not suitable present high errors for all constants, what, consequently, make them easy to be identified. Table 6 shows the errors computed to some suitable combinations. Using correct combinations very satisfactory results can be obtained.

Ortho III rectangular plate – suitable modes: 2, 3, 6, 8, 10, 12, and 13
Suitable; combinations	Errors (%)
A 11	A 12	A 22	A 66	A 16	A 26
3-6-12	0.25	0.51	0.10	1.05	0.69	0.99
8-10-12	4.26	5.42	3.67	3.43	4.60	5.62
2-3	0.85	1.72	1.22	1.39	0.10	0.96
3-6	0.63	0.32	0.04	0.22	0.53	0.14
6-12	0.71	1.06	0.05	1.39	1.19	1.30
10-13	4.42	2.05	0.71	4.07	5.21	0.63
Ortho III square plate – suitable modes: 2, 4, 6, 8, 9, 12, and 13
Suitable; combinations	Errors (%)				Differences (N/m)
A 11	A 12	A 22	A 66	A 16	A 26
2-4-6	2.62	3.94	3.18	0.59	2.02	3.23
2-4-8	1.20	2.40	1.64	1.27	0.11	1.19
2-4	2.35	1.11	0.02	5.46	4.26	2.90
2-6	0.76	2.76	3.52	0.01	0.53	2.90
2-13	1.75	2.75	1.51	1.92	0.33	1.25
12-13	3.08	2.41	2.63	0.53	1.85	2.54

Table 6.

Errors computed for the ortho III plate using some suitable combinations.

[30 -30 30 -30] _S orthotropic plate (ortho IV)

For the rectangular plate: it was found the following suitable modes: 2, 4, 5, 8, 11, 12 and 14. For this plate, the majority of combinations among the suitable modes are suitable combinations.

For the square plate: it was found the following suitable modes:1, 4, 6, 9, 11, 12, and 14. Similar to rectangular plate, here the most of combinations among the suitable modes are suitable combinations.

Table 7 shows errors computed to some suitable combinations. Using correct combinations, very satisfactory results can be obtained.

Ortho IV rectangular plate – suitable modes: 2, 4, 5, 8, 11, 12, and 14
Suitable; combinations	Errors (%)				Differences (N/m)
A 11	A 12	A 22	A 66	A 16	A 26
2-4-5	1.52	2.15	1.28	1.54	-1.7 x 10^4	-0.6 x 10^4
2-4-8	1.50	2.09	0.83	1.14	1.4 x 10^4	0.9 x 10^4
5-11-12	0.27	0.09	0.04	1.11	-1.2 x 10^4	1.9 x 10^4
2-4	1.18	1.83	1.02	2.87	-1.2 x 10^4	-0.4 x 10^4
2-5	0.49	1.01	1.01	1.41	-1.8 x 10^4	-0.7 x 10^4
5-11	0.90	0.37	0.04	1.05	3.1 x 10^4	0.2 x 10^4
Ortho IV square plate – suitable modes: 1, 4, 6, 9, 11, 12, and 14
Suitable; combinations	Errors (%)				Differences (N/m)
A 11	A 12	A 22	A 66	A 16	A 26
1-6-12	0.16	0.30	0.13	3.07	1.7 x 10^4	2.5 x 10^4
6-9-12	0.18	0.35	0.30	2,68	1.3 x 10^4	4.8 x 10^4
9-12-14	0.18	1.87	0.40	0.68	-6.2 x 10^4	8.1 x 10^4
1-12	0.59	0.26	0.74	2.79	2.0 x 10^4	2.9 x 10^4
4-11	0.94	1.55	2.77	0.23	2.2 x 10^4	0.8 x 10^4
12-14	0.26	2.17	0.19	0.90	-9.7 x 10^4	9.4 x 10^4

Table 7.

Errors computed for the ortho IV plate using some suitable combinations.

[90 45 0 -45] _S quasi-isotropic plate

For the rectangular plate: it was found the following suitable modes: 1, 4, 6, 8, 11, 12, and 14. For this plate the most of combinations among the suitable modes are suitable combinations.

For the square plate: it was found the following suitable modes: 2, 3, 7, 8, 10, and 11. Similar to rectangular plate, here the majority of combinations among the suitable modes are suitable combinations.

Table 8 shows errors computed to some suitable combinations. Using correct combinations, very satisfactory results can be obtained.

Quasi-iso rectangular plate – suitable modes: 1, 4, 6, 8, 11, 12, and 14
Suitable; combinations	Errors (%)				Differences (N/m)
A 11	A 12	A 22	A 66	A 16	A 26
1-4-6	0.60	1.47	0.19	0.19	9.1 x 10^3	4.7 x 10^3
1-4-8	0.09	0.80	0.11	0.96	7.4 x 10^3	9.5 x 10^3
1-4-11	0.33	2.04	0.26	0.08	1.4 x 10^3	8.1 x 10^3
1-4	0.25	0.54	0.29	1.93	5.7 x 10^4	7.1 x 10^4
1-6	0.28	0.91	1.77	0.00	1.0 x 10^4	0.1 x 10^4
4-8	0.73	0.39	0.05	0.95	1.1 x 10^4	0.3 x 10^4
Quasi-iso square plate – suitable modes: 2, 3, 7, 8, 10, and 11
Suitable; combinations	Errors (%)				Differences (N/m)
A 11	A 12	A 22	A 66	A 16	A 26
2-3-7	0.37	1.33	0.35	0.91	4.7 x 10^5	4.7 x 10^5
2-3-8	0.32	1.36	0.33	0.97	-4.4 x 10^5	-4.4 x 10^5
2-10-11	0.07	3.29	0.45	0.81	-2.3 x 10^5	-1.1 x 10^5
2-3	0.33	0.79	0.34	2.68	2.2 x 10^4	1.3 x 10^4
3-7	0.56	1.73	0.43	0.80	6.1 x 10^5	5.2 x 10^5
8-11	0.89	3.05	0.06	1.34	1.1 x 10^5	-1.2 x 10^5

Table 8.

Errors computed for the quasi-isotropic plate using some suitable combinations.

Figs. 3 to 14 show the fifteen first in-plane mode shapes to all the analyzed plates: rectangular and square geometry.

Figure 3.

Fifteen first in-plane mode shapes and natural frequencies to anisotropic rectangular plate.

Figure 4.

Fifteen first in-plane mode shapes and natural frequencies to anisotropic square plate.

Figure 5.

Fifteen first in-plane mode shapes and natural frequencies to ortho I rectangular plate.

Figure 6.

Fifteen first in-plane mode shapes and natural frequencies to ortho I square plate.

Figure 7.

Fifteen first in-plane mode shapes and natural frequencies to ortho II rectangular plate.

Figure 8.

Fifteen first in-plane mode shapes and natural frequencies to ortho II square plate.

Figure 9.

Fifteen first in-plane mode shapes and natural frequencies to ortho III rectangular plate.

Figure 10.

Fifteen first in-plane mode shapes and natural frequencies to ortho III square plate.

Figure 11.

Fifteen first in-plane mode shapes and natural frequencies to ortho IV rectangular plate.

Figure 12.

Fifteen first in-plane mode shapes and natural frequencies to ortho IV square plate.

Figure 13.

Fifteen first in-plane mode shapes and natural frequencies to quasi-iso rectangular plate.

Figure 14.

Fifteen first in-plane mode shapes and natural frequencies to quasi-iso square plate.

4. Conclusions

The identification of elastic properties using VFM has shown to be a very efficient technique since the correct combinations among mode shapes and weighing functions are used. This factor is the key point to find the correct results. However, the identification of these suitable combinations is not so simple in some situations, mainly to the extensional elastic stiffness identification method. Fortunately, there are some characteristics that can help to find such combinations, as it was shown here. A great advantage of this method is related to the large number of possibilities to make combinations able to give very satisfactory results.