A Damage Detection Algorithm Utilizing Dynamic Displacement of Bridge under Moving Vehicle

Abstract

A damage detection method is proposed, which utilizes dynamic displacement of bridge structures under moving vehicle. The problem is first elaborated with closed-form solution of dynamic displacement, which is decomposed into quasi-static component and dynamic component. Dynamic curvature is defined as second derivative of the dynamic displacement for detecting damage location and estimating damage extent. Damage is modeled by local reduction of stiffness in this paper. Numerical study was conducted on a simply supported beam to verify the proposed method. Vehicle model is analyzed with Newmark's method using Matlab to obtain the contact force acting on the bridge. Beam model is established in commercial finite element software ABAQUS. The effects of road surface roughness and vehicle-bridge interaction are both considered in the analysis. In order to identify damage location and extent, dynamic curvature was calculated with midspan displacement. Parametric study on measurement noise level, damage location, damage extent, and multiple damage cases is performed, and the analysis results show both reliability and efficacy of this method in damage detection of bridge structures. At last, conclusions are drawn for its application to bridges in engineering practice.

Zhen Sun 1,2 and Tomonori Nagayama 1 and Di Su 1 and Yozo Fujino 3

Academic Editor:Matteo Aureli

1, Department of Civil Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
2, Jiangsu Transportation Institute, 2200 Chengxin Street, Nanjing 211112, China
3, Institute of Advanced Sciences, Yokohama National University, 79-1 Tokiwadai, Hodogaya-ku, Yokohama 240-8501, Japan

Received 30 July 2015; Revised 21 September 2015; Accepted 4 October 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

During the service life of a bridge, factors such as material aging and harsh environmental conditions as well as extreme events such as natural disasters may cause structural damage to varying degrees [1]. In order to avoid severe accident of bridge collapse, detection of any potential damage is important for serviceability and safety of bridges.

In recent years, vibration based methods have been widely studied for bridge assessment. Detailed review of the damage detection methods can be found in [2, 3]. Most of these methods utilize mode frequency or mode shapes as the identification parameter. However, there are several difficulties for these frequency domain methods: (1) baseline of previous test data of undamaged structure is needed, which is often unavailable; (2) frequency, as a global modal parameter, is not sensitive to local damage; (3) frequency is often affected by environmental conditions such as temperature or humidity. These difficulties have prevented practical use of vibration based methods.

Besides the frequency domain methods, some efforts are devoted to investigating the problem with time domain damage detection methods. In the time domain methods, a moving vehicle is often used as external excitation to obtain the dynamic response of the bridge. The basic concept is that the moving vehicle will pass every section of the bridge, which includes the damage part. When the vehicle is right over the damage range, singularity will be caused in the response compared with other parts. Through data-processing techniques to find singularity in the signal, damage location could be detected. Huang [4] applied Hilbert-Huang Transform (HHT) to a model of simply supported damaged bridge under loading of a moving truck. The undamaged model is taken as baseline, and variation of local frequency and deformation is employed to detect damage. Some other researchers used wavelet analysis to detect singularity in bridge response signal, hence detecting damage in bridges. Zhu and Law [5] conducted wavelet analysis on bridge response under a moving load. Wavelet coefficient of dynamic displacement is obtained to detect damage in a bridge model. Nguyen and Tran [6] conducted similar study, but, instead of using bridge response, they used vehicle response for wavelet analysis, which makes it unnecessary for sensor installation on a bridge. Hester and Gonzalez [7] performed wavelet transform on response of a damaged bridge under moving load and investigated feasibility to detect damage location. They concluded that multiple measurement points were required to achieve higher accuracy. In the wavelet-based damage detection techniques, there is a major inherent difficulty. In order to localize damage, narrow window of bridge response, that is, short time series, is needed. However, frequency resolution is impaired when short time series are chosen to calculate the frequency. Due to this contradiction, effective detection of damage is difficult.

