Analytical and Mathematical Modeling and Optimization of Fiber Metal Laminates (FMLs) Subjected to Low-Velocity Impact via Combined Response Surface Regression and Zero- One Programming

This paper presents analytical and mathematical modeling and optimization of the dynamic behavior of the fiber metal laminates (FMLs) subjected to low-velocity impact. The deflection to thickness (w/h) ratio has been identified through the governing equations of the plate that are solved using the first-order shear deformation theory as well as the Fourier series method. With the help of a two degrees-of-freedom system, consisting of springs-masses, and the Choi’s linearized Hertzian contact model the interaction between the impactor and the plate is modeled. Thirty-one experiments are conducted on samples of different layer sequences and volume fractions of Al plies in the composite Structures. A reliable fitness function in the form of a strict linear mathematical function constructed. Using an ordinary least square method, response regression coefficients estimated and a zero-one programming technique proposed to optimize the FML plate behavior subjected to any technological or cost restrictions. The results indicated that FML plate behavior is highly affected by layer sequences and volume fractions of Al plies. The results also showed that, embedding Al plies at outer layers of the structure significantly results in a better response of the structure under low-velocity impact, instead of embedding them in the middle or middle and outer layers of the structure.

small mass and the large-mass impact, which is based on the ratio of the impactor mass to the target mass.Vlot (1997) showed that the damaged zone of FMLs after the impact is smaller than the FRPs.He also demonstrated that the energy to create the first crack lamia interfaces between each Aluminum and fiber layer for FMLs is relatively low compared to traditional FRPs.One knows that achieving to the complete successful bounding in interfaces among lamia interfaces is the main technology challenge in production of FML structures.So, the cracks usually initiates among lamia interfaces in FML structures having lower energy relative to crack initiation energy of fiber composite laminates.It is clear that the materials of each two bonded laminas in FMLs are not identical.Caprino et al. (2004) showed the correlation between the maximum force and displacement, residual displacement, and energy can be described by simple empirical laws, provided that the material is far enough from penetration.They also demonstrated that the overall impact force-displacement curve of FMLs under low-velocity impact only depends on the impact energy, rather than mass and speed of the impactor separately.Caprino et al. (2007) presented a mechanistic model with neglecting the macroscopic behavior of the structure (i.e.neglecting local deformation due to indentation, overall deflection, damage initiation and development) to study the low-velocity impact.They also showed that the residual displacement after impact, mainly due to the plastic deformation of aluminum, is linearly dependent on the square root of the impact energy.Atas (2007) has done an experimental investigation to carry out the damage process of FMLs under low-velocity impact.He demonstrated that the plastic deformation and shear fracture in the aluminum layers as well as fiber fracture and delamination in the composite plies were the primary energy-absorbing mechanisms.In addition, the delamination between the aluminum layers and the composite plies should not be ignored.Abdullah and Cantwell (2006) showed the positive effect of the FMLs in comparison with the FRP laminates in high-velocity impact too.They presented that the membrane stretching, plastic deformation and shear fracture in the aluminum layers as well as plastic drawing and delamination in the nonmetal fiber reinforced layers were the primary sources of energy-absorbing mechanisms in FMLs.They investigated that the layer sequence is an important parameter on the perforation resistance of FMLs structures.
It is seen that there is no effort in optimizing the dynamic behavior of FML structures under low-velocity impact.The paper presents an analytical and a mathematical model to study the dynamic behavior of the FMLs under low-velocity impact in detail.In addition, effect of layer sequences and volume fraction of Al layers are studied.Many approaches for multidisciplinary optimization of composite laminates were proposed to save the number of fitness function evaluations.In most of them genetic algorithms have been widely applied to laminate design optimization problems, since they are suitable for integer programming and are more likely to find global optima (Soremekun, 2001).They do not require gradient or sensitivity coefficient evaluations; however, a large number of analyses must be done to obtain the optimal solution.Here, using an experimental design technique in binary space (Montgomery, 2008), we first offered a fitness function to deflecting ratio as a function of layer sequences and then an integer programming technique to seek the optimum solution to decide the layer sequences proposed.In the present approach, any technological or financial restrictions could be accessible.

