# Vibration Response Evaluation under Shock-Type Loading with Emphasis on Finite Element Model Updating

Authors: Jalal Akbari; Leila Nazari; Samaneh Mirzaei
Journal: Shock and Vibration
 (2020)

---

## Abstract

In some cases, impulse- or shock-type excitations as the dynamic loading are inevitable, and obtaining proper response with the well-known numerical methods is not easy. This paper focuses on dynamic response estimation against short-time loading with an updated finite element model using frequency response functions (FRF) and particle swarm optimization (PSO) technique. Because there is not an analytical method for assessing the numerical responses under shock-type excitations, in this paper, experimental tests are designed on a laboratory scale to evaluate the numerical responses. The vibration responses of the system against shock loading are compared with the Newmark average acceleration scheme and also with experimental data. The results reveal that the unconditionally stable Newmark method against regular loads has an appropriate performance. Still, under short-time loading, it faces numerical damping error, and this method should not be blindly applied under shock-type loads.


---

## Body

## 1. Introduction

Direct numerical integration technique is commonly used for response estimation of any structural systems. The effective parameters of any dynamic problems are mass, damping, and stiffness matrixes and also applied load in the single-degree-of-freedom (SDF) or multi-degrees-of-freedom (MDF) systems. Usually, these matrixes are obtained from finite element modeling (FEM) of any structural systems, and modeling errors include idealization in FEM, size of the finite element (FE) mesh, material properties, and inaccuracies in boundary conditions. To reduce the mentioned errors and obtain positive responses from numerical modeling, updating the mass and stiffness matrixes is a crucial process because those are inherent properties of the system and they affect structural responses and behavior. The modal parameters, e.g., natural frequencies, mode shapes, and also modal damping, are required for updating the numerical models, and obtaining these parameters is called system identification. Various useful approaches have been presented to extract the modal parameters of a structure in the time or frequency domain [1]. In addition to the numerical errors, the measurement errors occur because of the environmental effects, sensor errors, random errors, and errors in signal processing. As the modeling errors and measurement errors can lead to poor structural parameter identification, it is required to validate the model updating experimentally. Model updating methods simultaneously utilize the structural response obtained by the finite element method (FEM) and the measured structural response to calibrate mathematical modeling. However, in this research, it is assumed that the experimental data are clean due to the appropriate controls during the experimental tests.Model updating could be classified into modal-based and response-based methods. The modal-based model updating technique relies on the modal characteristics data obtained from an experimental modal analysis that is indirectly extracted from the measured FRF data. D’Ambrogio and Sestieri [2] emphasize that the numerical procedures used for modal identification using experimental vibration data can introduce errors exceeding the level of required accuracy to update FE models. In response-based finite element model (FEM) updating methods, the measured FRF data are directly utilized to identify the unknown structural parameters. Most of the FRF-based model updating techniques are used to minimize a residual error between the analytical and experimental input force and output response [3–7]. Many researchers utilized vibration data of structures for finite element (FE) model updating. Baghchi [8] used vibration data for updating the mathematical model. Lin and Zho [9] updated the FE model using vibration data under base excitation. Pradhan and Modak [10] modified the RFM and presented a normal FRF concept to update mass and stiffness matrices of structure. For the damped system model updating, several investigations are conducted. Lin and Zho [11], Arora et al. [12], and Yuan and Yu [13] applied the FRF approach to update the FE model in a damped system based on vibration information. Garcia and Santini [14] updated damped structural systems with a two-step model-updating algorithm.Many kinds of researches are focused on model updating from a mathematical point of view. Sipple and Sanayei [15] used numerical sensitivity techniques for model updating using FRF. They used mathematical-based sensitivity method for gradient calculations. Wei-Ming and Jia-Zhen [16] updated the FE model based on model reduction and iterative approach. Weng et al. [17] applied the inverse substructure method for model updating. In their research studies, a system was divided into several substructures, and updating is applied for each subsystem. Papadimitrious and Papdioti [18] used the component mode synthesis technique for the FE model updating. Wan and Ren [19] utilized a residual-based Gaussian process for model updating. Sarmadi et al. [20] implemented an iterative least-squares minimal residual technique for model updating, and Wei et al. [21] updated multi-degree-of-freedom systems based on the intrinsic chirp component decomposition method.Also, model-updating techniques are usually involved with optimization methods or algorithms. In this regard, there are much related researches. Ntotsios and Papadimitrio [22] applied the multiobjective optimization technique to FE model updating based on eigenvalue strain energy residuals. Christodoulou et al. [23] used the Pareto optimization technique for structural model updating; Jung and Kim [24] used a hybrid genetic algorithm for finite element model updating. Shabbir and Omenzetter [25] utilized particle swarm optimization with a sequential niche technique for dynamic finite element model updating.Various researchers studied vibration responses of the systems under impact loading. Tagarielli et al. [26] predicted the dynamic response of composite sandwich beams under shock loading. Wali et al. [27] obtained the vibration response of the sandwich plate under low-velocity impact loading. Wang et al. [28] estimated the dynamic responses of clamped sandwich beams subjected to impact loading with analytical, numerical, and also experimental studies. Mazurkiewicz et al. [29] compared the numerical responses with experimental testing data of individual structural elements under impulse loading.To the knowledge of the authors, numerical studies were not conducted to calculate the vibration responses of MDF systems with the commonly used Newmark method against impact or impulse loads [30]. Therefore, an experimental test is designed on a laboratory scale to evaluate and validate the numerical responses of the Euler–Bernoulli-type beam with experimental data. Thus, vibration responses of a cantilever aluminum beam were studied with analytical, updated numerical and also experimental methods. For this purpose, at first, the analytical formulation of the free vibration of the beam was derived. Then the numerical model with the finite element technique was created, and modal data were extracted from the FRF response. The beam was excited with a hammer, and the responses were measured by using the acceleration sensors at four reference points. The updating method was based on the particle swarm optimization (PSO) technique. After updating the finite element model, the dynamic responses of the beam against impact loading were evaluated with the Newmark method and experimental data.