In this paper, a curvature based damage detection method is proposed, which uses dynamic displacement of a bridge subject to loading of a moving vehicle. It does not require baseline of undamaged state of bridge, and both damage location and damage extent can be detected. Firstly, the physical mechanism of dynamic displacement of damaged bridge under loading of a moving vehicle is elaborated. Dynamic displacement is decomposed into quasi-static component and dynamic component, and curvature is used to detect damage location and estimate damage extent. Secondly, the proposed method is verified with a finite element model using beam elements, in which vehicle model and beam model are established in Matlab and ABAQUS, respectively. Vehicle-bridge interaction and roughness are both considered, and dynamic displacement is used to calculate curvature and detect damage. The effect of parameters of measurement noise level, damage location, damage extent, and multiple damage cases is studied, which illustrates reliability and efficacy of the proposed method. Finally, conclusions are made for application of this method in engineering practice.

2. Damage Detection with Dynamic Displacement

The dynamic problem of vehicle-bridge system has been a subject of interest for a long time. In the earlier times, most of the studies are aimed at predicting dynamic response of bridges and providing a dynamic amplification factor for engineering design [8]. Detailed review of the early analytical and numerical studies could be found in previous literature such as [9, 10]. In the following section, physical background is elaborated, and it is shown that dynamic curvature is a good indicator for damage detection.

2.1. Analysis on Dynamic Displacement under Moving Vehicle

In simplified case of an undamaged uniform simply supported beam, vehicle is modeled as a moving load with constant speed v , as illustrated in Figure 1. The closed-form solution was deduced by Weaver et al. [11], and the displacement of the bridge is expressed as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] ; [figure omitted; refer to PDF] ; [figure omitted; refer to PDF] is the length of the bridge; [figure omitted; refer to PDF] is the amplitude of the moving load; [figure omitted; refer to PDF] is the mass per unit length; [figure omitted; refer to PDF] is the flexural rigidity of the bridge; [figure omitted; refer to PDF] is the velocity of the moving force; [figure omitted; refer to PDF] is the natural frequency of the bridge.

Figure 1: Schematic figure of simply supported beam under moving load.

[figure omitted; refer to PDF]

It is noticed that the first part is related to the vehicle speed, which is referred to as the quasi-static component; the second part contains bridge natural frequency and is referred to as dynamic component, which comes from the vibration of the bridge. Figure 2 illustrates the dynamic displacement at midspan, which is composed of these two parts: quasi-static component and dynamic component. The proposed method concentrates on calculating curvature from the quasi-static component of the dynamic displacement, and an innovative technique is used in curvature calculation to cancel out the dynamic displacement.

Figure 2: Quasi-static and dynamic components in deflection signal.

[figure omitted; refer to PDF]

2.2. Damage Detection with Curvature

The relation between bending moment and curvature of a beam is expressed in the following equation [12]: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] denotes Young's modulus of the material, [figure omitted; refer to PDF] denotes the second moment of area of beam cross section, [figure omitted; refer to PDF] denotes deflection of the beam, [figure omitted; refer to PDF] denotes the beam coordinate, and [figure omitted; refer to PDF] denotes the bending moment.

Curvature can be calculated from the second derivative of deflection measurements, which is expressed as [figure omitted; refer to PDF] If there is damage (stiffness loss) in some location of the beam, abrupt change is expected at the corresponding location on the curvature plot.

In this study, damage extent is defined as [figure omitted; refer to PDF] , which indicates loss in stiffness [figure omitted; refer to PDF] : [figure omitted; refer to PDF] where the subscripts [figure omitted; refer to PDF] and [figure omitted; refer to PDF] denote "intact" and "damaged," respectively.

The stiffness in intact and damaged states is defined, respectively, as [figure omitted; refer to PDF] Defining [figure omitted; refer to PDF] and substituting (5) and (6) into (4) yields [figure omitted; refer to PDF] Hence, damage extent can be calculated as [figure omitted; refer to PDF] Curvature can be approximated with second derivative of displacement measurements: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] denotes displacement. Illustration of the proposed method is shown in Figure 3. In order to obtain dense deflection measurement, Maxwell-Betti Reciprocal Theorem is utilized and it is explained as follows. Consider a beam in two states under different loads. Maxwell-Betti Reciprocal Theorem states that the work of the forces of the first state on the displacements of the second state is equal to the work of forces of the second state on the displacements of the first state. Based on this theorem, the deflection of a beam at point A under a load at point B is the same as the deflection at point B under the same load at point A. By extending this equivalence relationship to a moving load, the static deflection w(x ) of a beam under a load at point A, P (A ), is obtained with deflection measurement at point A, w (A ), under a moving load P (x ). The load is assumed to move slowly so that the dynamic effect is negligible. In application to engineering practice, the vehicle load is assumed to be point load, and the bridge is modeled as a beam.