Governing Equations
In this paper, the developed plate equations by Whitney and Pagano (1970) are used as follow as: (1) Latin American Journal of Solids and Structures 10(2013) 391 -408 u 0 , v 0 and w 0 are the plate displacements in x, y and z directions at the plate mid-plane and and are the shear rotations in the x and y directions.For the specially orthotropic form (B ij = 0, A 16 = A 26 = D 16 = D 26 = 0) results in: (2) k sh is the shear correction factor introduced by Mindlin (1951), usually equals to π 2 /12 and q is the dynamic normal load over the plate and also: ; (3) which ρ 0 shows the density of each layer and ρ is the total density of the plate, I is the moment of inertia and h is the thickness of the plate.
are the reduced in-plane stiffness components and are the reduced transverse shear stiffness components as described by Whitney and Pagano (1970).
Here, a simply supported rectangular plate is chosen to study with the dimensions of a and b, with the following boundary conditions: (4) A FML plate that is in contact with a spherical mass is shown in Fig. 1 (Payeganeh et al., 2010).
Figure 1 Schematic view of the FML impacted by a spherical mass on its center (Payeganeh et al., 2010)

Constitutive Equations
Constitutive equations of stress-strain relationship for a FML are as follows (Reddy, 1997): Latin American Journal of Solids and Structures 10(2013) 391 -408 = ; (5) which {σ} represents the stresses in the principle directions, the matrix {ε} shows the strains in the principle directions, and denotes the reduced stiffness matrices for the FML structure.
The constitutive equation could be determined by considering the force-couple resultants in terms of stresses, using integration of Eq. ( 5) through the plate thickness, which yields (Shokuhfar et al., 2008): where N and S are vectors of forces and M is the vector of moments respectively., , and are the components of extensional and shear, coupling, and bending stiffness matrices respectively.Also and are the midplane and shear strains respectively, κ is the curvatures, and kR sh R is the shear correction factor.

Contact Force Analysis
Fig. 2 (Payeganeh et al., 2010) shows a two-degrees of freedom springs-masses model is used to determine the impact force (Fig. 2).The motion equation is equal to: (7) Figure 2 A two-degrees of freedom springs-masses model (Payeganeh et al., 2010) which the index 1 indicates the plate and the index 2 belongs to the impactor.F is the contact force, m 1 =m p and m 2 =m i represents respectively the mass of the impactor and the FML plate, z 1 and z 2 are respectively the relative displacements of the structure and the impactor masses, Latin American Journal of Solids and Structures 10(2013) 391 -408 K 1 = K bs is the bending-shear stiffness the plate.Using of the Choi's linearized model (Choi, and Lim, 2004), instead of nonlinear Hertzian contact law, the contact force can be obtained as: In the above equations, K 2 = K l represents the linearized contact coefficient in Choi's linearized contact law, F m is the maximum predicted contact force, and K c represents the contact stiffness in the improved Hertzian contact law (Choi, and Lim, 2004), where R 2 is radius of the curvature, ν 2 is the Poisson's ratio, E and E 22 are the elastic modulus of the isotropic impactor and the transverse elastic modulus of the top layer of the FML structure respectively.
Replacing the value of F in Eq. ( 7) with the values presented by Choi's differential equations after some simplifications (Payeganeh et al., 2010): The system of ordinary equations could be solved by the Runge-Kutta method, using of the MATLAB ODE solver.The Runge-Kutta method is very good known numerical method for solution of ordinary differential equations particularly for linear simple set of ODEs (Eq.9).These equations are homogeneous equations with initial value conditions.The initial values of displacements and velocities are as follow: These initial conditions must be considered in the MATLAB ODE solver.