## 2. Problem Definitions

### 2.1. Analytical Free Vibration Formulation

In this section, the analytical free vibration formulations for modal data extraction are presented based on the differential equation of vibration of the Euler–Bernoulli-type beams. For this purpose, according to Figure1, a cantilever beam with length L, a rectangular cross section with width b and height h, a concentrated mass M at its free end, and stiffness equal to EI is considered.Figure 1
Cantilever beam with a mass concentrated at the free end.In the analytical method, by solving the equilibrium equation of the beam for free vibration, the frequencies and the mode shapes are determined. The dynamic equation of the Euler–Bernoulli beams with a uniform cross section in the forced-vibration case is obtained as follows:(1)EI∂4wx,t∂x4+ρA∂2wx,t∂t2=px,t,where ρA is equal to the mass of unit length or m¯ and p(x, t) is the applied force to the arbitrary point. The parameter w(x, t) also refers to the beam displacement response equation along the z direction at any point in the time domain. For free vibration analysis, the right side of equation (1) should be equal to zero. As a result, equation (1) is rewritten as follows:(2)a2∂4wx,t∂x4+∂2wx,t∂t2=0;a=EIρA.To solve equation (2), the method of separation variable is used as follows:(3)wx,t=ψxzt,a2ψIVxzt+ψtz¨t=0.This equation is transformed into two separate equations in the time domain (t) and the space domain (x) as follows:(4)z¨t+ω2zt=0,ψIVx−β4ψx=0.The new variables are defined as follows:(5)−a2ψIVxψx=z¨tzt=−ω2,β4=ω2a2.By solving equation (4), the responses are obtained according to the following equation:(6)zt=Acosωt+Bsinωtψx=G1cosβx+G2sinβx+G3coshβx+G4sinhβx⟶wx,t=ψxzt,where zt is the time-domain response and ψx is the eigenvector of the beam. The boundary and initial conditions of the beam are required to determine the coefficients Gi and A, B. According to Figure 1 and equation (7), in the left boundary of the beam, the displacement and slope are zero:(7)w0,t=0∂w∂x0,t=0⟶G1+G3=0G2+G4=0.Based on equation (7), the mode shape is rewritten as follows:(8)ψx=G1cosβx−coshβx+G2sinβx−sinhβx.The boundary conditions of the cantilever beam with a concentrated massM at the free end are as follows:(9)Ml,t=0∂∂xEI∂2wl,t∂x2=M∂2wl,t∂t2⟶EI∂2wl,t∂x2=0EI∂3wl,t∂x3−M∂2wl,t∂t2=0.The resultEIw″l,t=0 and EIw‴l,t=Mw¨l,t could be rewritten as equations (10) and (11), respectively:(10)G1cosβl−coshβl+G2sinβl−sinhβl=0,(11)EIMψ‴l=ψlz¨tzt⟶−EIMω2ψ‴l=ψl.By replacing the valuesψl,ψ″l into equation (11), the following equation is obtained:(12)EIβ3Mω2−G1sinβl+sinhβl+G2cosβl+coshβl=G1cosβl−coshβl+G2sinβl−sinhβl.By replacingλ=βl, α=m¯l/Mλ, and β4/ω2=EI/m¯ in equation (10), equation (12) is expressed as follows:(13)cosλ+coshλsinλ+sinhλα−sinλ+sinhλ−cosλ+coshλαcosλ+coshλ−sinλ+sinhλG1G2=00.To obtain the eigenvalues or frequencies (λ), the determinant of equation (13) should be zero as follows:(14)αcosλ+coshλ2+cosλ+coshλ−sinλ+sinhλ+αsinh2λ−sin2λ+sinλ+sinhλ−cosλ+coshλ=0.The simplified form of equation (14) is expressed as equation (15). By solving this equation, the eigenvalues (λn) are available. This nonlinear equation requires special numerical techniques to obtain all roots of eigenvalues [31]:(15)sinλncoshλn−cosλnsinhλn1+cosλncoshλn−m¯lMλn=0.According to the following equation, the frequency (ωn) in each mode is available:(16)λn2=ωnl2m¯EI⟶ωn=λnl2EIm¯.To obtain the eigenvectors of the beam, one should calculate the coefficientsG1andG2. In this study, G1 is set to the unit, and the general unknown factor Gn is computed from the following equation:(17)Gn=αncosλn+αncoshλn−sinλn+sinhλnαnsinλn−αnsinhλn+cosλn−coshλn,λn=βnl:αn=m¯lMλn=m¯βnM.Therefore, the eigenvectors or equivalently mode shapes of the studied beam in Figure1 is available according to the following equation:(18)ψnx=sinβnx−sinhβnx+Gncosβnx−coshβnx,ψnx=sinλnxl−sinhλnxl+Gncosλnxl−coshλnxl.


