- Email: [email protected]
An inverse solution for low-velocity impact in composite plates
Computers and Structures 79 (2001) 2607±2619
www.elsevier.com/locate/compstruc
An inverse solution for low-velocity impact in composite plates Andreas P. Christoforou *, Abdallah A. Elsharkawy, Lot® H. Guedouar Department of Mechanical and Industrial Engineering, College of Engineering and Petroleum, Kuwait University, P.O. Box 5969, Safat 13060, Kuwait Received 2 June 2000; accepted 28 June 2001
Abstract An inverse model for estimating the mass and velocity of the impactor in low-velocity impact of composite plates is presented. In the model, a least-squares optimization technique is used to optimally reconstruct the force±time history by comparing direct model simulations of the impact response with experimental measurements. Guidelines for better initial guesses and faster direct models that are based on an impact characterization procedure are provided to the inverse algorithm. Three sets of experimental data covering a wide range of impact response are used to validate the model with excellent results. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Impact; Inverse analysis; Composites
1. Introduction A major problem in composite structures is their low resistance to damage due to impact of foreign bodies. A key step in understanding the impact phenomenon is to be able to estimate the response and impact force through analytical models and/or experiments. It is known that the impact response depends on the geometry and material of the structure, the boundary conditions, the mass, the geometry, and velocity of the impactor. Usually in direct analyses and experimental measurements, all these parameters are known. In real cases however, when impact occurs, the mass, the geometry, and velocity of the impactor are not known and an inverse model is needed to accurately estimate them. This requires that the impact response of the structure is known. By estimating the characteristics of the impactor, the impact response (i.e., force, displacement, and strain histories) of the structure can be reconstructed. This is very important since for a given structure, the state of
* Corresponding author. Tel.: +965-481-1188, ext. 5790; fax: +965-481-7131. E-mail address: [email protected] (A.P. Christoforou).
damage depends on the impact characteristics. Furthermore, the identi®cation, estimation, and reconstruction of the response are essential for controlling the severity of impact and health monitoring of the structure. For these reasons, the identi®cation of the impact force has received much attention in recent years. In early papers, Doyle [1,2] used a frequency-domain deconvolution technique to estimate the impact force on beams and plates from strain gauge measurements. The fast Fourier transform (FFT) was used to convert the measured strain history into its frequency components, which were used to obtain the frequency components of the impact force. Then the inverse fast Fourier transform (IFFT) was implemented to reconstruct the impact force history in the time-domain. Chang and Sun [3], on the other hand, used experimentally generated Green's functions and time-domain deconvolutions to estimate the impact force history from strain measurements. A least-square ®t method was used in their numerical scheme. However, the amplitudes of the Green's functions had to be scaled to account for dierent impactor characteristics. Hollandsworth and Busby [4] used acceleration measurements to identify the impact force on aluminum beams by solving a dynamic inverse problem. Wu et al. [5] presented two methods for the purpose of detecting the impact force history. In the ®rst method,
0045-7949/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 1 ) 0 0 1 1 3 - 4
2608
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
the classical plate theory, and the normal-mode superposition method using Bessel functions were used to construct Green's functions in the time-domain. A constrained optimization technique was used to reconstruct the impact force history from measured strain histories. This method was limited to axisymmetric impact problems because of the restriction of the Green's functions. In the second method, the impact force history was reconstructed using the gradient projection method by comparing the response to recorded force and strain histories of a controlled impact. Tsai et al. [6] also used a Green's function approach to estimate impact force from experimental data on sandwich panels. With the advent of smart structure technology, Choi and Chang [7] presented an impact load identi®cation method to estimate the location and history of an impact force using piezoelectric sensor measurements. Their identi®cation process consisted of a system model and a response comparator, employing a quasi-Newton procedure to estimate the location of impact, and the leastsquare ®t method to predict the force history. Lihua and Baoqi [8] used neural networks and piezoelectric sensors to identify the location and amplitude of a static force on a square plate. The scope of this work is to provide an easy and reliable methodology to estimate the mass and velocity of a foreign object impacting a structure. The inverse model presented below, employs the direct model developed by Christoforou and Yigit [9,10] and a least square optimization technique to reconstruct the impact force history from experimental data. It is worth noting that in this study, the impactor was taken to be spherical in shape with ®xed diameter but not a sphere as shown in Fig. 1. The problems that are considered here are low velocity impact events, e.g., tools dropped on ¯exible composite structures. Usually the range of impactor diameters is between 10 and 25 mm. In this range, which is the case in the present study, the contact is not visible and has no major in¯uence on the impact response.
Previous studies showed that the relative stiness of the structure (which is a measure of ¯exibility), and the mass ratio of the impactor to the structure signi®cantly aect the impact response [9,10]. The in¯uence of the diameter of the impactor will be in the contact zone, which is taken care by the contact law (especially in impacts of rigid bodies or very sti structures). For smaller diameters the impact will be visible in the form of a damage in the contact zone, and identi®cation is not necessary. Therefore, the mass of the impactor rather than its diameter was considered to be one of the parameters to be estimated by the proposed inverse analysis. Furthermore, the present study is limited to the case of small displacements and with no damage due to impact. However, the proposed scheme is simple, fast and accurate in estimating the impactor characteristics as demonstrated in the experimental results used to validate the model. Certainly, and this is beyond the scope of this initial work, the direct model may be improved by including not only geometrical non-linear eects but also eects due to structural damage, which is more important in impact of composite structures. Also, the scheme may be applied not only to force measurements, but also to other measurements such as from strain gauges or piezoelectric sensors. 2. Theory 2.1. Direct impact model The equations of motion of a rectangular, specially orthotropic, symmetrically laminated, elastic and shear deformable composite plate, subject to lateral dynamic loading, are given as D11
o2 wy o2 wx o2 wx D D D 66 12 66 ox2 oy 2 ox oy ow qh3 o2 wx kA55 wx ox 12 ot2
Fig. 1. Geometry of the impact problem.
1
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
o2 wy o2 wy o2 wx D66 2 D22 2 D12 D66 ox oy ox oy ow qh3 o2 wy kA44 wy 12 ot2 oy kA55
owx o2 w 2 ox ox
qh
kA44
owy o2 w 2 oy oy
o2 w ot2
1 1 m2i 1 E Ez Ei 2
p x; y; t 3
where Dij and Aij are the bending and shear stinesses de®ned as usual, h is the plate thickness, q is the material density, w is the transverse de¯ection, wx and wy are the shear rotations, p is the lateral load per unit area, k is the Mindlin shear correction factor, x and y are the space variables, and t is the time. Assuming a concentrated impact loading, the lateral load per unit area can be written as p x; y; t F t d x
xc d y
yc
4
where F t is the impact force, d is the Dirac-delta function, and xc ; yc is the contact point. In the case of lateral impact of the plate by a foreign object, the motion of the impactor is described by m wi
F t
5
where m and wi are the mass and the displacement of the impactor, respectively. The impact force can be expressed as follows: F t Fc
if Fc P 0
F t 0 if Fc < 0
6
where Fc is to be obtained from an appropriate contact model. In this study a linear version of the elasto-plastic contact model proposed by Yigit and Christoforou [11,12] is used, and is given as Fc a Ky a
7
where a is the local indentation which is the relative approach between the impactor and the target of the contact point given as a t wi t
w xc ; yc ; t
Ky is the contact stiness given as p Ky 1:5Kh ay
2609
8
9
where Kh is the Hertzian stiness for the case of spherical body of radius R, in contact with a ¯at surface, given as p Kh 43 RE 10 where R is the radius of the spherical impactor and E is given as
11
where mi and Ei are the Poisson's ratio and Young's modulus for the impactor, respectively, Ez is the transverse Young's modulus of the laminate, which is usually taken to be equal to that of the ply, i.e., Ez E22 . The onset of ``yielding'' or damage was obtained by combining the classical Hertzian analysis and the elastic± plastic contact theory. By assuming the material is elastic before damage and perfectly plastic after damage occurs, the critical indentation ay was obtained using the maximum shear failure criterion and is given as [11] ay
0:68Sy2 p2 R E2
12
where Sy is the yield strength of the softer material for metals, and Sy 2Su , where Su is the shear strength of the laminate for the ®ber composites. This paper is concerned with a simply supported rectangular composite plate of uniform thickness with dimensions a and b, impacted by a spherical impactor as shown in Fig. 1. The boundary conditions in this case can be expressed as follows: w
owx 0 at x 0; a ox
13a
w
owy 0 at y 0; b ox
13b
The initial conditions of the impact problem are: w x; y; 0 0; w_ i 0 v0
w_ x; y; 0 0;
wi 0 0;
and 14
where v0 is the initial velocity of the impactor. The equations of motion are discretized by a standard modal analysis procedure, and numerically solved to obtain the response. 2.2. Inverse impact model The direct impact analysis shows that the impact response depends on the properties of the structure and the impactor. Furthermore, for a given structure and impact energy, the mass of the impactor is the parameter that governs the impact response (i.e., from locally dominated to globally dominated response, [11,12]). In direct simulations of the response (either numerical or experimental), all impact parameters are known. In the inverse analysis, on the other hand, the response is known usually from experimental data, with some of the impact parameters unknown. In this study, the unknown parameters to be estimated are the mass and the velocity of the impactor, and the predicted values will be denoted m and v0 , respectively. All other plate and
2610
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
Fig. 2. Flowchart of the proposed inverse impact algorithm.
impactor parameters are known. The inverse problem is solved by an optimization technique, where the experimental impact response is compared to the response predicted by the direct solution for an m; v0 pair. The optimal solution is achieved by minimizing the following objective function: f m; v0
N X F m; v0 ; ti i1
subject to
Fe ti 2
15
0 < m 6 mL
16
0 < v0 6 vL
17
where N is the total number of readings, F m; v0 ; ti and Fe ti are the predicted and experimental force±time histories, respectively, ti is the time at ith reading, and mL and vL are limiting values for the mass and velocity of the impactor, respectively. In this study, a M A T L A B routine (version 5.2) called ``leastsq'' which performs least squares optimization was
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
used to obtain the optimal solution. In this routine, both the Gauss±Newton method and the Levenberg±Marquart method are included as options. Basically the dierence between these two methods is the way of the search direction for the optimal solution. For the results presented in this paper Gauss±Newton method was used. Starting from an initial guess for an m; v0 pair, the optimizer will progress towards the solution by solving the direct impact problem to obtain F m; v0 ; ti . The design variables m and v0 are varied using the Gauss±Newton optimization algorithm. This process is guided by objective function gradient evaluations in the design space, and will continue until the objective function reaches a minimum value based on the speci®ed convergence criterion. Fig. 2 shows the ¯owchart for the proposed inverse algorithm. For practical reasons, and without loss of generality, it should be noted here, that the range of the mass and velocity of the impactor is restricted to the range of low velocity impact. A simple adjustment is implemented to force the optimizer to reverse the search direction when any one of the design variables crosses the boundary of the allowable space. This is achieved by deliberately assigning a value to the objective function higher than its previous value.
3. Results and discussion As mentioned earlier, the design parameters of the identi®cation problem are the mass m and the velocity v0 of the impactor. Although the eect of these two parameters on the impact response is well known, at this point of presentation, it is important to show these effects, with the following set of simulations. The material properties, geometry of the composite plate, and the impact conditions are given in Table 1. Fig. 3 shows the eect of the mass of the impactor on the impact force-time response. The velocity v0 of the
Table 1 Properties of the composite plate and the impact parameters Plate 0=90=0=90=0s , T300/934 carbon-epoxy Plate size: a 200 mm; b 200 mm Plate thickness: h 2:69 mm (0.269 mm/layer) E11 120 GPa; E22 7:9 GPa; G12 G23 5:5 GPa m12 m23 0:3; q 1580 kg/m3 Impactor 12.7 mm diameter spherical; q 7960 kg/m3 , steel Contact parameters Kh 8:89 108 N/m1:5 ; ay 2:5 10 Ky 6:65 106 N/m
5
m
2611
impactor was held ®xed at 3 m/s. For small mass m 0:005 kg the duration of the impact is very short (see Fig. 3a). As the mass of the impactor increases, multiple impacts occur, and the duration and the amplitude of the impact force increase. Fig. 4 shows the eect of the velocity v0 of the impactor on the force-time response. The mass of the impactor was held ®xed at 0.5 kg. The results show that the impactor velocity has signi®cant in¯uence on the impact force amplitude. However, it has no in¯uence on the duration of the impact. From the results presented in Figs. 3 and 4, it is clear that the impact response is very sensitive to the mass and velocity of the impactor. Therefore, it will be easy to predict these variables accurately using the inverse technique. In order to validate the inverse model proposed in this paper, two experimental results reported by Pierson and Vaziri [13], were used as inputs to predict the mass and velocity of the impactor. In the experiments, simply supported 45=90= 45=03 s T800/TH/3900-2 CFRP plates with nominal dimensions of 127 76:2 4:65 mm3 were used. In the ®rst experiment, the exact mass and the velocity of the impactor were 6.15 kg, and 1.76 m/s, respectively. Fig. 5 shows a comparison between the impact force responses from the experimental data, the direct solution based on the exact impactor mass and velocity values, and the inverse model prediction. As can be seen there is a small discrepancy between the direct solution and the experiment. As expected, this discrepancy in¯uences the inverse model prediction, and the estimated mass and velocity of the impactor were 5.41 kg and 1.81 m/s, respectively. The eect of the initial guess on the inverse solution was also investigated. It was found that the proposed inverse algorithm is independent of the initial guess. As it will be discussed later, a good initial guess, however, could speed up convergence. Fig. 6 displays the convergence history of the inverse model for two initial guesses, (m 1 kg and v0 5 m=s) and (m 0:5 kg and v0 3 m=s). As can be seen only 30 iterations were needed for convergence to the optimal solution. In the second experiment, the mass and the velocity of the impactor were 0.314 kg, and 7.7 m/s, respectively. Fig. 7 shows a comparison between the impact force responses from the experimental data, the direct solution based on the exact values of the impactor mass and velocity, and the inverse model prediction. The estimated mass and velocity of the impactor were 0.28 kg and 8.098 m/s, respectively. Fig. 8 displays the convergence history of the inverse model for two initial guesses, (m 0:5 kg and v0 3 m=s) and (m 1 kg and v0 5 m=s). These results con®rm further the conclusions drawn from Fig. 6, i.e., the inverse solution is unique and independent of the initial guess. There was no local minimum solution observed in all cases which is in agreement with the physics of the impact problem.
2612
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
Fig. 3. Eect of the mass m of the impactor on the impact force±time response F t.
Again only 30 iterations were needed for convergence to the optimal solution. At this point, it is useful to explore ways for providing guidelines for better initial guesses and/or computationally faster direct solutions in the inverse algorithm. To this end, the characterization procedure introduced by Christoforou and Yigit [9,10] can be used.
In their procedure, it was found that the non-dimensional maximum impact force F max is governed by two non-dimensional parameters, the relative stiness k, and the loss factor f. These parameters are given as F max
v0
F pmax mKy
18
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
2613
Fig. 4. Eect of the velocity v0 of the impactor on the impact force±time response F t.
k
f
Kst Ky 1 16
where
mKy qhD
p D11 D22
19
D 12 D12 2D66
20
and Kst is the global static stiness of the plate. Applying the procedure, the relative stiness of the plate used in the experiments was calculated to be k 0:258. The ®rst experiment is considered as a heavy mass impact (m 6.15 kg), and f and F max were calculated to be 10.61 and 0.45, respectively. Using the estimated
1=2
21
2614
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
Fig. 5. Comparison between the predicted and the actual impact force histories for the ®rst set of experimental data [13].
mass and velocity values, f and F max were calculated to be 9.95, and 0.47, respectively. The second experiment is considered as a medium mass impact (m 0.314 kg) and f and F max were calculated to be 2.4 and 0.513, respectively. The predicted f and F max were 2.26 and 0.515, respectively. Both cases are displayed on a characterization curve as shown in Fig. 9. As can be seen, the agreement between the experiments and inverse prediction is excellent. Also, the ®rst case is located in the so-called fully global region (quasi-static), where the response can be predicted by using a simple one-degree of freedom system. The second case is located in the socalled transition region in which the response is dynamic in nature. Fig. 9 shows the functional relationship between F max , f and k. In the inverse solution f and F max are not
Fig. 6. Convergence history of the inverse model for two initial guesses, for the ®rst set of experimental data [13].
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
2615
known because they depend on the mass and velocity of the impactor. The only parameter that is known is the relative stiness k. Examination of the ®gure reveals that, for practical purposes and especially for sti structures, many impact situations may be considered fully global for a range of values of f. In fully global situations, the response is a function of k only, and is given in a normalized form as
F s
Fig. 7. Comparison between the predicted and the actual impact force histories for the second set of experimental data [13].
r! r ! k k sin s k1 k1
22
In this case F max and the non-dimensional response frequency xr can be written as
Fig. 8. Convergence history of the inverse model for two initial guesses, for the second set of experimental data [13].
2616
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
Fig. 9. Characterization curve for the impact response. Fig. 10. Comparison between the predicted and the actual impact force histories for the third set of experimental data [2].
r k k1
23
r p=tr k xr p k 1 Ky =m
24
F max
F 2fxn F_ x2n F 0;
where tr is the contact duration. Therefore, only Fmax and tr , which can be determined from the experimental response, are needed to estimate the mass and velocity of the impactor. Applying Eqs. (23) and (24) to the ®rst experimental set, the mass and velocity of the impactor were estimated to be 5.45 kg and 1.86 m/s, respectively, which compare well with the full inverse model prediction. Although the second experimental set is dynamic in nature, the same procedure can be applied to obtain a good initial guess. In this case the mass and velocity of the impactor were estimated to be 0.324 kg and 8.59 m/s, respectively. In many cases, this kind of estimations may be adequate. If not, they may serve as initial guesses to the full inverse model. As a rule of thumb, in most impacts having time duration in the milliseconds region, the quasi-static approach can be used with a reasonable accuracy. In the transition region, however, the quasistatic model can provide a good initial estimate. Further study of the ®gure, in the in®nite plate and the local regions, the impact response is a function of f only. In the local region, there are only single short duration impacts (microseconds) and the non-dimensional impact force can be expressed as F s sin s
25
This means that F max 1, and xr 1. As before, the mass and velocity of the impactor can be estimated easily. In the case of impact of an in®nite plate, the response is governed by
xn 1
26
The solution of Eq. (26) will be the direct solution in the inverse algorithm. To validate this model, an experimental set reported by Doyle [2] was used. In the experiment, a seven layer [0/90] Schotchply glass/epoxy composite plate with nominal dimensions of 406 406 1:78 mm3 was used. The mass and velocity of the impactor were 8.537 g and 2.3 m/s, respectively. Fig. 10 shows a comparison between the impact force responses from the experimental data, the direct solution based on the exact values of the impactor mass and velocity, and the inverse model prediction using Eq. (26). The estimated mass and velocity of the impactor were 9.9 g and 2.1 m/s, respectively. The values of f and F max were calculated to be 2.085 and 0.22, respectively for the experimental data, and 2.25 and 0.197, respectively for the inverse model prediction. This set is also displayed on the characterization curves of Fig. 9. As can be seen, the agreement is excellent for both the values and the type of response.
4. Conclusions An inverse model for the low-velocity-impact of composite plates has been presented. An ecient direct solution was utilized in a least squares optimization technique to predict the mass and the velocity of the impactor. A M A T L A B routine called ``Leastsq'' which performs least squares optimization was used to obtain the optimal solution. Three sets of experimental data were used to validate the accuracy of the proposed inverse algorithm. The results showed that the convergence was fast and independent of the initial guess. Also, using a characterization procedure indicating the type of
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
response, excellent guidelines can be provided for better initial guesses, and for simpler and computationally faster direct solutions in the proposed inverse algorithm.
2617
Load experimental data ®le consisting of 2 columns.The ®rst column corresponds to the experimental time data points and the second column represents the corresponding measured impact force load piers1.dat
Appendix A. tion
MATLAB
computer code for inverse solu-
This program solves the inverse impact problem to predict the the impactor mass and velocity from experimentally measured impact force time history. Variables de®nitions: time_exp experimental time data points force_exp experimental force measurements ndata_exp number of experimentally measured data points t®nal_exp time of last experimentally measured data point fact1, fact2 parameters to control the least square search steps error_last last computed error value. Used to force the optimizer to reverse the search direction when the boundaries of the allowable space are reached itr iteration counter ntimes number of times the size of the time vector is enlarged in order to accommodate the largest expected impact duration that may be considered by the optimizer conv_history a matrix to record the convergence history vel0, mass0 initial guesses for the velocity and mass, respectively initail_guess scaled initial guess vector options vector of parameters to control the optimization routines (see Matlab help). sol predicted optimum solution vector (scaled) v0star, mstar predicted optimum mass and velocity, respectively vel_mass Optimum solution vector time, force time and force vectors resulting from the direct solution of the impact problem using the predicted velocity and mass Note: The ®les `impact.exe' and `impact_inv.m' have to be available in the default directory. Clear the workspace and close all ®gures clear; close all; Declaration of global variables global time_exp force_exp t®nal_exp ndata_exp fact1 fact2 error_last global itr ntimes conv_history
ntimes 4; Assign the experimental data points to their corresponding vectors of time and force time_exp piers 1(:,1); force_exp piers 1(:,2); Plot the experimental impact response ®gure(1); plot(time_exp, force_exp,`+'), hold on determine ndata_exp and t®nal_exp ndata_exp length (time_exp); t®nal_exp time_exp (ndata_exp); Initialize the convergence history matrix conv_history [ ] Assign appropriate values to fact1 and fact2 fact1 0.000001; fact2 0.00001; Assign and scale the initial guesses for the velocity and mass vel0 3; mass0 .5; initial_guess [vel0*fact1,mass0*fact2]; Assign the optimization parameters options(2) 1e-7; options(3) 1e-7; options(5) 1; options(14) 200; Initialize the iteration counter itr 0; Solve the optimization problem using the least squares `leastsq' Matlab routine.`impact_inv' is a Matlab function (see listing below) that evaluates the objective function for any given set of the parameters to be predicted, i.e.the velocity and mass.It is used by the least squares optimization routine. [sol, options] leastsq (`impact_inv', initial_guess, options); Compute the predicted optimum velocity and mass v0star sol(1)/fact1; mstar sol(2)/fact2; fprintf(`Predicted velocity and mass: %g %g n', v0star, mstar) vel_mass [v0star mstar];
2618
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
Save the predicted parameters into a data ®le named `param.dat'.This ®le will be read by the executable Fortran code `impact.exe' designed to solve the direct impact problem save param.dat vel_mass -ascii Solve impact problem for the predicted velocity and mass using the executable Fortran code `impact.exe' in order to reconstruct the corresponding force time history.Results will be available in the ®le `Result.dat' !impact Load results load result.dat Assign the results to their corresponding time and force vetors time result (:,1); force result (:,2); Plot the predicted impact response plot(time,force,`g'), hold on Save the convergence history into a data ®le named `history.dat' save history.dat conv_history -ascii function err impact_inv(params): `impact_inv' is a Matlab function that evaluates the objective function (norm of the force error relative to the experimentally measured impact force response) for any given set of the parameters to be predicted, i.e.the velocity and mass.It is used by the least squares optimization routine. Declaration of global variables global time_exp force_exp t®nal_exp ndata_exp fact1 fact2 error_last global itr ntimes conv_history Calculate the velocity and mass velocity params (1)/fact1; mass params (2)/fact2; Check for out of range parameters and take appropriate action if velocity<1 velocity>10 mass<0.00001 mass>10 err error_last*1.001;return end Solve the direct impact problem vel_mass [velocity mass]; save param.dat vel_mass -ascii !impact load result.dat
time result(:,1); force result(:,2); n length(time); t®nal time(n); Interpolate the predicted force response at the experimental time data points.Then adjust the sizes of the experimental and predicted force and time vectors so that they have the same size. if t®nalt®nal); it it (1); force_inter interp1(time,force,time_exp(1:it-1)); force_inter_ mod [force_inter;zeros([ndata_exp-it+1, 1])]; force_exp_mod force_exp; elseif t®nal>t®nal_exp it ®nd (>t®nal_exp); it it (1); force_inter_mod interp 1(time,force,time_exp); force(it:n); force_exp_mod [force_exp;zeros([n-it+1,1])]; else force_inter_mod [interp1(time,force,time_exp)]; force_exp_mod force_exp; end Compute the average error vector err [force_inter_mod-force_exp_mod]/ndata_exp; Enlarge the average error vector err [err;zeros([ndata_exp*ntimes-length(err),1])]; Update the last error value error_last err; Compute the second norm error err_norm sqrt (sum(err.^2)); Update the convergence history matrix conv_history [conv_history; [itr velocity mass err_norm]]; itr itr+1;
References [1] Doyle JF. Further developments in determining the dynamic contact law. Exp Mech 1984;24:265±70. [2] Doyle JF. Experimentally determining the contact Force during the transverse impact of an orthotropic plate. J Sound Vibrat 1987;118(3):441±8. [3] Chang C, Sun CT. Determining transverse impact force on a composite laminate by signal deconvolution. Exp Mech 1989;29:414±9. [4] Hollandsworth PE, Busby HR. Impact force identi®cation using the general inverse technique. Int J Impact Engng 1989;8(4):315±22. [5] Wu E, Tsai TD, Yen CS. Two methods for determining impact-force history on elastic plates. Exp Mech 1995;35: 11±8. [6] Tsai CZ, Wu E, Luo BH. Forward and inverse analysis for impact on sandwich panels. AIAA J 1998;36(11): 2130±6.
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619 [7] Choi K, Chang FK. Identi®cation of impact force and location using distributed sensors. AIAA J 1996;34(1):136± 42. [8] Lihua S, Boaqi T. Identi®cation of external force acting on a plate by using neural networks. In: Chang F-K, editor. Proceedings of the International Workshop on Structural Health Monitoring, Stanford University, Stanford, CA, 18±28 September, 1997. p. 229±37. [9] Christoforou AP, Yigit AS. Characterization of impact in composite plates. Compos Struct 1998;43:15±24.
2619
[10] Christoforou AP, Yigit AS. Eect of ¯exibility on low velocity impact response. J Sound Vibrat 1998;217(3):563±78. [11] Yigit AS, Christoforou AP. On the impact of a spherical indenter and an elastic±plastic transversely isotropic halfspace. Compos Engng 1994;4:1143±52. [12] Yigit AS, Christoforou AP. Impact dynamics of composite beams. Compos Struct 1995;32:187±95. [13] Pierson MO, Vaziri R. Analytical solution for low-velocity impact response of composite plates. AIAA J 1996;34: 1633±40.
www.elsevier.com/locate/compstruc
An inverse solution for low-velocity impact in composite plates Andreas P. Christoforou *, Abdallah A. Elsharkawy, Lot® H. Guedouar Department of Mechanical and Industrial Engineering, College of Engineering and Petroleum, Kuwait University, P.O. Box 5969, Safat 13060, Kuwait Received 2 June 2000; accepted 28 June 2001
Abstract An inverse model for estimating the mass and velocity of the impactor in low-velocity impact of composite plates is presented. In the model, a least-squares optimization technique is used to optimally reconstruct the force±time history by comparing direct model simulations of the impact response with experimental measurements. Guidelines for better initial guesses and faster direct models that are based on an impact characterization procedure are provided to the inverse algorithm. Three sets of experimental data covering a wide range of impact response are used to validate the model with excellent results. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Impact; Inverse analysis; Composites
1. Introduction A major problem in composite structures is their low resistance to damage due to impact of foreign bodies. A key step in understanding the impact phenomenon is to be able to estimate the response and impact force through analytical models and/or experiments. It is known that the impact response depends on the geometry and material of the structure, the boundary conditions, the mass, the geometry, and velocity of the impactor. Usually in direct analyses and experimental measurements, all these parameters are known. In real cases however, when impact occurs, the mass, the geometry, and velocity of the impactor are not known and an inverse model is needed to accurately estimate them. This requires that the impact response of the structure is known. By estimating the characteristics of the impactor, the impact response (i.e., force, displacement, and strain histories) of the structure can be reconstructed. This is very important since for a given structure, the state of
* Corresponding author. Tel.: +965-481-1188, ext. 5790; fax: +965-481-7131. E-mail address: [email protected] (A.P. Christoforou).
damage depends on the impact characteristics. Furthermore, the identi®cation, estimation, and reconstruction of the response are essential for controlling the severity of impact and health monitoring of the structure. For these reasons, the identi®cation of the impact force has received much attention in recent years. In early papers, Doyle [1,2] used a frequency-domain deconvolution technique to estimate the impact force on beams and plates from strain gauge measurements. The fast Fourier transform (FFT) was used to convert the measured strain history into its frequency components, which were used to obtain the frequency components of the impact force. Then the inverse fast Fourier transform (IFFT) was implemented to reconstruct the impact force history in the time-domain. Chang and Sun [3], on the other hand, used experimentally generated Green's functions and time-domain deconvolutions to estimate the impact force history from strain measurements. A least-square ®t method was used in their numerical scheme. However, the amplitudes of the Green's functions had to be scaled to account for dierent impactor characteristics. Hollandsworth and Busby [4] used acceleration measurements to identify the impact force on aluminum beams by solving a dynamic inverse problem. Wu et al. [5] presented two methods for the purpose of detecting the impact force history. In the ®rst method,
0045-7949/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 1 ) 0 0 1 1 3 - 4
2608
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
the classical plate theory, and the normal-mode superposition method using Bessel functions were used to construct Green's functions in the time-domain. A constrained optimization technique was used to reconstruct the impact force history from measured strain histories. This method was limited to axisymmetric impact problems because of the restriction of the Green's functions. In the second method, the impact force history was reconstructed using the gradient projection method by comparing the response to recorded force and strain histories of a controlled impact. Tsai et al. [6] also used a Green's function approach to estimate impact force from experimental data on sandwich panels. With the advent of smart structure technology, Choi and Chang [7] presented an impact load identi®cation method to estimate the location and history of an impact force using piezoelectric sensor measurements. Their identi®cation process consisted of a system model and a response comparator, employing a quasi-Newton procedure to estimate the location of impact, and the leastsquare ®t method to predict the force history. Lihua and Baoqi [8] used neural networks and piezoelectric sensors to identify the location and amplitude of a static force on a square plate. The scope of this work is to provide an easy and reliable methodology to estimate the mass and velocity of a foreign object impacting a structure. The inverse model presented below, employs the direct model developed by Christoforou and Yigit [9,10] and a least square optimization technique to reconstruct the impact force history from experimental data. It is worth noting that in this study, the impactor was taken to be spherical in shape with ®xed diameter but not a sphere as shown in Fig. 1. The problems that are considered here are low velocity impact events, e.g., tools dropped on ¯exible composite structures. Usually the range of impactor diameters is between 10 and 25 mm. In this range, which is the case in the present study, the contact is not visible and has no major in¯uence on the impact response.
Previous studies showed that the relative stiness of the structure (which is a measure of ¯exibility), and the mass ratio of the impactor to the structure signi®cantly aect the impact response [9,10]. The in¯uence of the diameter of the impactor will be in the contact zone, which is taken care by the contact law (especially in impacts of rigid bodies or very sti structures). For smaller diameters the impact will be visible in the form of a damage in the contact zone, and identi®cation is not necessary. Therefore, the mass of the impactor rather than its diameter was considered to be one of the parameters to be estimated by the proposed inverse analysis. Furthermore, the present study is limited to the case of small displacements and with no damage due to impact. However, the proposed scheme is simple, fast and accurate in estimating the impactor characteristics as demonstrated in the experimental results used to validate the model. Certainly, and this is beyond the scope of this initial work, the direct model may be improved by including not only geometrical non-linear eects but also eects due to structural damage, which is more important in impact of composite structures. Also, the scheme may be applied not only to force measurements, but also to other measurements such as from strain gauges or piezoelectric sensors. 2. Theory 2.1. Direct impact model The equations of motion of a rectangular, specially orthotropic, symmetrically laminated, elastic and shear deformable composite plate, subject to lateral dynamic loading, are given as D11
o2 wy o2 wx o2 wx D D D 66 12 66 ox2 oy 2 ox oy ow qh3 o2 wx kA55 wx ox 12 ot2
Fig. 1. Geometry of the impact problem.
1
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
o2 wy o2 wy o2 wx D66 2 D22 2 D12 D66 ox oy ox oy ow qh3 o2 wy kA44 wy 12 ot2 oy kA55
owx o2 w 2 ox ox
qh
kA44
owy o2 w 2 oy oy
o2 w ot2
1 1 m2i 1 E Ez Ei 2
p x; y; t 3
where Dij and Aij are the bending and shear stinesses de®ned as usual, h is the plate thickness, q is the material density, w is the transverse de¯ection, wx and wy are the shear rotations, p is the lateral load per unit area, k is the Mindlin shear correction factor, x and y are the space variables, and t is the time. Assuming a concentrated impact loading, the lateral load per unit area can be written as p x; y; t F t d x
xc d y
yc
4
where F t is the impact force, d is the Dirac-delta function, and xc ; yc is the contact point. In the case of lateral impact of the plate by a foreign object, the motion of the impactor is described by m wi
F t
5
where m and wi are the mass and the displacement of the impactor, respectively. The impact force can be expressed as follows: F t Fc
if Fc P 0
F t 0 if Fc < 0
6
where Fc is to be obtained from an appropriate contact model. In this study a linear version of the elasto-plastic contact model proposed by Yigit and Christoforou [11,12] is used, and is given as Fc a Ky a
7
where a is the local indentation which is the relative approach between the impactor and the target of the contact point given as a t wi t
w xc ; yc ; t
Ky is the contact stiness given as p Ky 1:5Kh ay
2609
8
9
where Kh is the Hertzian stiness for the case of spherical body of radius R, in contact with a ¯at surface, given as p Kh 43 RE 10 where R is the radius of the spherical impactor and E is given as
11
where mi and Ei are the Poisson's ratio and Young's modulus for the impactor, respectively, Ez is the transverse Young's modulus of the laminate, which is usually taken to be equal to that of the ply, i.e., Ez E22 . The onset of ``yielding'' or damage was obtained by combining the classical Hertzian analysis and the elastic± plastic contact theory. By assuming the material is elastic before damage and perfectly plastic after damage occurs, the critical indentation ay was obtained using the maximum shear failure criterion and is given as [11] ay
0:68Sy2 p2 R E2
12
where Sy is the yield strength of the softer material for metals, and Sy 2Su , where Su is the shear strength of the laminate for the ®ber composites. This paper is concerned with a simply supported rectangular composite plate of uniform thickness with dimensions a and b, impacted by a spherical impactor as shown in Fig. 1. The boundary conditions in this case can be expressed as follows: w
owx 0 at x 0; a ox
13a
w
owy 0 at y 0; b ox
13b
The initial conditions of the impact problem are: w x; y; 0 0; w_ i 0 v0
w_ x; y; 0 0;
wi 0 0;
and 14
where v0 is the initial velocity of the impactor. The equations of motion are discretized by a standard modal analysis procedure, and numerically solved to obtain the response. 2.2. Inverse impact model The direct impact analysis shows that the impact response depends on the properties of the structure and the impactor. Furthermore, for a given structure and impact energy, the mass of the impactor is the parameter that governs the impact response (i.e., from locally dominated to globally dominated response, [11,12]). In direct simulations of the response (either numerical or experimental), all impact parameters are known. In the inverse analysis, on the other hand, the response is known usually from experimental data, with some of the impact parameters unknown. In this study, the unknown parameters to be estimated are the mass and the velocity of the impactor, and the predicted values will be denoted m and v0 , respectively. All other plate and
2610
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
Fig. 2. Flowchart of the proposed inverse impact algorithm.
impactor parameters are known. The inverse problem is solved by an optimization technique, where the experimental impact response is compared to the response predicted by the direct solution for an m; v0 pair. The optimal solution is achieved by minimizing the following objective function: f m; v0
N X F m; v0 ; ti i1
subject to
Fe ti 2
15
0 < m 6 mL
16
0 < v0 6 vL
17
where N is the total number of readings, F m; v0 ; ti and Fe ti are the predicted and experimental force±time histories, respectively, ti is the time at ith reading, and mL and vL are limiting values for the mass and velocity of the impactor, respectively. In this study, a M A T L A B routine (version 5.2) called ``leastsq'' which performs least squares optimization was
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
used to obtain the optimal solution. In this routine, both the Gauss±Newton method and the Levenberg±Marquart method are included as options. Basically the dierence between these two methods is the way of the search direction for the optimal solution. For the results presented in this paper Gauss±Newton method was used. Starting from an initial guess for an m; v0 pair, the optimizer will progress towards the solution by solving the direct impact problem to obtain F m; v0 ; ti . The design variables m and v0 are varied using the Gauss±Newton optimization algorithm. This process is guided by objective function gradient evaluations in the design space, and will continue until the objective function reaches a minimum value based on the speci®ed convergence criterion. Fig. 2 shows the ¯owchart for the proposed inverse algorithm. For practical reasons, and without loss of generality, it should be noted here, that the range of the mass and velocity of the impactor is restricted to the range of low velocity impact. A simple adjustment is implemented to force the optimizer to reverse the search direction when any one of the design variables crosses the boundary of the allowable space. This is achieved by deliberately assigning a value to the objective function higher than its previous value.
3. Results and discussion As mentioned earlier, the design parameters of the identi®cation problem are the mass m and the velocity v0 of the impactor. Although the eect of these two parameters on the impact response is well known, at this point of presentation, it is important to show these effects, with the following set of simulations. The material properties, geometry of the composite plate, and the impact conditions are given in Table 1. Fig. 3 shows the eect of the mass of the impactor on the impact force-time response. The velocity v0 of the
Table 1 Properties of the composite plate and the impact parameters Plate 0=90=0=90=0s , T300/934 carbon-epoxy Plate size: a 200 mm; b 200 mm Plate thickness: h 2:69 mm (0.269 mm/layer) E11 120 GPa; E22 7:9 GPa; G12 G23 5:5 GPa m12 m23 0:3; q 1580 kg/m3 Impactor 12.7 mm diameter spherical; q 7960 kg/m3 , steel Contact parameters Kh 8:89 108 N/m1:5 ; ay 2:5 10 Ky 6:65 106 N/m
5
m
2611
impactor was held ®xed at 3 m/s. For small mass m 0:005 kg the duration of the impact is very short (see Fig. 3a). As the mass of the impactor increases, multiple impacts occur, and the duration and the amplitude of the impact force increase. Fig. 4 shows the eect of the velocity v0 of the impactor on the force-time response. The mass of the impactor was held ®xed at 0.5 kg. The results show that the impactor velocity has signi®cant in¯uence on the impact force amplitude. However, it has no in¯uence on the duration of the impact. From the results presented in Figs. 3 and 4, it is clear that the impact response is very sensitive to the mass and velocity of the impactor. Therefore, it will be easy to predict these variables accurately using the inverse technique. In order to validate the inverse model proposed in this paper, two experimental results reported by Pierson and Vaziri [13], were used as inputs to predict the mass and velocity of the impactor. In the experiments, simply supported 45=90= 45=03 s T800/TH/3900-2 CFRP plates with nominal dimensions of 127 76:2 4:65 mm3 were used. In the ®rst experiment, the exact mass and the velocity of the impactor were 6.15 kg, and 1.76 m/s, respectively. Fig. 5 shows a comparison between the impact force responses from the experimental data, the direct solution based on the exact impactor mass and velocity values, and the inverse model prediction. As can be seen there is a small discrepancy between the direct solution and the experiment. As expected, this discrepancy in¯uences the inverse model prediction, and the estimated mass and velocity of the impactor were 5.41 kg and 1.81 m/s, respectively. The eect of the initial guess on the inverse solution was also investigated. It was found that the proposed inverse algorithm is independent of the initial guess. As it will be discussed later, a good initial guess, however, could speed up convergence. Fig. 6 displays the convergence history of the inverse model for two initial guesses, (m 1 kg and v0 5 m=s) and (m 0:5 kg and v0 3 m=s). As can be seen only 30 iterations were needed for convergence to the optimal solution. In the second experiment, the mass and the velocity of the impactor were 0.314 kg, and 7.7 m/s, respectively. Fig. 7 shows a comparison between the impact force responses from the experimental data, the direct solution based on the exact values of the impactor mass and velocity, and the inverse model prediction. The estimated mass and velocity of the impactor were 0.28 kg and 8.098 m/s, respectively. Fig. 8 displays the convergence history of the inverse model for two initial guesses, (m 0:5 kg and v0 3 m=s) and (m 1 kg and v0 5 m=s). These results con®rm further the conclusions drawn from Fig. 6, i.e., the inverse solution is unique and independent of the initial guess. There was no local minimum solution observed in all cases which is in agreement with the physics of the impact problem.
2612
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
Fig. 3. Eect of the mass m of the impactor on the impact force±time response F t.
Again only 30 iterations were needed for convergence to the optimal solution. At this point, it is useful to explore ways for providing guidelines for better initial guesses and/or computationally faster direct solutions in the inverse algorithm. To this end, the characterization procedure introduced by Christoforou and Yigit [9,10] can be used.
In their procedure, it was found that the non-dimensional maximum impact force F max is governed by two non-dimensional parameters, the relative stiness k, and the loss factor f. These parameters are given as F max
v0
F pmax mKy
18
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
2613
Fig. 4. Eect of the velocity v0 of the impactor on the impact force±time response F t.
k
f
Kst Ky 1 16
where
mKy qhD
p D11 D22
19
D 12 D12 2D66
20
and Kst is the global static stiness of the plate. Applying the procedure, the relative stiness of the plate used in the experiments was calculated to be k 0:258. The ®rst experiment is considered as a heavy mass impact (m 6.15 kg), and f and F max were calculated to be 10.61 and 0.45, respectively. Using the estimated
1=2
21
2614
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
Fig. 5. Comparison between the predicted and the actual impact force histories for the ®rst set of experimental data [13].
mass and velocity values, f and F max were calculated to be 9.95, and 0.47, respectively. The second experiment is considered as a medium mass impact (m 0.314 kg) and f and F max were calculated to be 2.4 and 0.513, respectively. The predicted f and F max were 2.26 and 0.515, respectively. Both cases are displayed on a characterization curve as shown in Fig. 9. As can be seen, the agreement between the experiments and inverse prediction is excellent. Also, the ®rst case is located in the so-called fully global region (quasi-static), where the response can be predicted by using a simple one-degree of freedom system. The second case is located in the socalled transition region in which the response is dynamic in nature. Fig. 9 shows the functional relationship between F max , f and k. In the inverse solution f and F max are not
Fig. 6. Convergence history of the inverse model for two initial guesses, for the ®rst set of experimental data [13].
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
2615
known because they depend on the mass and velocity of the impactor. The only parameter that is known is the relative stiness k. Examination of the ®gure reveals that, for practical purposes and especially for sti structures, many impact situations may be considered fully global for a range of values of f. In fully global situations, the response is a function of k only, and is given in a normalized form as
F s
Fig. 7. Comparison between the predicted and the actual impact force histories for the second set of experimental data [13].
r! r ! k k sin s k1 k1
22
In this case F max and the non-dimensional response frequency xr can be written as
Fig. 8. Convergence history of the inverse model for two initial guesses, for the second set of experimental data [13].
2616
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
Fig. 9. Characterization curve for the impact response. Fig. 10. Comparison between the predicted and the actual impact force histories for the third set of experimental data [2].
r k k1
23
r p=tr k xr p k 1 Ky =m
24
F max
F 2fxn F_ x2n F 0;
where tr is the contact duration. Therefore, only Fmax and tr , which can be determined from the experimental response, are needed to estimate the mass and velocity of the impactor. Applying Eqs. (23) and (24) to the ®rst experimental set, the mass and velocity of the impactor were estimated to be 5.45 kg and 1.86 m/s, respectively, which compare well with the full inverse model prediction. Although the second experimental set is dynamic in nature, the same procedure can be applied to obtain a good initial guess. In this case the mass and velocity of the impactor were estimated to be 0.324 kg and 8.59 m/s, respectively. In many cases, this kind of estimations may be adequate. If not, they may serve as initial guesses to the full inverse model. As a rule of thumb, in most impacts having time duration in the milliseconds region, the quasi-static approach can be used with a reasonable accuracy. In the transition region, however, the quasistatic model can provide a good initial estimate. Further study of the ®gure, in the in®nite plate and the local regions, the impact response is a function of f only. In the local region, there are only single short duration impacts (microseconds) and the non-dimensional impact force can be expressed as F s sin s
25
This means that F max 1, and xr 1. As before, the mass and velocity of the impactor can be estimated easily. In the case of impact of an in®nite plate, the response is governed by
xn 1
26
The solution of Eq. (26) will be the direct solution in the inverse algorithm. To validate this model, an experimental set reported by Doyle [2] was used. In the experiment, a seven layer [0/90] Schotchply glass/epoxy composite plate with nominal dimensions of 406 406 1:78 mm3 was used. The mass and velocity of the impactor were 8.537 g and 2.3 m/s, respectively. Fig. 10 shows a comparison between the impact force responses from the experimental data, the direct solution based on the exact values of the impactor mass and velocity, and the inverse model prediction using Eq. (26). The estimated mass and velocity of the impactor were 9.9 g and 2.1 m/s, respectively. The values of f and F max were calculated to be 2.085 and 0.22, respectively for the experimental data, and 2.25 and 0.197, respectively for the inverse model prediction. This set is also displayed on the characterization curves of Fig. 9. As can be seen, the agreement is excellent for both the values and the type of response.
4. Conclusions An inverse model for the low-velocity-impact of composite plates has been presented. An ecient direct solution was utilized in a least squares optimization technique to predict the mass and the velocity of the impactor. A M A T L A B routine called ``Leastsq'' which performs least squares optimization was used to obtain the optimal solution. Three sets of experimental data were used to validate the accuracy of the proposed inverse algorithm. The results showed that the convergence was fast and independent of the initial guess. Also, using a characterization procedure indicating the type of
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
response, excellent guidelines can be provided for better initial guesses, and for simpler and computationally faster direct solutions in the proposed inverse algorithm.
2617
Load experimental data ®le consisting of 2 columns.The ®rst column corresponds to the experimental time data points and the second column represents the corresponding measured impact force load piers1.dat
Appendix A. tion
MATLAB
computer code for inverse solu-
This program solves the inverse impact problem to predict the the impactor mass and velocity from experimentally measured impact force time history. Variables de®nitions: time_exp experimental time data points force_exp experimental force measurements ndata_exp number of experimentally measured data points t®nal_exp time of last experimentally measured data point fact1, fact2 parameters to control the least square search steps error_last last computed error value. Used to force the optimizer to reverse the search direction when the boundaries of the allowable space are reached itr iteration counter ntimes number of times the size of the time vector is enlarged in order to accommodate the largest expected impact duration that may be considered by the optimizer conv_history a matrix to record the convergence history vel0, mass0 initial guesses for the velocity and mass, respectively initail_guess scaled initial guess vector options vector of parameters to control the optimization routines (see Matlab help). sol predicted optimum solution vector (scaled) v0star, mstar predicted optimum mass and velocity, respectively vel_mass Optimum solution vector time, force time and force vectors resulting from the direct solution of the impact problem using the predicted velocity and mass Note: The ®les `impact.exe' and `impact_inv.m' have to be available in the default directory. Clear the workspace and close all ®gures clear; close all; Declaration of global variables global time_exp force_exp t®nal_exp ndata_exp fact1 fact2 error_last global itr ntimes conv_history
ntimes 4; Assign the experimental data points to their corresponding vectors of time and force time_exp piers 1(:,1); force_exp piers 1(:,2); Plot the experimental impact response ®gure(1); plot(time_exp, force_exp,`+'), hold on determine ndata_exp and t®nal_exp ndata_exp length (time_exp); t®nal_exp time_exp (ndata_exp); Initialize the convergence history matrix conv_history [ ] Assign appropriate values to fact1 and fact2 fact1 0.000001; fact2 0.00001; Assign and scale the initial guesses for the velocity and mass vel0 3; mass0 .5; initial_guess [vel0*fact1,mass0*fact2]; Assign the optimization parameters options(2) 1e-7; options(3) 1e-7; options(5) 1; options(14) 200; Initialize the iteration counter itr 0; Solve the optimization problem using the least squares `leastsq' Matlab routine.`impact_inv' is a Matlab function (see listing below) that evaluates the objective function for any given set of the parameters to be predicted, i.e.the velocity and mass.It is used by the least squares optimization routine. [sol, options] leastsq (`impact_inv', initial_guess, options); Compute the predicted optimum velocity and mass v0star sol(1)/fact1; mstar sol(2)/fact2; fprintf(`Predicted velocity and mass: %g %g n', v0star, mstar) vel_mass [v0star mstar];
2618
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619
Save the predicted parameters into a data ®le named `param.dat'.This ®le will be read by the executable Fortran code `impact.exe' designed to solve the direct impact problem save param.dat vel_mass -ascii Solve impact problem for the predicted velocity and mass using the executable Fortran code `impact.exe' in order to reconstruct the corresponding force time history.Results will be available in the ®le `Result.dat' !impact Load results load result.dat Assign the results to their corresponding time and force vetors time result (:,1); force result (:,2); Plot the predicted impact response plot(time,force,`g'), hold on Save the convergence history into a data ®le named `history.dat' save history.dat conv_history -ascii function err impact_inv(params): `impact_inv' is a Matlab function that evaluates the objective function (norm of the force error relative to the experimentally measured impact force response) for any given set of the parameters to be predicted, i.e.the velocity and mass.It is used by the least squares optimization routine. Declaration of global variables global time_exp force_exp t®nal_exp ndata_exp fact1 fact2 error_last global itr ntimes conv_history Calculate the velocity and mass velocity params (1)/fact1; mass params (2)/fact2; Check for out of range parameters and take appropriate action if velocity<1 velocity>10 mass<0.00001 mass>10 err error_last*1.001;return end Solve the direct impact problem vel_mass [velocity mass]; save param.dat vel_mass -ascii !impact load result.dat
time result(:,1); force result(:,2); n length(time); t®nal time(n); Interpolate the predicted force response at the experimental time data points.Then adjust the sizes of the experimental and predicted force and time vectors so that they have the same size. if t®nal
References [1] Doyle JF. Further developments in determining the dynamic contact law. Exp Mech 1984;24:265±70. [2] Doyle JF. Experimentally determining the contact Force during the transverse impact of an orthotropic plate. J Sound Vibrat 1987;118(3):441±8. [3] Chang C, Sun CT. Determining transverse impact force on a composite laminate by signal deconvolution. Exp Mech 1989;29:414±9. [4] Hollandsworth PE, Busby HR. Impact force identi®cation using the general inverse technique. Int J Impact Engng 1989;8(4):315±22. [5] Wu E, Tsai TD, Yen CS. Two methods for determining impact-force history on elastic plates. Exp Mech 1995;35: 11±8. [6] Tsai CZ, Wu E, Luo BH. Forward and inverse analysis for impact on sandwich panels. AIAA J 1998;36(11): 2130±6.
A.P. Christoforou et al. / Computers and Structures 79 (2001) 2607±2619 [7] Choi K, Chang FK. Identi®cation of impact force and location using distributed sensors. AIAA J 1996;34(1):136± 42. [8] Lihua S, Boaqi T. Identi®cation of external force acting on a plate by using neural networks. In: Chang F-K, editor. Proceedings of the International Workshop on Structural Health Monitoring, Stanford University, Stanford, CA, 18±28 September, 1997. p. 229±37. [9] Christoforou AP, Yigit AS. Characterization of impact in composite plates. Compos Struct 1998;43:15±24.
2619
[10] Christoforou AP, Yigit AS. Eect of ¯exibility on low velocity impact response. J Sound Vibrat 1998;217(3):563±78. [11] Yigit AS, Christoforou AP. On the impact of a spherical indenter and an elastic±plastic transversely isotropic halfspace. Compos Engng 1994;4:1143±52. [12] Yigit AS, Christoforou AP. Impact dynamics of composite beams. Compos Struct 1995;32:187±95. [13] Pierson MO, Vaziri R. Analytical solution for low-velocity impact response of composite plates. AIAA J 1996;34: 1633±40.