Figure 3: Illustration of damage identification from dynamic displacement.

[figure omitted; refer to PDF]

For dynamic displacement of a damaged bridge under vehicle loading, it contains three parts: quasi-static component, dynamic component, and damage-induced component. Therefore (9) can be expressed as [figure omitted; refer to PDF] Noticing that dynamic component is a periodic function (Figure 2), in curvature calculation formula, calculation interval [figure omitted; refer to PDF] can be chosen as the length of the natural period: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] denotes fundamental natural period of the bridge; [figure omitted; refer to PDF] denotes sampling frequency of data acquisition. In this way, the three terms in the numerator of dynamic component are canceled out, while contribution from damage is still kept in the calculated curvature. For a beam, the fundamental natural period is integer times of higher bending mode periods. Therefore, higher bending mode components are also canceled out by the option of [figure omitted; refer to PDF] .

In this analysis, [figure omitted; refer to PDF] is dynamic displacement instead of static deflection, and [figure omitted; refer to PDF] is chosen as the length of fundamental natural period. [figure omitted; refer to PDF] is defined as dynamic curvature to distinguish with the curvature calculated with static deflection.

In this proposed method, as shown in Figure 3, a displacement transducer is installed at one position, while a vehicle passes through a bridge. If there is discontinuity in recorded dynamic displacement signal, it will correspond to damage on the bridge. Dynamic curvature will highlight the damage location, and damage extent can also be calculated from curvature values.

3. Investigation with a Finite Element Model

In order to verify the proposed damage detection method, investigation is conducted on a vehicle-bridge model using Matlab and finite element software ABAQUS. In the analysis, the vehicle model is established using Matlab, and the contact force at the wheel is calculated with Newmark's method. The bridge is modeled as a simply supported beam in ABAQUS, and contact force is imported from Matlab calculation result. The analysis scheme can be used to compute the bridge behavior with good accuracy and high efficiency.

3.1. Numerical Example

3.1.1. Vehicle Model

Equation of motion for the vehicle is obtained using the principle of virtual work. The vehicle loading is simulated with a quarter-car model shown in Figure 4 composed of a spring with stiffness coefficient of [figure omitted; refer to PDF] N/m and a damper with damping coefficient of [figure omitted; refer to PDF] N·s/m. The mass of vehicle body is 100 kg. Considering the mass of tire is much smaller than the vehicle body, it is neglected in the calculation.

Figure 4: Quarter-car model.

[figure omitted; refer to PDF]

The equation of motion for the vehicle is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] denote mass, damping coefficient, and stiffness of the vehicle, respectively; [figure omitted; refer to PDF] denote acceleration, velocity, and displacement of vehicle motion, respectively; [figure omitted; refer to PDF] denotes the vertical force on the vehicle, which is contributed from road roughness function and bridge displacement at the contact points. Constant vehicle speed and straight vehicle path are assumed in this study. The dynamic equation is calculated in Matlab with Newmark's method, and the contact force on the bridge is calculated as follows: [figure omitted; refer to PDF]

3.1.2. Bridge Model

Steel beam model is established with three-dimensional beam element B31 in commercial finite element software ABAQUS shown in Figure 5. It is 15 m long and 1 m wide, with a thickness of 0.06 m, and the boundary condition is specified as simply supported. The beam is modeled totally with 200 beam elements.

Figure 5: Layout of bridge model.

[figure omitted; refer to PDF]

The bridge specifications are provided in Table 1. Damage is simulated by reduction of cross section, which is 11 m away from the left end, with a length of 0.75 m, and the thickness reduced to 0.05 m. Damage extent is 42% stiffness loss. In analysis, vehicle passes at a constant speed of 1 m/s, so totally it takes 15 s to pass the bridge. Rayleigh damping is defined for the model.