### 2.2. Numerical Modal Parameters

The finite element method was used to obtain the modal parameters in the numerical method. Firstly, using the finite element method, the studied beam in Figure1 was divided into twenty-five equal parts with 3 centimeters length. Each element has two nodes, and each node has a transitional degree of freedom (wi) and a rotation degree of freedom (θi). Stiffness matrix Ke and mass matrix Me of each element are determined as follows:(19)Ke=EIl31261−126161412−61212−12−6112−6161212−61412;Me=ρAl42015622154−131221412131−31254131156−221−131−312−221412,where E is the modulus of elasticity, I is the moment of inertia of the cross section, and l is the length of each element (l=L/Nelem), where L is the length of the beam, Nelem is the number of elements, ρ is the mass density of the material, and A refers to the cross section of the beam. It should be noted that the last element of the beam has a concentrated mass (M) at the second node. Therefore, during the assembling process of the total mass matrix, this concentrated mass should be added to the mass matrix as follows:(20)Mcon=0000000000M00000.The beam’s length was equal to 75 cm, and the number of elements was 25. To compare the results and update the finite element model, in the second stage, the number of elements reduced to four elements.


### 2.3. Experimental Modal Parameters

The aluminum beam, as mentioned in Section2.1, was also studied in an experimental study. As shown in Figure 2, the acceleration data were collected by installing unidirectional acceleration sensors at four points of the beam and applying the hammer force at the specified positions.Figure 2
Overall geometry of the beam and the locations of acceleration sensors (a) and hammer (b) in the experimental study.
(a)(b)According to Figure3, the input data (forces) and output data (accelerations) are recorded at positions P1, P2, and P3, and they are collected at 4 stations (CH1, CH2, CH3, and CH4). The entire data collection process was 30 seconds. Each test is repeated three times. Based on the data sampling frequency of the data acquisition system (400 Hz), for each test, 12,000 points are recorded.Figure 3
Location of exerted forces (P1, P2, and P3) and acceleration sensors (CH1, CH2, CH3, and CH4) in the studied beam.Extraction of modal parameters was accomplished using frequency response function (FRF). For example, FRFs; in the case that the hammer exerted at the free end (Input 1), the acceleration data (Outputs 1, 2, 3, and 4) are presented in Figure 4. The required experimental modal data are extracted from FRFs and written as codes in MATLAB [32].Figure 4
The frequency response functions (FRFs) for experimental tests in the studied beam. (a) Output 1, Input 1. (b) Output 2, Input 1. (c) Output 3, Input 1. (d) Output 4, Input 1.
(a)(b)(c)(d)(e)(f)(g)


### 2.4. Updating the Finite Element Model

The aim of updating the finite element model is correcting or improving the finite element model of a structure in terms of updating the global or local mass and stiffness matrices, boundary conditions, and geometrical and material properties. For this purpose, the updating process was accomplished with the modal assurance criterion (MAC) method [5] based on the recorded data at four stations in the experiment program to correcting the global mass and stiffness matrices. As well, the finite element model is reduced to coincide with the experimental setup with the SEREP algorithm [33]. Update has been conducted for mass and stiffness matrices using the particle swarm optimization (PSO) technique [34]. The objective function is defined based on the frequencies and the components of the mode shapes in the experimental and numerical data as follows:(21)f=ωinum−ωiexpωiexp2+∑j=1ndfϕi,jnum−ϕi,jexpϕi,jexp2,where ωinumandϕi,jnum are the frequencies and the jth components of mode shapes of the finite element model, respectively, and ωiexpandϕi,jexp are the frequencies and the jth mode shapes of the experimental data, respectively.