Deflection and Stress-Strain Analysis
The current problem could be converted to a system of ordinary differential equations of secondorder in time for the Fourier coefficients of the transverse deflection (Carvalho, 2000), neglecting the effects of rotary inertia (Mindlin, 1951).Therefore the motion equation of the FML plate subjected to a point load q(x, y, t) is equal to (Payeganeh et al., 2010, andKhalili et al. 2007a): where for a concentrated load located at the point (x c , y c ), P mn (t) are the terms of the Fourier series as: (12) which F(t) is the impact load (See Eq. ( 8)).The impact solution for a rectangular plate with simply supported boundary conditions is assumed to be in the following form (Christoforou, 1991, andKhalili, 2007b): where A mn (t) , B mn (t) and WR mn R(t) are the time dependent coefficients.Using the changes of variables method (Khalili, 2007c), Eq. ( 11) simplifies as follows: (14) where: is the square of the fundamental frequencies of the plate.For m = n = 1, which K 1 is equal to: (16) The value of W mn (t), would be easily calculated based on the Runge-Kutta method of 4 th and 5 th ranks and using a software like MATLAB and its ode 45 solver.Substituting the results in Eqs.
(13), the values of w, and could be calculated.
In many practical situations, two or more variables are related and there is an interest to develop a strict mathematical model to estimate the relationship between response variable respect to the relevant independent variables.For example, in composite materials, the w/h ratio of a composite plate is an unknown function of layer sequences.In order to predict the desired cases or tracking optimum process setting, fitness function must be obtained using regression method (based on a set of experiments).
In general, the response variable y may be related to k predictor variables named as X i (i=1,2,…,k).The surface regression model could be as (Montgomery, 2009): It is called a multiple linear regression model and may be written in matrix notation showed in bold characters as follow as: (18) or in expanded form: (19) where the least square estimator of fitness function is equal to: (20) where notations " -1 " and " T " show the inverse and transpose of the desired matrix, respectively.The significance of fitness function could be evaluated using a test to determine whether a linear relationship exist between the response variable and a subset of the predictor variables.The appropriate statistical hypotheses are equal to: Rejection of H 0 in equation implies that at least one of the predictor variables contributes significantly to the fitness function.The test procedure involves an analysis of variance, ANOVA, portioning of total sum of square due to the model (SS(R)) and the residuals or errors (SSE), as follow as: where the sum of squares may be computed by: The test procedure for H 0 must be computed as:  The accuracy is checked by comparison of force-time relationships obtained by the present solution with those generated using an analytical method of Pierson, and Vaziri (1996) and those obtained from the experimental results of Delfosse et al. (1993).It was shown (Payeganeh, 2010) that a good agreement happens in the results.The effect of number of terms of the Fourier series in the solution for the transverse deflection ratio (deflection to thickness) of the plate is also verified (Payeganeh, 2010).Good convergence of the present series solution is reached by using of only 9 terms, but the full convergence is demonstrated with 100 terms.
Here, the effect of layer sequences and volume fraction of the Al layers on the impact response of FMLs is studied.The FML plate used is cross ply and symmetric.The impactor is a heavy spherical object (impactor mass/plate mass>2).Material and geometrical properties of the FML plate as well as the impactor are presented in Table 2.The structure consists of 10 layers and they are numbered from top to bottom.
Table 2 Geometrical and material properties of the FML plate and the impactor (Khalili et al., 2007b, Malekzadeh et al. 2006, and Krimbalis, 2007).In statistical terms, composite production process could be called as a multivariate problem due to simultaneous dealings of more than one variable.In such cases, process setting and optimizing should be followed after interpreting response variables in terms of gradient or sensitivity coefficient evaluations.Usually a large number of analyses are required to obtain the optimal solution.

