Abstract
Using experimental data and numerical simulations, a new combined technique is presented for characterization of thin and thick orthotropic composite laminates. Four or five elastic constants, as well as ply orientation angles, are considered as the unknown parameters. The material characterization is first examined for isotropic plates under different boundary conditions to evaluate the method’s accuracy. The proposed algorithm, so-called CPAM (Combined Programs of ABAQUS and MATLAB), utilizes an optimization procedure and makes simultaneous use of vibration test data together with their corresponding numerical solutions. The numerical solutions are based on a commercial finite element package for efficiently identifying the material properties. An inverse method based on particle swarm optimization algorithm is further provided using MATLAB software. The error function to be minimized is the sum of squared differences between experimental and simulated data of eigenfrequencies. To evaluate the robustness of the model’s results in the presence of uncertainty and unwanted noises, a sensitivity analysis that employs Gaussian disorder model is directly applied to the measured frequencies. The results with high accuracy confirm the validity and capability of the present method in simultaneous determination of mechanical constants and fiber orientation angles of composite laminates as compared to prior methods.
1. Introduction
In many engineering applications, such as design or quality control of advanced composite structural systems, it is particularly important to have correct understanding of mechanical constants of composites [1]. It is more complicated to determine the mechanical constants of fiber-reinforced polymer composites than to determine those of homogenous and isotropic materials, because, in such materials, which may include a variety of constituent materials mixed together, the number of elastic constants is increased and also some additional complexity due to their nonhomogenous nature may occur. When static test methods such as simple tensile test are used for composites, they will be faced with some problems. Special problems such as effects of supports, dependence on samples size, and nonuniform fields of stress-strain in the experiments usually cause a wide scatter and unwanted errors in the test results. In addition, the destructive nature of these tests and inability to repeat the test for a specific sample are other disadvantages of conventional experimental methods for determining the material properties. Thus, a powerful and useful technique that is able to be utilized for characterizing the real structures needs to be introduced. To resolve aforementioned problems, among previously proposed methods, the nondestructive evaluations of material parameters based on inverse computational techniques look promising for composite materials [2].
In the inverse computational methods, there are complex relationships between material properties and structural behavior which allow us to characterize the material properties of composites. These relationships are usually represented by a typical computational and mathematical model known as direct simulated solution. In the direct problem, measureable data are determined from model parameters. So, if there is a certain set of measured experimental data for evolution of structural behavior, the composite plates’ properties can be determined by solving an inverse problem [3]. But, in general, it cannot be measured exactly and the measured data are often diverted from real measurable quantities. This deviation can be as small as rounding error that is created by computer estimations. In fact, the deviation can be completely innate in measuring process, which means it is due to existing errors in the measurement process. So, data from model parameters are obtained without any errors in the direct process but, in general, data are in the deviated form.
Extraction of natural frequencies through modal analysis has been considered as experimental data in the inverse methods. Literature reveals that, in order to investigate the stability of solution to small changes imported to the experimental data, a sensitivity analysis needs to be applied accompanied with the various inverse methods proposed for determination of parameters of composite plates. When laminated composites are designed, design variables for a specific kind of material are fibers directions, number of layers, and thickness of each layer. Because of restrictions on the manufacture and use of these composite plates, choosing of fibers direction is usually limited to a set of definite discrete values and also layers’ thickness is specific [4]. Moreover, among the various factors, stacking sequence of layers is one of the most important parameters that can increase structures’ resistance to the damage caused by different loads. The orientation of the fibers in layers highly affects the mechanical strength of laminated composites, so identification of the stacking sequence of layers for theoretical research and practical applications seems worthy [5]. Sepahvand and Marburg [6] attempted to determine elastic constants of orthotropic plates by expansion of Polynomial Chaos theory and also experimental data resulting from modal analysis. Parameters of their model are obtained by expansion of generalized Polynomial Chaos and the corresponding unknown constants are calculated using inverse stochastic problem. Ip et al. [7] determined elastic constants of cylindrical orthotropic shells by using Bayesian estimation based on natural frequencies obtained from a model with free boundary conditions. Their model is formulated by using Rayleigh-Ritz method that is based on the characteristic functions of the beam. Then they succeeded to realize the high dependence of circumferential wave pattern of vibrational mode and in-plane shear modulus on the natural frequencies of cylindrical shells. To simplify the modeling process and reduce the complexity of numerical modeling, Hwang and Chang [8] used an FEM model with an optimization process to determine just elastic constants of thin and thick composite laminates and also aluminum plates. Deobald and Gibson [9] determined four independent elastic constants of thin orthotropic plates. The Classical Lamination Theory and three-mode Rayleigh formulation optimized with objective function of sum of square difference of experimental and numerical frequencies are used in their work. Liu et al. [10] used an advanced neural network algorithm as optimization algorithm in their proposed inverse method to determine nonhomogenous properties of composite laminates. They used displacement response of surfaces due to linear loading as experimental data and Hybrid Numerical Method (HNM) as direct solution for evaluating the dynamic response of surfaces. Finally, by applying three disorder levels of 1, 3, and 10 percent to the measured experimental data, they could use these simulated data to evaluate the amount of solution stability. Mohan Rao and Arvind [11], in their first case study, determined stacking sequence of the layers of hybrid laminated composite panels in order to maximize critical buckling load under thermal loading using a metaheuristic optimization algorithm known as scatter search. They used a same algorithm in order to determine stacking sequence of the multilayer plates in their second case study and minimize the weight and manufacturing cost of these plates by considering constraints on the buckling load and the first natural frequency.
As it is evident from the literature, simultaneous identification of mechanical constants and stacking sequence of composite plates has not been reported and the presented study tries to perform this task. In the present study, natural frequencies extracted from vibration analysis of structure are considered as experimental results. Accuracy of the obtained results by inverse method depends commonly on some issues like sufficient accuracy of the method used in vibration analysis and solution of eigenproblem. In this study, the commercial FEM code of ABAQUS is used for vibration modeling and calculation of natural frequencies of thin and thick homogenous and orthotropic plates. In order to determine four elastic constants and stacking sequence of such composite plates by inverse method, the PSO algorithm in MATLAB software is employed. The sum of squared difference of experimental and numerical frequencies is used as the objective function. By linking MATLAB software and ABAQUS software, it is possible to perform this alternating optimization process. In order to apply sensitivity analysis, the Gaussian noise is added directly to the (simulated) measured frequencies and then robustness of the proposed inverse method to the unwanted disorders is evaluated. Outcomes show that the proposed inverse method is properly stable towards adding noise. Results obtained from the present method are validated by comparing those obtained from other previously proposed methods.
2. Vibration Analysis of the Composite Plates (Direct Solution)
A modeling method which assumes that the structural properties of the material are known and consequently predicts the physical behavior of the structure through them is usually referred to as direct modeling or forward solution [3]. As an example, for a transversely isotropic composite plate, direct modeling of a vibration problem can consist of a model which is able to provide natural frequencies and mode shapes of a two-dimensional thin composite plate with four independent elastic constants and also arbitrary fiber orientation angle for each layer. For this purpose, an available commercial finite element code can be used for undamped free vibration modeling of the plate.
In the present study, natural frequencies and modes shapes of the considered composite plates have been extracted using modal analysis in ABAQUS software (version 6.12-1). All models are meshed by eight-node plate element (S8R element) with six degrees of freedom for each node. This element can be used for both thin and thick composite plates and also can be utilized for homogenous and orthotropic materials [13].
3. Optimization Design Based on PSO Algorithm
The problem of determining the parameters of composite plates by inverse method can be discussed as an optimization problem. One of the most powerful algorithms that get noticed by researchers is PSO algorithm [14, 15]. In this algorithm, each design variable of the population that is known as particle possesses a velocity that enables it to move through the problem space and is replaced instead of crossover and mutation of other population-based algorithms (which search more than one initial solution of the search space for design variables) like genetic algorithms [14]. Velocity can be defined by a vector in the problem space. So, each particle is characterized by a location (solution) and velocity. Relocation of a particle is done by means of the present location of the particle and its velocity. It means that resultant of the current location (solution) and velocity vector denotes the location of each particle in the next iteration. Each particle maintains its current location, velocity, and a memory of its previous best location, as well as the global best location of the entire swarm that all particles have achieved so far at each iteration. The velocity vector and displacement of each particle in the next iteration can be estimated, respectively, through the following relations:where is the inertia weight, and are positive constants whose names are social and self-cognitive parameters, respectively, and in which m is the population size. and are uniformly distributed random values within , and indicates the repetition number and is the maximum allowable iteration. The first right-hand term in (1) is previous velocity of the particle which enables it to move through the problem space. The second and third terms show the velocity difference based on the best solution of each particle and the best solution of all particles, respectively [15]. Obviously, the greater social parameter value of the particle is, the more velocity vector approaches the best answer of all particles and also the greater self-cognitive parameter is, the more velocity vector approaches the best answer that the particle has experienced. In general, PSO algorithm can be summarized in three phases: evaluation, comparison, and imitation. Evaluation phase measures desirability of each particle error. Comparison phase determines the best particle among all particles and, in the imitation phase, the new location of the particles based on knowledge that has been obtained until now is determined. These three phases are allowed to continue until the end-up condition is satisfied. The ultimate goal is to find the particle that has the best answer for investigation problem.
FPSO (Fuzzy-Particle Swarm Optimization) is another type of PSO algorithm for identification of discrete parameters. In this case, the set of location vectors can be expressed in a fuzzy matrix where its elements represent fuzzy logic membership functions. How to use this algorithm for job scheduling problem is briefly explained here. is assumed to be a set of jobs that must be processed on a set of machines . All jobs of set must be carried out and so the total cost of jobs is minimized.
In FPSO algorithm, the relations between jobs and machines can be defined in a matrix form aswhere is fuzzy logic and represents dependency of on . So, amount of elements of matrix must be in range . In another view, shows the possibility that carries out job . In other words, sum of each column’s elements of the matrix must be equal to 1. In a similar way, the velocity matrix can be defined as
Equations (1) and (2) can be rewritten for FPSO algorithm as
where and are of order matrix and their elements are uniformly distributed random variables in range . Operator multiplies the corresponding elements of matrixes together. There is not any changes in other variables. According to the constraints on the elements of matrix , the illegal values of that may be generated after updating the fuzzy matrix must be modified. To this end, first, all elements of matrix must be positive. So, the negative elements generated in matrix are intentionally set to zero. Second, sum of each column’s elements must be equal to 1. So each element of this matrix must be divided to the sum of its column elements.
4. Definition of the Optimization Problem
Design variables considered in this study are four independent elastic constants of thin transversely isotropic plates and also fiber orientation angles. These constants are which represent Young’s modulus parallel and perpendicular to the fiber orientation, in-plane Poisson’s ratio, and the shear modulus, respectively. Constants are required for FE modeling in ABAQUS, which can be defined by the following equations for transversely isotropic material system [15]:
For thick plates, is conceded as fifth independent constant and (6) should not be used to evaluate this constant. The objective function that must be minimized by the optimization algorithm is defined as where represent the measured experimental frequencies and are the corresponding frequencies calculated from the numerical solution. The search space that stands between upper and lower bounds of elastic constants is represented as
It is to be noted that the stiffness matrix (constitutive relationships between stress and strain components) must be positive definite. Therefore, there are some relations between elastic constants of the composite material. These relations are defined as the following constraints for the algorithm [12]:
As it is mentioned, fibers’ orientation is limited to some definite values; so they are considered as discrete parameters in this study. The fiber orientation angles between which the algorithm is allowed to search are given in the following equation:
5. Settings for Optimization Algorithm
5.1. Initial Population
For generating the initial starting particles, the first value of is randomly generated within the specified constraints. Since range of is wider than and according to the fourth constraint of (9), should be larger than ; its value is randomly generated between and the upper limit of . In order to determine the reasonable amount of , according to the third constraint, it is also randomly determined within interval . Because there is not any constraint for orientation of fibers’ angles, the initial population of them is created randomly.
5.2. Evaluation of Solution and Error Analysis
According to initial solution of particles which has been generated within the feasible region of the third and fourth constraints, it is not necessary to check satisfaction of these two constraints in the initial evaluation. But because there is no guarantee that they remain valid in the next iterations, feasibility of them must be also checked alongside the two other constraints in the evaluation process. Then if any constraint is violated, error of solutions will accept a great value and thus the corresponding particle will be penalized. But if the solution is reasonable, it will be sent to ABAQUS solver for proceeding the calculations. Sum of squared deviation of outputs from the target values is shown as error of the desired solution.
6. Employing the Inverse Method for Determination of Parameters
For determination of parameters, it is essential that the direct modeling explained in the previous sections be linked with code written in MATLAB for optimization algorithm. In other words, there must be an ability to make a linking bridge between ABAQUS and MATLAB software in order to call both frequencies provided by ABAQUS in MATLAB and also parameters provided by MATLAB in ABAQUS. This alternating process is therefore repeated between ABAQUS and MATLAB software until the error function defined in (7) is minimized. In the following, how ABAQUS and MATLAB software work is explained. To have better perception of the optimization procedure, the inverse solution algorithm is explained here step by step, and its corresponding flowchart is represented in Figure 1; it does not necessarily mean the time priority of these steps’ operation during performing the algorithm.
In Step , parameters are generated by code written in MATLAB software. In Step , a text file is created by this code and these generated parameters are saved in this file according to commands written in this code. In Step , ABAQUS software is run (without GUI) by python language that is explained in the next step (Step ). Commands required for operating of this step are also utilized in code written by MATLAB software. In Step , a python code that is intelligible for ABAQUS is written beforehand for ABAQUS software. The required commands for calling parameters in the text file of Step , importing them into ABAQUS, getting the natural frequencies generated by the model in ABAQUS, creating another text file, and finally saving these frequencies in this file, are written in this code. In Step , the model that is previously created in ABAQUS to determine the natural frequencies of plates is run. In Step , the code that is written by python language calls natural frequencies that are obtained by the model created in ABAQUS in the previous step. In Step , the code that is written by python programming language saves frequencies called by ABAQUS (previous stage) in a second interface text file. In Step , the code that is written in MATLAB calls frequencies saved in the second interface text file and then the error function in (7) is checked. Whenever this error function is minimized in every iteration of this alternating process, the process will be stopped and the corresponding parameters that generate these frequencies and make the error function be minimized are chosen as the set of optimal solutions; otherwise, this process will continue unless the error function is minimized. So, this algorithm that is capable of linking PSO algorithm in MATLAB to the model in ABAQUS is briefly called CPAM.
7. Result and Discussion
7.1. Validation of FEM Modeling
Due to use of natural frequencies calculated from numerical solution in objective function and the importance of their accuracy in determination of the parameters, first it is necessary to be sure of accuracy of direct solution for providing desired natural frequencies that are obtained by modeling in ABAQUS software. In order to validate the FEM modeling, an aluminum plate has been considered and its first six natural frequencies have been compared with those obtained by Rayleigh-Ritz method presented by Deobald and Gibson [9]. Hwang and Chang [8] have extracted them again to validate their FEM modeling in ANSYS software. Dimensions of the considered square aluminum plate are 25.4 cm × 25.4 cm × 0.316 cm, density is 2.77 gr/cm3, Young’s modulus is 72.4 GPa, shear modulus is 28 GPa, and Poisson’s ratio is 0.33. In order to match mesh size of this plate with the mesh size reported in [8, 9], the 10 × 10 in-plane mesh sizes that totally give 100 elements are considered for the plate. Natural frequencies obtained from modal analysis and impact technique [9] are shown in Table 1. In addition, natural frequencies obtained by Rayleigh-Ritz method [9] are also displayed in this table in comparison with those extracted by ANSYS [8] and ABAQUS solvers, respectively. The actual elastic constants are used in all of these three methods.
As can be seen, the natural frequencies obtained from these three methods are very close together in average. However, there is still natural frequencies that are not matched well together and it can lead to creation of inherent errors for determination of parameters of plates inversely. Nevertheless, as it is obvious from Table 1, natural frequencies obtained by modeling in ABAQUS software are closer to the experimental frequencies resulting from modal analysis than those obtained by modeling by ANSYS solver and Rayleigh-Ritz method; this means more accuracy for the present method in modeling and direct solution and finally leads to providing a higher accuracy in identifying the parameters.
For more complete validation of the present method, several examples were examined here, referring to the research work of Hwang and Chang [8] to determine the natural frequencies of homogenous and orthotropic plates with various boundary conditions and stacking sequence. Dimensions, density, and stacking sequence of these plates are available in Table 2 and their mechanical properties are in Table 3. Presented values for the mechanical properties of aluminum plates are obtained from metals handbook [12] and those for composite panels from static tests available in the literature [8]. Natural frequencies obtained from vibration test are shown in Table 4 and those obtained from modeling in ABAQUS software are shown in Table 5. The composite plates chosen from Hwang and Chang’s paper have free boundary condition at all edges. In aluminum plate A1, boundary conditions of all edges are also free but, in aluminum plate A2, boundary conditions of bigger edge are clamped and other ones are free.
Moreover, in order to include experiment in the presented study, 4 thick and thin glass/epoxy plates are investigated in modal laboratory (Figure 2). To reduce the errors of production process, these plates are placed in a vacuum condition. Dimensions, density, and stacking sequence of these plates are shown in Table 2. Plates C3 and C4 are laid on the foam devices (Figure 3) and because natural frequencies are independent of position and rigidities of these devices, boundary conditions of these plates are considered completely free [9, 16]. For thick plates D3 and D4, one edge is clamped and other ones are free (Figure 4). Because, for this kind of boundary condition, frequencies of all modes shift down, frequencies of higher mode can be investigated more accurately. As can be seen from Figures 3 and 4, laser test is used in this experiment to reduce experimental errors. Plates in this study are excited by an impulse hammer equipped with force transducer. By using a signal analyzer, frequency response can be calculated. Natural frequencies obtained by modal test for these plates are listed in Table 4.
In order to promote validation of the present study, mechanical constants of these glass/epoxy plates are also obtained by investigating tensile test specimens in laboratory (Figure 2). Tensile test machine is attached by a force sensor and also a strain gauge in direction of the applied force. As it is obvious from Figure 5, two other strain gauges are installed on the specimens in direction of 90° and 45° relative to direction of applied force. So, by recording the amount of strains in 3 different directions and also amount of applied force, 4 independent constants are obtained by the tensile test which are listed in Table 3.
(a) Tensile test after failure
(b) Tensile test before failure
7.2. Validation of CPAM Method for Determination of Homogenous Plates (Aluminum) Constants
In order to validate CPAM method described in the previous sections, first, an attempt is made to determine mechanical properties of the discussed aluminum plates in the work of Deobald and Gibson [9]. Deobald and Gibson considered these plates as a transversely isotropic material and tried to determine their four elastic constants by modal analysis and Rayleigh-Ritz method. Four elastic constants have been investigated in their study, while aluminum plates are homogenous and it is enough to determine their two constants; therefore, Hwang and Chang [8] determined these two constants beside determining the four ones and compared their results with those reported by Deobald and Gibson. Results for determination of four elastic constants by using the presented approach and also comparison of them with the average results of Deobald and Gibson and Hwang and Chang’s study are shown in Table 6. Table 7 represents comparison between the results provided in this study for investigation of two elastic constants and the average of corresponding ones, obtained by Hwang and Chang.
It is obvious that considering just two elastic constants results in a more accurate estimation and is less time-consuming than considering four ones. So, if it is guaranteed that experimented material is homogenous, just two elastic constants should be determined and even if four elastic constants are investigated, the results will be reliable. If two constants are investigated, the following restrictions must be also applied to the optimization algorithm; however, the CPAM inverse method with four elastic constants can be usable for both homogenous and orthotropic materials.
In order to assure the accuracy of the CPAM inverse method for determination of parameters of isotropic and orthotropic plates with various boundary conditions and stacking sequence, more cases are investigated by Huang and Chang. Dimensions, density, and layout of considered plates, their mechanical properties, and obtained natural frequencies are presented in Tables 2, 3, and 4, respectively. Using the frequencies given in Table 4, two and four elastic constants are obtained by the present approach for plates A1 and A2; Table 8 shows these results in comparison with the corresponding ones reported by Huang and Chang.
Although four constants determined for both specimens A1 and A2 have a high accuracy like aluminum plates investigated by Deobald and Gibson [9], characterization of two constants is done faster and more accurately than characterization of four constants. Elastic constants of specimen A1 are obtained with higher accuracy than the ones obtained for specimen A2. This is mainly because, for the case of free boundary condition at all edges in comparison with other considered boundary conditions, the natural frequencies predicted from numerical solution are closer to the experiment frequencies. This viewpoint is exactly evident in Table 5. It can also be because implementation of free boundary conditions in practice (experiment) is definitely closer to the theoretical case (finite element modeling in ABAQUS) in comparison with the other conditions such as clamped or supported conditions.
Elastic constants obtained from inverse solution presented in this study have higher accuracy than elastic constants obtained by Huang and Chang [8] and Deobald and Gibson [9]. As it is clear from Table 1, natural frequencies obtained from ABAQUS are closer to the experimental frequencies than those obtained from Rayleigh-Ritz method and ANSYS. Also the physical constraints given in (9) cause the optimization algorithm to generate parameters closer to reality and therefore lead to more accurate elastic constants determined in this study.
7.3. Validation of CPAM Method for Determination of Orthotropic Plates (Composite) Parameters
As it has been already mentioned, carbon/epoxy composite plates studied in Hwang and Chang’s paper [8] with free boundary conditions on their all edges as well as glass/epoxy plates tested in the modal laboratory are considered here for characterization of composite plates. Parameters of these plates are elastic constants and fiber orientation angles. As mentioned before, these parameters are obtained simultaneously and the results of this simultaneous investigation for mechanical constants are presented in Table 9(a) and those for stacking sequence in Table 9(b). Specimen B1 has approximately square shape and B2 is a rectangular plate. These plates are approximately similar to plates C1, C2, C3, and C4 in terms of dimensions. These six specimens, as their length/thickness ratio is about more than 50, are considered as thin plates. Dimensions of D-type specimens are different; so, they are considered as thick plates [8]. As can be seen from Table 9(a), elastic constants determined by the presented method for carbon/epoxy plates are compared with those investigated by Hwang and Chang [8]. These constants are compared with results of static tensile test (Table 3) for glass/epoxy plates and amount of error percentage is also calculated for them.
As seen from Table 9(a) for carbon/epoxy plates, elastic constants obtained for square and rectangular B-type and C-type plates are not so different from each other. So, using such square or rectangular plates does not make much difference in outcomes for determination of these constants. Results for elastic constants obtained from CPAM inverse method show a good agreement with corresponding ones obtained from static tests available in the literature for B1, B2, C1, and C2 [8]. The differences between estimated constants from presented inverse method and those measured from experiment cannot necessarily just be due to error associated with vibration modeling in ABAQUS software or PSO algorithm. Hwang and Chang [8] also showed that such differences can be caused by the sources of noises created by test and measurement equipment. As it is explicitly evident from Table 9(a), results obtained from the presented method are relatively more accurate than those reported by Hwang and Chang [8]. As it is noted previously, this fact is because the natural frequencies predicted by ABAQUS solver are closer to the values measured experimentally than frequencies obtained by ANSYS. Moreover, the appropriate constraints for elastic constants given in (9) make optimization algorithm generate parameters that are closer to reality. For glass/epoxy plates, amount of errors of these constants obtained for C3 and C4 (in comparison with tensile test) shows high accuracy of this method and also experiment. There is no difference between results of C3 and C4 which means thta differences in stacking sequence of plates do not affect accuracy of the present method.
As it was also mentioned in previous sections, four independent elastic constants must be predicted for thin plates, but constant is considered as another independent constant that must be investigated for thick plates D1, D2, D3, and D4 and, therefore, the second part of (6) should be omitted from algorithm code. Because high-frequency modes affect more transverse shear deformation of plates than low modes, the first ten frequencies presented in Table 4 are used in error function in order to determine the elastic constants of specimens D1, D2, D3, and D4.
Stacking sequence identification of plates B1 and B2 is done by considering just five unknown angles. It means that stacking sequence is assumed as in which () must be determined and due to the presence of the repetition pattern of sublaminates in this lay-up configuration, it is enough to identify just these five unknown angles instead of twenty ones. In order to evaluate ability of CPAM method for identification of five angles, these five angles are considered unknown. Similarly, for plates C1, C2, C3, and C4, stacking sequence is considered as such that four angles must be evaluated. In another case of plates C2 and C4, all angles are considered to be known for identification system except 6 angles . Results obtained from this identification for C2 and C4 are indicated by F1 and F2, respectively. Results obtained from this stacking sequence identification for plates B1, B2, C1, C2, C3, C4, F1, and F2 are presented in Table 9(b). It can be seen from Table 9(b) that CPAM method can identify accurately fibers’ angle of the plates with eight and ten layers and provide acceptable results. However, for cases F1 and F2, all angles are identified accurately except one of them. This can be because of difference between experimental and theoretical frequencies. It is clearly noted that, according to the results of Tables 4 and 5, natural frequencies obtained from experiment deviated a little from numerical solution and there is a difference between these two types of frequencies. As it is mentioned previously, it can cause inborn error in the inverse solution system.
Here, five angles are considered unknown for thick plates D1 in addition to four ones for plates D2, D3, and D4 to investigate parameters of them and their results for stacking sequence identification for them are presented in Table 9(b). The way that the number of angles is chosen for both of these plates is similar to the way explained for B-type and C-type plates previously. In addition, symbols G1 and G2 are also used for calculation of six angles for plates D2 and D4 (like symbols F1 and F2 for thin plates C2 and C4). As it is clear from Table 9(b), CPAM method can identify angles of these eight-layer, ten-layer, and six-layer plates with high accuracy.
For thin plates, four independent elastic constants have been considered in the identification process. But, for thick plates D1, D2, D3, and D4, the constant is also considered as fifth independent variable in the identification algorithm. So, the second part of (6) is not used for thick plates. For thin plates, first six natural frequencies are used in the error function (see (7)). Since the frequencies of higher modes have more effects on the transverse shear deformation of plates in comparison with lower modes and also due to being more sensitive to this effect in thicker plates, instead of considering the first six frequencies, the first ten frequencies that are presented in Table 4 are used in construction of the error function; it can be another important reason for high accuracy of G-type cases in comparison with F-type ones.
7.4. Sensitivity Analysis
Study of stability of an inverse problem which results in smooth movement of all parameters to the real response of the problem is an especially important issue that must be considered in inverse methods. The inverse solution may diverge by moving from initial chosen points in the solution space. It means that the difference between measured experimental data (real responses) and those evaluated from direct solution (simulated model) is increased. It is essential that all parameters generated in solution space move towards the real responses when the process continues; under such circumstances, inversion of solution is stable [17]. As it is defined by Hadamard [17], an inverse problem will be termed well-posed if(1)there is a solution for data in solution space,(2)solution in model space (model parameters) is unique,(3)solution remains stable to the small changes in the input data.
If these three tests are not satisfied, the problem will be called ill-posed. If data are sufficient, in other words, data space is larger than model space, the first test will be satisfied. The proposed inverse method has used six and ten experimental natural frequencies for determination of the four and five elastic constants (along with fibers angles), respectively. Also, due to capability of the inverse solution for correct identification of the model space, the first two tests are passed well and the problem is expected to be well-posed for determination of parameters. But the third test of Hadamard is about the stability of inverse solution. In a well-posed problem, a small disturbance (disorder) in the model parameters causes small disturbance in the data space and vice versa. In other words, if a small disturbance in data space creates a large disturbance in model space, an ill-posed problem will be faced. One way to evaluate the stability of inverse solution system is to directly import some disorders to the measured experimental data and recalculate the model parameters (elastic constants and fibers’ angles) by using these noisy data comprising disorders (simulated data) [18].
One of the artificial disorders that may be used by inverse problems is creating Gaussian disorder. In order to create this disorder, a vector of some quasi-random numbers (using a Gaussian distribution with average and standard deviation ) is generated by Box-Muller method. For applying this disorder to the data, the average is set to zero and the standard deviation is obtained from following relation [18]:where is th natural frequency, is total number of measured frequencies for th mode, and quantifies the level of disorder; as an example, means disorder level is . Using a test by different levels of simulated disorder can provide an appropriate criterion for the stability of the inverse method against being polluted with disturbances.
Thus, frequencies polluted with disturbances and obtained by adding Gaussian disorder to the measured experimental values are utilized for simulating experimental data that are measured with noise. If Hadamard tests for simulating data with an appropriate level of disorder are passed well by an inverse method, then this method can be employed everywhere for scientific and practical usage. So, effect of different levels of disorder (1, 5, and 10 percent) is investigated for specimens A2, C1, C3, D1, and D3 and consequently the obtained results are presented in Tables 10, 11, 12, 13, and 14, respectively. Ten groups of frequencies are simulated for each level of disorder and the presented results in these tables are the average of this ten-time operation of inverse solution (solution must run once for each group of simulated frequency to obtain material properties). Percentage errors presented in these tables for different levels of disorder are calculated in comparison with the free disorder values obtained from the presented method (it is not calculated in comparison with reference values for elastic constants).
The results of Tables 10, 11, 12, 13, and 14 show that Young’s modulus has the least sensitivity to disturbances and disorders in the system and relatively remains with good accuracy, even in high levels of the disorder. Young’s modulus is a little sensitive to disorders in the system; however, results obtained for this modulus are again stable and satisfactory. Poisson’s ratio is more sensitive than the other constants to disorders appearing in the system because it has a smaller value than others; so a small change and disorder in the system cause it to be changed with high percentage due to its small amount. This factor is well evident in the results of other constants. As it is mentioned above, which is smaller than has higher percentage of error (and so does modulus in comparison with ). High percentage of error for out-of-plane shear modulus in thick plates D1 and D3 indicates that this material constant shows more sensitivity than other constants to changes in natural frequencies in transverse vibration. As it is mentioned previously, since high-mode frequencies influence more the transverse shear deformation than low-mode frequencies, this sensitivity in thick plates seems reasonable. The fibers’ angles of both thin and thick plates show acceptable resistance against disorders and the obtained results are reasonable. As can be seen from Tables 11, 12, 13, and 14, angles that are closer to the mid-planes (inner stacking sequence) are more sensitive than angles further away from them (outer stacking sequence). The reason is that outer stacking sequence changes more stiffness matrix (and, consequently, natural frequencies) of the plates in comparison with inner stacking sequence. According to the results presented in Tables 10, 11, 12, 13, and 14, because the deviations for the estimated parameters (that are obtained with disorders in the system) are in an acceptable domain and also in the same range of disorders, the proposed inverse method passes the third test of Hadamard successfully and is appropriate to use in practical applications.
8. Conclusion
In this paper, a new combined FEM optimization method is presented to determine the elastic constants and stacking sequence of composite plates concurrently by using natural frequencies measured experimentally. To determine the parameters of aluminum, carbon/epoxy, and glass/epoxy composite plates, this rapid, convergent, and accurate method is employed with applying proper constraints and choosing appropriate parameters. It is shown that the proposed method can predict well the parameters of such materials. Furthermore, implementation of free boundary condition at all edges of the plates in vibrational experiment leads to better results. Although four elastic constants can be investigated for homogenous materials, determination of two elastic constants increases the accuracy and quickness of results. No obvious dependence on the laminate dimensions has been observed for determination of parameters of carbon/epoxy and glass/epoxy composite laminates. By comparing the results obtained through this method with those obtained by similar methods available in the literature, the high accuracy of the proposed method is properly confirmed. This method can predict out-of-plane shear modulus with agreeable accuracy for thick orthotropic plates. Finally, by applying sensitivity analysis with three levels of disorder (1, 5, and 10 percent), the stability of the method is evaluated and, as a result, the proposed inverse solution shows high and good resistance against unwanted disorders.
Conflicts of Interest
The authors declare that they have no conflicts of interest.