## 2.1. Analytical Free Vibration Formulation

In this section, the analytical free vibration formulations for modal data extraction are presented based on the differential equation of vibration of the Euler–Bernoulli-type beams. For this purpose, according to Figure1, a cantilever beam with length L, a rectangular cross section with width b and height h, a concentrated mass M at its free end, and stiffness equal to EI is considered.Figure 1
Cantilever beam with a mass concentrated at the free end.In the analytical method, by solving the equilibrium equation of the beam for free vibration, the frequencies and the mode shapes are determined. The dynamic equation of the Euler–Bernoulli beams with a uniform cross section in the forced-vibration case is obtained as follows:(1)EI∂4wx,t∂x4+ρA∂2wx,t∂t2=px,t,where ρA is equal to the mass of unit length or m¯ and p(x, t) is the applied force to the arbitrary point. The parameter w(x, t) also refers to the beam displacement response equation along the z direction at any point in the time domain. For free vibration analysis, the right side of equation (1) should be equal to zero. As a result, equation (1) is rewritten as follows:(2)a2∂4wx,t∂x4+∂2wx,t∂t2=0;a=EIρA.To solve equation (2), the method of separation variable is used as follows:(3)wx,t=ψxzt,a2ψIVxzt+ψtz¨t=0.This equation is transformed into two separate equations in the time domain (t) and the space domain (x) as follows:(4)z¨t+ω2zt=0,ψIVx−β4ψx=0.The new variables are defined as follows:(5)−a2ψIVxψx=z¨tzt=−ω2,β4=ω2a2.By solving equation (4), the responses are obtained according to the following equation:(6)zt=Acosωt+Bsinωtψx=G1cosβx+G2sinβx+G3coshβx+G4sinhβx⟶wx,t=ψxzt,where zt is the time-domain response and ψx is the eigenvector of the beam. The boundary and initial conditions of the beam are required to determine the coefficients Gi and A, B. According to Figure 1 and equation (7), in the left boundary of the beam, the displacement and slope are zero:(7)w0,t=0∂w∂x0,t=0⟶G1+G3=0G2+G4=0.Based on equation (7), the mode shape is rewritten as follows:(8)ψx=G1cosβx−coshβx+G2sinβx−sinhβx.The boundary conditions of the cantilever beam with a concentrated massM at the free end are as follows:(9)Ml,t=0∂∂xEI∂2wl,t∂x2=M∂2wl,t∂t2⟶EI∂2wl,t∂x2=0EI∂3wl,t∂x3−M∂2wl,t∂t2=0.The resultEIw″l,t=0 and EIw‴l,t=Mw¨l,t could be rewritten as equations (10) and (11), respectively:(10)G1cosβl−coshβl+G2sinβl−sinhβl=0,(11)EIMψ‴l=ψlz¨tzt⟶−EIMω2ψ‴l=ψl.By replacing the valuesψl,ψ″l into equation (11), the following equation is obtained:(12)EIβ3Mω2−G1sinβl+sinhβl+G2cosβl+coshβl=G1cosβl−coshβl+G2sinβl−sinhβl.By replacingλ=βl, α=m¯l/Mλ, and β4/ω2=EI/m¯ in equation (10), equation (12) is expressed as follows:(13)cosλ+coshλsinλ+sinhλα−sinλ+sinhλ−cosλ+coshλαcosλ+coshλ−sinλ+sinhλG1G2=00.To obtain the eigenvalues or frequencies (λ), the determinant of equation (13) should be zero as follows:(14)αcosλ+coshλ2+cosλ+coshλ−sinλ+sinhλ+αsinh2λ−sin2λ+sinλ+sinhλ−cosλ+coshλ=0.The simplified form of equation (14) is expressed as equation (15). By solving this equation, the eigenvalues (λn) are available. This nonlinear equation requires special numerical techniques to obtain all roots of eigenvalues [31]:(15)sinλncoshλn−cosλnsinhλn1+cosλncoshλn−m¯lMλn=0.According to the following equation, the frequency (ωn) in each mode is available:(16)λn2=ωnl2m¯EI⟶ωn=λnl2EIm¯.To obtain the eigenvectors of the beam, one should calculate the coefficientsG1andG2. In this study, G1 is set to the unit, and the general unknown factor Gn is computed from the following equation:(17)Gn=αncosλn+αncoshλn−sinλn+sinhλnαnsinλn−αnsinhλn+cosλn−coshλn,λn=βnl:αn=m¯lMλn=m¯βnM.Therefore, the eigenvectors or equivalently mode shapes of the studied beam in Figure1 is available according to the following equation:(18)ψnx=sinβnx−sinhβnx+Gncosβnx−coshβnx,ψnx=sinλnxl−sinhλnxl+Gncosλnxl−coshλnxl.