Table 1: Properties of the beam model.

Density	7800 kg/m3
Young's modulus	[figure omitted; refer to PDF] N/m2
Poisson's ratio	0.3
Damage location	11 m from the left support

Road surface is continuously distributed in a random fashion, which affects dynamic behavior of both vehicle and bridge. It is usually assumed to be a stationary Gaussian random process [13, 14] and is generated from power spectral density functions. Consider [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is power spectral density function and [figure omitted; refer to PDF] denotes road condition and is chosen as [figure omitted; refer to PDF] ; [figure omitted; refer to PDF] is the frequency interval; [figure omitted; refer to PDF] ; [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the upper and lower cut-off frequencies, respectively; here they are chosen as 5 Hz and 0.1 Hz; [figure omitted; refer to PDF] is the wave number; [figure omitted; refer to PDF] is the random phase angle uniformly distributed from 0 to [figure omitted; refer to PDF] .

When calculating the dynamic behavior of the vehicle, the space function is transformed to time function according to the vehicle speed. The generated road profile is shown in Figure 6, with the largest amplitude of roughness below 2 mm.

Figure 6: Bridge road surface profile.

[figure omitted; refer to PDF]

3.1.3. Calculation Algorithm

In this analysis, interaction between the bridge and the vehicle is realized through contact force. At the contact points between the vehicle wheel and the bridge, displacement of the vehicle and the bridge should satisfy compatibility at the contact points, where the same deformation of both the vehicle and the bridge should be maintained. Convergence should be satisfied at every contact point.

There has been step by step algorithm for the computation [8], in which the vehicle moves through the bridge, and the contact force is calculated from the vehicle behavior at each contact point. Then the contact force is input to the bridge, which causes bridge response and the response is fed back to the vehicle model. After that the vehicle moves to the next point, and the process continues until the vehicle passes the bridge. In addition to step by step algorithm, there is also iterative procedure which treats the contact force as time history [15]. In this analysis, computation procedure is conducted in a global manner presented in Figure 7, and interaction between the vehicle and the bridge is realized through the contact force time history. First, road profile is used as the initial input for the vehicle, and the contact force history is calculated for each contact point. The contact force history is then applied to the bridge model in ABAQUS, which leads to bridge response. Then the bridge response is incorporated in the vehicle dynamic equation as input, together with the road profile, and a new contact force time history is generated. This procedure is repeated until user defined convergence criteria are satisfied. Particularly the convergence criterion is set that the maximum deflection difference between successive iterations is smaller than 0.01% of the dynamic displacement at midspan. Once the convergence criterion is satisfied, displacement response of the bridge is extracted for damage detection.

Figure 7: Calculation algorithm for vehicle-bridge model.

[figure omitted; refer to PDF]

3.1.4. Analysis Result

Analysis is conducted until convergence is satisfied. The difference of dynamic displacement at the midspan point is shown in Figure 8. It indicates that, after three iterations, the maximum difference between two iterations is [figure omitted; refer to PDF] m, around 0.0025% of the midspan deflection of the bridge. It satisfies the convergence criterion that the displacement difference is smaller than 0.01% of the midspan deflection.

Figure 8: Convergence of displacement at contact points.

[figure omitted; refer to PDF]

Dynamic displacement at the midspan point is extracted for damage detection given in Figure 9. The natural frequencies of the bridge were also extracted, which shows that the fundamental natural frequency is 0.60765 Hz. Hence, the fundamental natural period is obtained as 1.648 s. In data processing, the second derivatives of the dynamic displacement data are taken at an interval of the natural period of the bridge. The vehicle passes the bridge in 15 s, and displacement is recorded. Therefore, [figure omitted; refer to PDF] is calculated as 13.4 Hz. For the calculation of dynamic curvature, [figure omitted; refer to PDF] is chosen as [figure omitted; refer to PDF] . [figure omitted; refer to PDF] is then rounded to an integer value, which is 22.

Figure 9: Dynamic displacement at the center of the bridge span.

[figure omitted; refer to PDF]

The dynamic curvature is illustrated in Figure 10. When the vehicle-bridge interaction is not considered, the dynamic curvature is also presented for comparison, given by the red dashed line. In the curvature figure, there are no calculated curvature values at both ends of the bridge span, because the second derivative was taken at an interval [figure omitted; refer to PDF] . It can be indicated from the figure that when interaction is not considered, curvature plot is smoother, and peak caused by damage is more apparent. In contrast, there are more fluctuations on the curvature result which incorporated interaction between the vehicle and the bridge. But, even considering the interaction effect, the damage location is clearly identified.

Figure 10: Curvature of bridge at midpoint of the span.

[figure omitted; refer to PDF]

Apparent peak at the location of 11 m on the dynamic curvature plot indicates the damage location. Using the damage extent calculation equation (9), the damage extent is obtained as 39.98%. Compared with the actual damage extent of 42%, the difference is around 2%. It demonstrates that, even with the effects from road roughness and vehicle-bridge interaction, the damage detection method can detect damage location and identify damage extent with good accuracy.

3.2. Effect of Measurement Noise

In order to investigate the effect of measurement noise on efficacy of the proposed method, random noise is generated and added to the calculated displacement response of the bridge. The displacement with measurement noise is calculated as [5, 6] [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the polluted bridge displacement with measurement noise. [figure omitted; refer to PDF] is the noise level and [figure omitted; refer to PDF] is a standard normal distribution with zero mean value and unit standard deviation. [figure omitted; refer to PDF] is the calculated displacement at midpoint of bridge span, and [figure omitted; refer to PDF] is its standard deviation.

Three different noise levels of 0.1%, 0.3%, and 0.5% are investigated, and curvature result is shown in Figure 11. It is shown that damage location can be successfully determined for all three cases. However, damage extent is less accurate, which is identified as 31.3%, 30.1%, and 28.2%, respectively. This is due to the fact that curvature values are distorted by measurement noise. The proposed method can both successfully locate the damage and provide a quantitative estimation of the damage extent.

Figure 11: Curvature at midspan with different noise level.

[figure omitted; refer to PDF]

3.3. Effect of Damage Location and Extent

In order to investigate the effect of damage location and extent, different damage locations and extents are analyzed. Damage is assumed at location of 3 m (L/5), 5 m (L/3), and 12 m (4L/5), respectively, and damage length is 0.75 m with thickness reduced to 0.05 m. Noise level is assumed as 0.1%, and it is added to calculated displacement in the same way as in the last subsection. All other conditions are defined the same as in Section 3.1. The curvature results are illustrated in Figure 12. It is shown that damage can be detected reliably regardless of its location. Damage extent is identified as 36.8%, 32.9%, and 33.9%, respectively.

Figure 12: Curvature at midspan with different damage location.

[figure omitted; refer to PDF]

Effect of damage extent is also studied, by assuming different damage extent of 20%, 40%, and 60%. Noise level is assumed as 0.1%, and all other conditions are defined the same as in Section 3.1. The curvature results are plotted in Figure 13, which shows that the peak caused by damage increases as damage extent increases from 20% to 60%. Damage extent is identified as 17.9%, 32.2%, and 49.6%. The discrepancy can be attributed to the road surface roughness and measurement noise.

Figure 13: Curvature at midspan with different damage extent.

[figure omitted; refer to PDF]

3.4. Effect of Multi-Damage Case

In this subsection, the beam with multiple damage is studied. The bridge model is the same as in Section 3.1 except that multiple damage cases are artificially created. Totally damage is generated at three locations, which are 5 m (1/3L), 10 m (2/3L), and 12 m (4/5L) of the bridge. The damage extent is 42% for all cases, and 0.1% noise level is added to calculated displacement.

Figure 14 shows results of the curvature at midspan of the bridge. It is observed that all three damage cases can be detected at corresponding locations. Damage extent is also identified using (9), which is 33.7%, 34.8%, and 32.3%, respectively.

Figure 14: Curvature at midspan with multiple damage cases.

[figure omitted; refer to PDF]

A finite element bridge model is used to investigate reliability and efficacy of the proposed method. Effects of measurement noise, damage location and extent, and a multi-damage case are studied, which show that the proposed method can be used to detect damage location and estimate damage extent successfully. It is noted that the purpose of this paper is to verify the proposed damage detection algorithm numerically, and a simple supported beam model is adopted. Further application to real bridges should take into account high local modes and torsion modes.

4. Conclusions

Safety of bridges during operation is of much concern for the society, which requires more reliable and efficient damage detection methods. In this paper, a damage detection method is proposed based on dynamic displacement of bridge under the loading of a moving vehicle. Through this study, several conclusions are drawn as follows:

(1) Closed-form solution of the vehicle-bridge problem shed light on decomposing dynamic displacement into quasi-static component and dynamic component, and curvature of the quasi-static component can be used to detect damage location and estimate damage extent.

(2) In curvature calculation, using dynamic displacement, the calculation interval is chosen as natural period of the beam, by which dynamic component can be canceled out, while damage is detected with quasi-static component.

(3) Numerical simulation of a beam model is investigated, in which road surface roughness and vehicle-bridge interaction are considered. Parametric study on measurement noise, damage location, damage extent, and multiple damage cases is performed, which verified the reliability and efficacy of this method.

In light of damage detection of bridges, this method utilizes dynamic displacement of bridges under moving vehicle loading, which can be operated with low cost and simple implementation. Moreover, baseline from undamaged state is not required, and both damage location and damage extent can be assessed. In field test measurements of real bridges, there could be higher local modes or torsional modes, which should be considered in application of the proposed method. Further investigations will be conducted to study these issues in future.

Acknowledgments

The authors would like to express gratitude to Professor Dionysius from Yokohama National University for his advices on this research. The first author is also grateful for the scholarship from MEXT (Ministry of Education, Culture, Sports, Science and Technology, Japan) during his Ph.D. study.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] Y. Fujino, D. M. Siringoringo, "Bridge monitoring in Japan: the needs and strategies," Structure and Infrastructure Engineering , vol. 7, no. 7-8, pp. 597-611, 2011.

[2] O. S. Salawu, "Detection of structural damage through changes in frequency: a review," Engineering Structures , vol. 19, no. 9, pp. 718-723, 1997.

[3] S. W. Doebling, C. R. Farrar, M. B. Prime, "A summary review of vibration-based damage identification methods," The Shock and Vibration Digest , vol. 30, pp. 91-105, 1998.

[4] K. Huang US 6192758 B1 United States, 2001

[5] X. Q. Zhu, S. S. Law, "Wavelet-based crack identification of bridge beam from operational deflection time history," International Journal of Solids and Structures , vol. 43, no. 7-8, pp. 2299-2317, 2006.

[6] K. V. Nguyen, H. T. Tran, "Multi-cracks detection of a beam-like structure based on the on-vehicle vibration signal and wavelet analysis," Journal of Sound and Vibration , vol. 329, no. 21, pp. 4455-4465, 2010.

[7] D. Hester, A. Gonzalez, "A wavelet-based damage detection algorithm based on bridge acceleration response to a vehicle," Mechanical Systems and Signal Processing , vol. 28, pp. 145-166, 2012.

[8] J.-C. Wyss, D. Su, Y. Fujino, "Prediction of vehicle-induced local responses and application to a skewed girder bridge," Engineering Structures , vol. 33, no. 4, pp. 1088-1097, 2011.

[9] S. P. Timoshenko History of Strength of Materials , D. Van Nostrand Company, New York, NY, USA, 1953.

[10] L. Fryba Vibration of Solids and Structures Under Moving Loads , Academia Prague, Prague, Czech Republic, 1999.

[11] J. R. W. Weaver, S. P. Timoshenko, D. H. Young Vibration Problems in Engineering , John Wiley & Sons, New York, NY, USA, 1990.

[12] S. P. Timoshenko Strength of Materials , D. Van Nostrand Company, New York, NY, USA, 1955.

[13] C. J. Dodds, J. D. Robson, "The description of road surface roughness," Journal of Sound and Vibration , vol. 31, no. 2, pp. 175-183, 1973.

[14] Y.-B. Yang, B.-H. Lin, "Vehicle-bridge interaction analysis by dynamic condensation method," Journal of Structural Engineering , vol. 121, no. 11, pp. 1636-1643, 1995.

[15] M. David Finite Element Analysis , Sciyo, 2010.