Geometrical properties of FML plate
Several variables affect dynamic behavior of FMLs subjected to low-velocity impact.Here, the effects of changing the volume fraction and layer sequences of the Al plies in FML plates are studied in detail.After embedding of some Al plies in 31 simulated random selected layer sequences with some volume fractions k s of the Al plies in the FML plate (here, k s = 0.2, 0.4, 0.6, 0.8 and 1), the w/h (deflection to thickness) ratio calculated through the analytical method described (Eq.( 14)).If Ks=1 then, the FMLs plate convert to a laminated plate with identical Al layers.In other hand, the plate changes to identical isotropic layered plate without any composite laminate.
Table 3 shows arrangement of the layer sequences and calculated results of the w/h ratio of the structure.
In order to model layer sequences, we have used five binary variables of X 1 , X 2 , X 3 , X 4 and X 5 .
These variables just get 1 if symmetrical FRP layers numbered as (1,10), (2,9), (3,8), (4,7) or (5,6) replaced with Al plies.For example, any structure embedded with Al plies in its outer symmetrical layers (1, 10) could be modeled by X 1 =1, X 2 =0, X 3 =0, X 4 =0 and X 5 =0.Embedding Al plies in the both inner layers (4, 7) and (5, 6) could be modeled as X 1 =0, X 2 =0, X 3 =0, X 4 =1 and X 5 =1.Using such arrangement, Table 3 could be explained in another technical terms through Table 4.  Figure 3 graphically shows the effect of changing layer sequences on w/h ratio.One sees that the biggest amount derived from the 5 th , 15 th , 25 th and 30 th sequence numbers and the lowest ones belong to 1 st , 6 th , 16 th, 26 th and the 31 st ones.This means that the worst way of embedding the Al layers is embedding them in the middle layers and the best way is embedding them in the outer ones.Using the outer layer sequences (1 st , 6 th , 16 th, 26 th and the 31 st ones) result in an improve of the stiffness of the outer layers of the structure, the layers that should stand the most of amount of the impact load.
Readers could easily reveal that such a composite plate having no Al plies (it means FRPs) has a w/h ratio amount of 0.971.In order to investigate fitness function of w/h ratio on the five binary layer sequences variables, a linear response surface regression method applied on the calculated data sets.Tables 5 and 6 show two standard output of ''SPSS'' statistical software nominated as parameters estimation and analysis of variance (ANOVA) tables which show a reliable response surface model fitted to the measurements as: (28) where shows the errors term.Standard deviation and percent error were investigated to validate the experiments.Errors between predicted and actual values were calculated according to:  Minimization of the w/h ratio could gain by focusing on the signs and coefficients of Eq. ( 28) easily.If the Al plies are embedded in the outer layers of the structure, w/h ratio will be reduced more than embedding them in the inner ones.This means that embedding the Al plies in layers (1, 10) or (2, 9) results in the minimum value of the w/h ratio of the structure.Whereas embedding them in layers (3,8) or (4,7) does not lead to the same results.To compare the value of embedding Al plies in each layer, coefficient of fitness function could be employed.From Eq. ( 28), one also sees that the effective coefficient of X 1 (it means, embedding the Al plies in layers (1,10)) and the effective coefficient of X 2 (in layers (2, 9)) are both -0.17.Whereas the effective coefficient of X 3 (embedding the Al plies in layers (3, 8)) equals to -0.081 (about one half of X 1 or X 2 ), and the effective coefficient of X 4 (in layers (4, 7)) equals to -0.090 (more than one half of X 1 or X 2 ).The more interesting result is that the effective coefficient of X 5 (embedding the Al plies in layers (5, 6)) is equal to +0.034.The positive sign in Eq. ( 28) means that embedding Al plies in layers (5, 6) of the structure, have a full negative effect on dynamic response of these structures.
Fitness function also shows that the minimum of achieves while X 1 =1, X 2 =1, X 3 =1, X 4 =0 and X 5 =0.So in order to achieve minimum w/h ratio, layers could be embedded by Al plies in all layers except the four inner layers (it means layers (4, 7) and (5, 6)).It is also seen that embedding Al plies in outer layers of a FRP plate subjected to low-velocity impact results in a better response, instead of embedding them in the middle or middle and outer layers of the structure.
The impact event usually is the contact local event.In the contact phenomena the outer layers play the main role.Indentation usually occurs on outer contacted layer.Therefore, in order to increase the ductility of contacted area, using of Al layers in outer layer is sufficient.Because, in these plies we have minimum w/h ratio with respect to the other ply sequence.It is also seen that embedding Al plies in outer layers of a FRP plate subjected to low-velocity impact results in a better response, instead of embedding them in the middle or middle and outer layers of the structure.If there are any restriction due to cost or any other constraints, Table 7 could conduct engineers to the optimum layer sequences subjected to their desired volume fractions.Therefore, embedding the Al plies only in 4 layers seems to be the best choice.Embedding the Al plies in 6 layers, only reduces the w/h ratio about 4 percent, whereas the weight of the structures increases about 9 percent that is not an advantage at least in aerospace industries.
Generally, by solving the following integer programming problem via any operational research computer packages, the optimum layer sequences could be attained subjected to different layers embedding cost and threshold on maximum layers enrichment cost: (30) where i C and i W show the cost and weights of the layer i, respectively.In addition, TC shows the total amount of cost and W explains total allowed weights of composite plates.Interested readers could easily add or remove any required restrictions (such as technological restrictions e.g., time consuming) to follow zero-one programming and find optimum solution via tradeoff among all technological, economical and available resource limitations.