## 2.2. Numerical Modal Parameters

The finite element method was used to obtain the modal parameters in the numerical method. Firstly, using the finite element method, the studied beam in Figure1 was divided into twenty-five equal parts with 3 centimeters length. Each element has two nodes, and each node has a transitional degree of freedom (wi) and a rotation degree of freedom (θi). Stiffness matrix Ke and mass matrix Me of each element are determined as follows:(19)Ke=EIl31261−126161412−61212−12−6112−6161212−61412;Me=ρAl42015622154−131221412131−31254131156−221−131−312−221412,where E is the modulus of elasticity, I is the moment of inertia of the cross section, and l is the length of each element (l=L/Nelem), where L is the length of the beam, Nelem is the number of elements, ρ is the mass density of the material, and A refers to the cross section of the beam. It should be noted that the last element of the beam has a concentrated mass (M) at the second node. Therefore, during the assembling process of the total mass matrix, this concentrated mass should be added to the mass matrix as follows:(20)Mcon=0000000000M00000.The beam’s length was equal to 75 cm, and the number of elements was 25. To compare the results and update the finite element model, in the second stage, the number of elements reduced to four elements.


## 2.3. Experimental Modal Parameters

The aluminum beam, as mentioned in Section2.1, was also studied in an experimental study. As shown in Figure 2, the acceleration data were collected by installing unidirectional acceleration sensors at four points of the beam and applying the hammer force at the specified positions.Figure 2
Overall geometry of the beam and the locations of acceleration sensors (a) and hammer (b) in the experimental study.
(a)(b)According to Figure3, the input data (forces) and output data (accelerations) are recorded at positions P1, P2, and P3, and they are collected at 4 stations (CH1, CH2, CH3, and CH4). The entire data collection process was 30 seconds. Each test is repeated three times. Based on the data sampling frequency of the data acquisition system (400 Hz), for each test, 12,000 points are recorded.Figure 3
Location of exerted forces (P1, P2, and P3) and acceleration sensors (CH1, CH2, CH3, and CH4) in the studied beam.Extraction of modal parameters was accomplished using frequency response function (FRF). For example, FRFs; in the case that the hammer exerted at the free end (Input 1), the acceleration data (Outputs 1, 2, 3, and 4) are presented in Figure 4. The required experimental modal data are extracted from FRFs and written as codes in MATLAB [32].Figure 4
The frequency response functions (FRFs) for experimental tests in the studied beam. (a) Output 1, Input 1. (b) Output 2, Input 1. (c) Output 3, Input 1. (d) Output 4, Input 1.
(a)(b)(c)(d)(e)(f)(g)


## 2.4. Updating the Finite Element Model

The aim of updating the finite element model is correcting or improving the finite element model of a structure in terms of updating the global or local mass and stiffness matrices, boundary conditions, and geometrical and material properties. For this purpose, the updating process was accomplished with the modal assurance criterion (MAC) method [5] based on the recorded data at four stations in the experiment program to correcting the global mass and stiffness matrices. As well, the finite element model is reduced to coincide with the experimental setup with the SEREP algorithm [33]. Update has been conducted for mass and stiffness matrices using the particle swarm optimization (PSO) technique [34]. The objective function is defined based on the frequencies and the components of the mode shapes in the experimental and numerical data as follows:(21)f=ωinum−ωiexpωiexp2+∑j=1ndfϕi,jnum−ϕi,jexpϕi,jexp2,where ωinumandϕi,jnum are the frequencies and the jth components of mode shapes of the finite element model, respectively, and ωiexpandϕi,jexp are the frequencies and the jth mode shapes of the experimental data, respectively.


## 3. Results

### 3.1. Frequencies and the Mode Shapes

Modal parameters from analytical, numerical, and also experimental methods were extracted and presented in this section. Table1 summarizes the geometrical and material properties of the studied aluminum beam.Table 1
Geometrical and material properties of the beam.
ρ (kg/m3)b (mm)h (mm)L (mm)m¯ (kg/m)M (kg)E (GPa)270040107501.081.58569In the numerical method, firstly, the beam is divided into 25 elements; after that, frequencies and mode shapes were extracted. Then, to compare the numerical and experimental results due to the limitations of acceleration sensors, the numerical model was also reduced from 25 elements to four elements. The results of the frequencies are presented in Table2. As well, in Figure 5, the mode shapes of the first four modes are illustrated with experimental [35], analytical, and also numerical methods.Table 2
First four natural frequencies (rad/s) with analytical, numerical, and experimental methods.
Mode numberAnalytical frequenciesNumerical frequencies (25 elements)Numerical frequencies (4 elements)Experimental frequencies130.3430.3430.3429.342411.78411.72411.98341.1231309.291309.061317.161151.4142715.072717.522777.812462.11Figure 5
The first four mode shapes of the beam in the analytical, numerical (4 elements), and experimental methods. (a) Mode no. 1. (b) Mode no. 2. (c) Mode no. 3. (d) Mode no. 4.
(a)(b)(c)(d)As can be seen from Figure5, the analytical and numerical mode shapes of the beam have good agreements. However, the experimental results in some modes differ from the other methods. The reason for the differences in some modes is related to the scaling of mode shapes. Therefore, the numerical model has been updated to modify the mass and stiffness matrices of the finite element model [5], which is performed with the particle swarm algorithm based on the objective function defined in equation (21). The convergence of the objective is shown in Figure 6. As well, Table 3 presents the values of the numerical frequencies before and after the updating process. In this study, numerical model updating was conducted only with frequency information.Figure 6
Convergence of the cost function for the numerical model updating with the PSO algorithm.Table 3
Numerical frequencies of the studied beam before and after updating.
Mode numberNumerical frequencies (before updating)Numerical frequencies (after updating)130.3428.802411.98344.8331317.161176.9242777.812404.61The updated values of mass and stiffness matrices are summarized in Table4. These values are applied to calculate the numerical responses of the beam under impact loading.Table 4
Mass and stiffness matrices of the beam before and after updating based on frequency data.
Kupdate=585230−343390131300−27830−343390387810−25951074370131300−259510259630−94480−2783074370−9448039910K=702420−456600182920−33420−456600546340−36461099230182920−364610366120−133460−3342099230−1334605750Μupdate=0.1820.014−0.010.0110.0140.1820.019−0.017−0.010.0190.1860.0320.011−0.0170.0321.643M=0.1750.015−0.0070.010.0150.1890.008−0.017−0.0070.0080.1960.0350.01−0.0170.0351.639Unit of mass is kg and stiffness is N/mAfter calculating the updated mass and stiffness matrices, the damping matrix is obtained using the Rayleigh damping matrix as follows:(22)C=αMupdate+βKupdate.The coefficientsα and β are obtained, according to equation (23), where ωj and ωi are the first and last available vibration frequencies, respectively, and also, ξ is a critical damping ratio that was obtained equal to 0.4% from the experimental data:(23)α=2ξωiωjωi+ωj,β=2ξωi+ωj.


### 3.2. Dynamic Analysis Using the Newmark Method (Harmonic Loading)

The unconditionally stable Newmark average acceleration method is recommended for the calculation of the response by modal superposition when filtering the high frequencies is naturally carried out by using a truncated modal space. It is the most accurate method, and the absence of numerical damping is not a problem [36]. In this paper, the capabilities of the Newmark method with the direct integration method, instead of the modal superposition method, are investigated by studying an SDF system. For this purpose, consider the damped SDF model with the specification that is summarized in Table 5. The natural frequency of the SDF system is obtained equal to 12 rad/s.Table 5
Specification of a damped SDF system.
m (kg)k (kN/m)ξω¯ (rad/s)p(t) (N)30004320.1%, 0.5%1.2, 60−mgsinω¯tThe structural response of the SDF system is presented in Figure7 for harmonic loading with a frequency equal to 1.2 rad/s and damping ratio (ξ) equal to 0.4%. As seen in this figure, the numerical method has excellent capability to predict the response. In this case, the natural frequency of the SDF system is 12 rad/s, and from the structural dynamics point of view, an excitation with ω¯ = 1.2 rad/s could be taken into account as a forced vibration-type loading (Figure 7).Figure 7
Response of the beam with the Newmark and analytical methods forξ=0.4% and ω¯ = 1.2 rad/s.The structural response of the SDF system is presented in Figure8 for harmonic loading with a frequency equal to 60 rad/s and damping ratio (ξ) equal to 0.5%. As seen in this figure, the numerical method has good capability to predict the responses. From structural dynamics points of view, an excitation with ω¯ = 60 rad/s and ω = 12 rad/s could be taken into account as a free vibration.Figure 8
Response of the beam with the Newmark and analytical methods forξ=0.4% and ω¯ = 60 rad/s.Figures7 and 8 reveal that for regular loading, the Newmark average acceleration method has an appropriate performance. The free end displacement responses of the studied beam were investigated to show the abilities of the Newmark method under impact or very short-time loadings, and the results were compared with the experimental data.