CONCLUSION
The presented work presents analytical and mathematical modeling and optimization of the behavior of the fiber-metal laminates (FML's) to low-velocity impact.Thirty-one fully random-Latin American Journal of Solids and Structures 10(2013) 391 -408 ized simulated experiments are conducted on samples of varying layer sequences and deflection ratio, w/h, has been identified through the governing equations of the plate that are solved using the first-order shear deformation theory as well as the Fourier series method.The interaction between the impactor and the plate is modeled with the help of a two degrees-of-freedom system, consisting of springs-masses, and the Choi's linearized Hertzian contact model.
The proposed fitness function showed that FMLs are highly affected by layer sequences.The results show that, in general, embedding the Al plies in outer layers of the structure results in smaller w/h ratios respect to embedding them in inner layers.One also sees that embedding about 40 percent of volume fraction of traditional composites with Al plies leads to the best response of the structures under low-velocity impact.The place of embedding the Al plies is important too.It is shown that embedding the Al plies instead of the outer layers of the FRP structure (here, layers (1, 10) and (2, 9)) results in the best response of the structure under lowvelocity impact.Fitness function also makes known that minimum deflection ratio achieves while layers embedded by Al plies in all layers except to the four inner layers.This is a brilliant result not also for the researchers who works on this subject, but also for the manufacturers especially in aerospace industries.
Finally, the proposed general set of standard 0-1 programming based on five decision variables on one minimization object function subject to two less than or equal constraints could conduct process engineers to carry out their design in right manners.
(10) where the L ij coefficients are introduced by Payeganeh et al. (2010) as: Latin American Journal of Solids and Structures 10(2013) 391 -408 (11) addition, to reject H 0 if calculated F exceeds .The test is usually summarized in an ANOVA

Figure 3
Figure 3 Amount of calculated deflection to w/h ratio versus the different layer sequences.
within 5 percent beside their normal distribution; indicated that process optimization by regression were capable and reliable to optimize w/h ratio.By substituting the measured amounts in the model, some errors were achieved.The realized regression model shows normal random distributed errors with a mean about zero and a low amount of constant standard deviation.ofSolids and Structures 10(2013) 391 -408

Table 1
ANOVA for significance of multiple regressions table such as Table 1.Latin American Journal of Solidsand Structures 10(2013) 391 -408

Table 3
Calculated w/h ratio in 31 set of layer sequences with different volume fractions ks of Al plies in the FML structure Latin American Journal of Solids andStructures 10(2013) 391 -408

Table 5
Regression parameter estimation based on 31 described samples

Table 7
Optimum layer sequences based on the desired volume fractions of Al plies