### 3.3. Dynamic Analysis Using the Newmark Method (Impact Loading)

According to Figure9, the studied beam is subjected to very short-time impact loading at its free end in the laboratory. As seen from Figure 10, the acceleration data at four stations (CH1, CH2, CH3, and CH4) are presented.Figure 9
History of the exerted impact load by using the hammer at the free end of the beam (experimental data).Figure 10
The history of experimental acceleration at four stations when the hammer is exerted at the free end of the beam. (a) Acceleration CH1. (b) Acceleration CH2. (c) Acceleration CH3. (d) Acceleration CH4.The Newmark-beta average acceleration technique [36] was used to determine the numerical response of the beam, and the results are compared with recorded experimental data at stations 1–4 (CH1 to CH4). As shown in Figure 11, the numerical and experimental values have no appropriate agreements even after updating the finite element model. The only possibility that comes to the authors’ mind is that the Newmark is not suitable for shock-type loadings, and the dynamic responses of the low damping MDF structures should be obtained from other methods [37, 38].Figure 11
Comparison of numerical and experimental accelerations of the beam at stations 1–4 (CH1, CH2, CH3, and CH4) under impact loading at the free end.


## 3.1. Frequencies and the Mode Shapes

Modal parameters from analytical, numerical, and also experimental methods were extracted and presented in this section. Table1 summarizes the geometrical and material properties of the studied aluminum beam.Table 1
Geometrical and material properties of the beam.
ρ (kg/m3)b (mm)h (mm)L (mm)m¯ (kg/m)M (kg)E (GPa)270040107501.081.58569In the numerical method, firstly, the beam is divided into 25 elements; after that, frequencies and mode shapes were extracted. Then, to compare the numerical and experimental results due to the limitations of acceleration sensors, the numerical model was also reduced from 25 elements to four elements. The results of the frequencies are presented in Table2. As well, in Figure 5, the mode shapes of the first four modes are illustrated with experimental [35], analytical, and also numerical methods.Table 2
First four natural frequencies (rad/s) with analytical, numerical, and experimental methods.
Mode numberAnalytical frequenciesNumerical frequencies (25 elements)Numerical frequencies (4 elements)Experimental frequencies130.3430.3430.3429.342411.78411.72411.98341.1231309.291309.061317.161151.4142715.072717.522777.812462.11Figure 5
The first four mode shapes of the beam in the analytical, numerical (4 elements), and experimental methods. (a) Mode no. 1. (b) Mode no. 2. (c) Mode no. 3. (d) Mode no. 4.
(a)(b)(c)(d)As can be seen from Figure5, the analytical and numerical mode shapes of the beam have good agreements. However, the experimental results in some modes differ from the other methods. The reason for the differences in some modes is related to the scaling of mode shapes. Therefore, the numerical model has been updated to modify the mass and stiffness matrices of the finite element model [5], which is performed with the particle swarm algorithm based on the objective function defined in equation (21). The convergence of the objective is shown in Figure 6. As well, Table 3 presents the values of the numerical frequencies before and after the updating process. In this study, numerical model updating was conducted only with frequency information.Figure 6
Convergence of the cost function for the numerical model updating with the PSO algorithm.Table 3
Numerical frequencies of the studied beam before and after updating.
Mode numberNumerical frequencies (before updating)Numerical frequencies (after updating)130.3428.802411.98344.8331317.161176.9242777.812404.61The updated values of mass and stiffness matrices are summarized in Table4. These values are applied to calculate the numerical responses of the beam under impact loading.Table 4
Mass and stiffness matrices of the beam before and after updating based on frequency data.
Kupdate=585230−343390131300−27830−343390387810−25951074370131300−259510259630−94480−2783074370−9448039910K=702420−456600182920−33420−456600546340−36461099230182920−364610366120−133460−3342099230−1334605750Μupdate=0.1820.014−0.010.0110.0140.1820.019−0.017−0.010.0190.1860.0320.011−0.0170.0321.643M=0.1750.015−0.0070.010.0150.1890.008−0.017−0.0070.0080.1960.0350.01−0.0170.0351.639Unit of mass is kg and stiffness is N/mAfter calculating the updated mass and stiffness matrices, the damping matrix is obtained using the Rayleigh damping matrix as follows:(22)C=αMupdate+βKupdate.The coefficientsα and β are obtained, according to equation (23), where ωj and ωi are the first and last available vibration frequencies, respectively, and also, ξ is a critical damping ratio that was obtained equal to 0.4% from the experimental data:(23)α=2ξωiωjωi+ωj,β=2ξωi+ωj.


## 3.2. Dynamic Analysis Using the Newmark Method (Harmonic Loading)

The unconditionally stable Newmark average acceleration method is recommended for the calculation of the response by modal superposition when filtering the high frequencies is naturally carried out by using a truncated modal space. It is the most accurate method, and the absence of numerical damping is not a problem [36]. In this paper, the capabilities of the Newmark method with the direct integration method, instead of the modal superposition method, are investigated by studying an SDF system. For this purpose, consider the damped SDF model with the specification that is summarized in Table 5. The natural frequency of the SDF system is obtained equal to 12 rad/s.Table 5
Specification of a damped SDF system.
m (kg)k (kN/m)ξω¯ (rad/s)p(t) (N)30004320.1%, 0.5%1.2, 60−mgsinω¯tThe structural response of the SDF system is presented in Figure7 for harmonic loading with a frequency equal to 1.2 rad/s and damping ratio (ξ) equal to 0.4%. As seen in this figure, the numerical method has excellent capability to predict the response. In this case, the natural frequency of the SDF system is 12 rad/s, and from the structural dynamics point of view, an excitation with ω¯ = 1.2 rad/s could be taken into account as a forced vibration-type loading (Figure 7).Figure 7
Response of the beam with the Newmark and analytical methods forξ=0.4% and ω¯ = 1.2 rad/s.The structural response of the SDF system is presented in Figure8 for harmonic loading with a frequency equal to 60 rad/s and damping ratio (ξ) equal to 0.5%. As seen in this figure, the numerical method has good capability to predict the responses. From structural dynamics points of view, an excitation with ω¯ = 60 rad/s and ω = 12 rad/s could be taken into account as a free vibration.Figure 8
Response of the beam with the Newmark and analytical methods forξ=0.4% and ω¯ = 60 rad/s.Figures7 and 8 reveal that for regular loading, the Newmark average acceleration method has an appropriate performance. The free end displacement responses of the studied beam were investigated to show the abilities of the Newmark method under impact or very short-time loadings, and the results were compared with the experimental data.


## 3.3. Dynamic Analysis Using the Newmark Method (Impact Loading)

According to Figure9, the studied beam is subjected to very short-time impact loading at its free end in the laboratory. As seen from Figure 10, the acceleration data at four stations (CH1, CH2, CH3, and CH4) are presented.Figure 9
History of the exerted impact load by using the hammer at the free end of the beam (experimental data).Figure 10
The history of experimental acceleration at four stations when the hammer is exerted at the free end of the beam. (a) Acceleration CH1. (b) Acceleration CH2. (c) Acceleration CH3. (d) Acceleration CH4.The Newmark-beta average acceleration technique [36] was used to determine the numerical response of the beam, and the results are compared with recorded experimental data at stations 1–4 (CH1 to CH4). As shown in Figure 11, the numerical and experimental values have no appropriate agreements even after updating the finite element model. The only possibility that comes to the authors’ mind is that the Newmark is not suitable for shock-type loadings, and the dynamic responses of the low damping MDF structures should be obtained from other methods [37, 38].Figure 11
Comparison of numerical and experimental accelerations of the beam at stations 1–4 (CH1, CH2, CH3, and CH4) under impact loading at the free end.


## 4. Conclusions

In this study, the modal parameters of the aluminum beam were determined by analytical, numerical, and also experimental methods. Then, updating the finite element model is conducted based on the modal data and particle swarm optimization techniques. Then, the updated FE model was used to estimate the response of the beam under a very short-time impact load. Based on the results of this study, the following conclusions could be drawn:(i)
In the numerical model (finite element), increasing the number of DOF leads to an accurate estimation of higher modes but does not improve the low mode values. As well, because of the technical limitations in measurement of the rotations, using model reduction techniques to eliminate the rotational DOF does not lead to a reduction of accuracy of the numerical model.(ii)
The results show that if the objective function is only defined based on frequency data, the accuracy of the model is slightly reduced. However, simultaneous use of mode shapes and frequencies data improve the accuracy, but extracting the mode shapes from the laboratory data is not an easy task.(iii)
The results reveal that the Newmark average acceleration method has excellent capabilities to accurately predict the dynamic responses under regular type loading even for low damping systems without facing numerical damping error.(iv)
The results show that the Newmark average acceleration method that inherently is unconditionally stable is not able to accurately estimate the dynamic responses under shock- or impact-type loading, and other numerical or analytical methods should be utilized.



---

*Source: 8861827-2020-09-15.xml*