Nonlinear Dynamic Analysis of Plates Subjected to Explosive Loads

Abstract

Blast loads have been increasingly studied in the past decades, especially regarding civil structures. Until recently, the negative phase of these loads has been disregarded, but studies concluded that the effect of suction must be included. In the case of plates, nonlinearity plays an important role, and the membrane effect should also be considered. This work focuses on the influence of nonlinearity in plates subjected to blast loads. Equations were developed for the calculation of blast load parameters, considering that positive and negative phases are approximated by the Friedlander equation and cubic polynomial, respectively. The plate is modeled as a SDOF system using von Karman's theory of large displacements. The development of the nonlinear dynamic differential equation is reviewed, considering a simply supported plate, and its solution is based on fourth order Runge-Kutta numerical method. A reference example is used as a benchmark and then parametric studies are conducted, in which the influence of scaled distance, mass of explosive, and the consideration or not of the negative phase is analyzed.

Keywords:
Plates; Dynamic Analysis; Nonlinearity; von Karman; Blast Load; Runge-Kutta

Graphical abstract

Graphical abstract

1 INTRODUCTION

Evaluation and search for criteria concerning the blast phenomenon have been studied recently after numerous events arising from both intentional and accidental explosions took place in various locations. In case structures are not designed considering this possibility, collapse may occur, generating great losses of buildings and lives (Kinney and Graham, 1985).

The first studies on blast loads were carried out by Friedlander (1940), showing that explosions create a rise in pressure above the atmosphere, also known as the positive phase, until it returns to ambient pressure after a certain interval of time. Thus, the Friedlander Equation was formulated to represent the behavior pressure vs. time in a shock wave. By then, only the positive phase of the blast load was considered.

Granström (1956) observed in that explosions can occur in two forms: spherical and hemispheric. The difference is the location of the explosive source, which can be triggered while still suspended in the air (spherical) or resting on a surface (hemispherical). In addition, a breakthrough in studies occurred when it was found that suction occurred after the end of the positive phase. The author used a cubic polynomial approximation to represent this negative phase.

Several studies on the negative phase of the blast phenomenon were carried on, such as Gantes and Pnevmatikos (2004), Dharani and Wei (2004) with their contributions to the Friedlander Extended Method; Teich & Gebbeken (2010) also presented valuable insights concerning the modelling of the negative phase. Rigby (2014) compared different equations presented in the literature for the case of the negative phase and based on experimental analyses, concluded that the cubic approximation is the one that best describes the experimental results.

Further studies by the US Department of Defense (2008) generalized the computation of load parameters by using several curves which depend on the scaled distance and mass of the explosive source. Later, Rigby (2014) presented equations for those curves, thus enabling their use in structural analyses.

In the case of dynamic loads acting on structures, there is the need for the consideration of large displacements which, in addition, demand the inclusion of membrane effects in the case of plates, for example. These effects consider in-plane behavior which is normally disregarded in the analyses of plates and one of the possible approaches are the well-known Von Karman (1910) equations. When energy-based methods like Galerkin are applied considering a Single Degree of Freedom (SDOF) system, the resulting equation is highly nonlinear, presenting a cubic term. Yamaki (1961) presents the deduction of the dynamic equations for thin plates, considering three types of membrane boundary conditions (immovable, movable and stress-free). The solution was based on Fourier series and was intended for dynamic loads in general, not focusing on explosions.

Houlston et al. (1985) concentrated on the prediction and measurement of the structural response of ship panels to free field air-blast explosions. Experimental and numerical results were presented, considering steel plates and full-scale stiffened panels. Chandrasekharappa and Srirangarajan (1987) presented a study about elastic plates submitted the positive phase of blast load. Gupta et al. (1987) presented a study considering a single degree of freedom (SDOF) simply supported plate, using the Friedlander equation to model the blast load. An analytical model was developed which considered the Lagrangian approach. The behavior of plates subjected to blast load was presented by Dharaneedpathy and Sudhesh (1990), in which three types of stiffener configurations were analyzed in a linear dynamics study, including damping. Krauthammer and Altenberg (2000) focused on the assessment of blast waves negative phase effects on glass panels. Librescu et al. (2004) studied the dynamic response of sandwich panels subjected to blast load, considering geometric non-linearity, initial imperfections, material anisotropy and the contribution of shear strains. Two types of blast loading were studied: underwater and air explosion, both considering spherical explosion. Feldgun et al. (2016) used the concept of large displacements in plates with membrane effect and explosive charge exclusively for the positive phase, to evaluate the behavior of plates in terms of displacements and stresses at its center. Kang et al. (2016) presented a simulation using blast loads, considering positive and negative phases.

This work focuses on the analysis of plates subjected to blast loads (positive and negative phases) considering membrane effects. The development of the nonlinear dynamic differential equation is reviewed, considering a simply supported plate. Equations were also developed for the obtention of blast load parameters, aiming at using closed-form expressions. The solution of the non-homogeneous differential equation is based on the numerical method of Runge-Kutta. A reference example of a square plate is used as a benchmark and then a parametric study is conducted, in which the influence of scaled distance, mass of explosive and the consideration or not of the negative phase is analyzed.

2 BLAST LOADS

2.1 Definition and load function

According to Granström (1956), an air blast consists of a mass of compressed air moving in the direction of propagation of the wave. It is a rapid release of energy that modifies the pressure of the environment and is dissipated into shock waves. When either an incident or a reflected shock wave reaches a target point, pressure instantly increases to a peak overpressure value. During the short period of time of this positive phase, pressure decays towards zero, after which a negative phase takes place. Then, the suction pressure increases up to a peak underpressure value and it is subsequently restored to zero at the end of the negative phase (Needham, 2010; Khaledy et al., 2018).

Figure 1 illustrates a representative blast pressure load p(t), where p_max is the peak overpressure, p_min is the peak underpressure, t_d is the duration of the positive phase, and t⁻_d is the duration of the negative phase. The positive and negative impulse values (i_d and i⁻_d, respectively) are related to the area under the pressure curve for each phase. All these parameters characterize a given blast load, and it will be addressed later.

Figure 1
Blast load at a target point due to a shock wave (Adapted from Rigby et al., 2013)

In literature, several studies have been carried out to formulate blast load equations based on empirical and semi-empirical methods. This work adopts the well-established exponential equation for the positive phase, proposed by Friedlander (1940), and a cubic polynomial equation for the negative phase, suggested by Granström (1956). The resulting blast load expression is

$p\left(t\right)=\left\{\begin{array}{c}{p}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(1-\frac{t}{t_d}\right)e^{-\frac{a^{\mathrm{\text{'}}}t}{t_d}},t\le t_d\\ -p_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(\frac{6.75\left(t-t_d\right)}{t_d^{-}}\right){\left(1-\frac{\left(t-t_d\right)}{t_d^{-}}\right)}^2,t_d\le t\le t_d+t_d^{-}\\ 0,t\ge t_d+t_d^{-}\end{array}\right.$ (1)

where a' is the decay coefficient of the positive phase, which can be obtained by enforcing the area under de curve to particular values of the positive impulse. It should be noticed that the evaluation of a' involves solving a nonlinear equation, which was accomplished by employing the Newton-Raphson method in the present study. The use of a numerical method is mandatory in this context because the blast wave behavior is approximated by imposing a continuous derivative (between positive and negative phases).

In contrast, the cubic equation is defined by only two parameters out of three (p_min, t⁻_d and i⁻_d) and it implies that the following relation holds

$$\frac{i_d^{-}}{p_{\min }{t}_d^{-}}=\frac{9}{16}$$ (2)

Also, the peak underpressure is placed at the first one third of the duration of the negative phase, and the first derivative at t= t_d is not imposed. Nevertheless, Equation (1) leads to an approximately smooth function when calibrated parameters from real blast loads are employed.

Most of the investigations on the response of structures under blast loads have applied only the positive phase of the load with Friedlander approximation. In these cases, the decay coefficient may be related to the potential damage caused to the structure, since it measures the positive impulse, which may be related to the amount of energy released during the explosion.

Nevertheless, it was just recently that the effect of the negative phase on the structure response was highlighted. One may find several approaches to approximate this portion of the pressure curve, namely linear (Krauthammer and Altenberg, 2000), extended Friedlander (Gantes and Pnevmatikos, 2004; Wei and Dharani, 2006 ), Teich extended Friedlander (Teich and Gebbeken, 2010), and cubic (Granström, 1956). However, Rigby (2014) found that the cubic polynomial function best approximates experimental blast pressure measurements, and so it was adopted here. Moreover, the use of Friedlander equation together with a cubic polynomial curve is recommended for non-confined explosions by the US Blast Resistant Structures Design Manual (NFEC, 1986).

2.2 Blast wave characterization

The values of the six parameters presented in Figure 1 depend on several factors: type of explosion (spherical or hemispherical explosion), type and angle of incidence of the blast wave (incident or reflected), type and amount of explosive material, and distance between explosive source and target. The present work focuses on reflected blast waves from both spherical (air) and hemispherical (surface) explosions, with a normal angle of incidence to the plane of the plate. For convenience, all parameters are usually expressed as a function of the scaled distance:

$$Z=\frac{R}{W_{TNT}^{1/3}}$$,(3)

where R is the distance between explosive source and target, and W_TNT is the equivalent mass of TNT.

Figure 2 shows the curves provided in the blast design manual “Structures to Resist the Effects of Accidental Explosions” by the US Department of Defense (2008) for spherical blasts, also reproduced by Rigby (2014) in SI units. Their corresponding expressions could be obtained by fitting curves to digitized data points for the parameters of the negative phase of hemispherical blasts. Equations (4) to (8) present expressions used in this work for spherical blasts. Equations (7) and (8) were developed by Rigby (2014) and Equations (4) to (6) are presented in Reis et al. (2019). Expressions for hemispherical blasts were also developed (Reis, 2019), but they will be omitted here since they are not used in this work. As the relation given by Eq. (2) is implied for the cubic approximation of the negative phase, the curves and expressions regarding the parameter t⁻_d could be suppressed. In addition, it is worth mentioning that the three parameters of the negative phase obey Eq. (2) only for Z>5, as observed by Granström (1956).

Figure 2
Parameters for spherical explosive charge (Adapted from RIGBY, 2014)

$$p_{\max }(Z)=\left\{\begin{array}{l}26303Z^{-1.2}\mathrm{,0.07}<Z\le 0.189\\ 13015Z^{-1.6}\mathrm{,0.189}<Z\le 0.487\\ 8100Z^{-2.3}\mathrm{,0.487}<Z\le 1.09\\ 8095Z^{-2.9}\mathrm{,1.09}<Z\le 2.41\\ 4941Z^{-2.4}\mathrm{,2.41}<Z\le 5.44\\ 864Z^{-1.4}\mathrm{,5.44}<Z\end{array}\right.\text{(kPa)}$$ (4)

$$\frac{t_d(Z)}{W_{TNT}^{1/3}}=\left\{\begin{array}{l}2.2Z^3-0.3Z^2-0.5Z+0.3\mathrm{,0.209}<Z\le 1.03\\ 8.2Z^3-34.8Z^2+48.9Z-20.4\mathrm{,1.03}<Z\le 1.49\\ -0.2Z^3+2.3Z^2-6.1Z+6.8\mathrm{,1.49}<Z\le 4.0\\ 1.5\log (Z)+1.2,Z>4.0\end{array}\right.\left(\frac{\text{ms}}{{\text{kg}}^{\text{1/3}}}\right)$$ (5)

$$\frac{i_d(Z)}{W_{TNT}^{1/3}}=\left\{\begin{array}{l}435Z^{-1.9}\mathrm{,0.071}<Z\le 0.198\\ 870Z^{-1.5}\mathrm{,0.198}<Z\le 1.28\\ 780Z^{-1.1},Z>1.28\end{array}\right.\left(\text{kPa}\frac{\text{ms}}{{\text{kg}}^{\text{1/3}}}\right)$$ (6)

$$p_{\min }(Z)=\left\{\begin{array}{l}101\mathrm{,0.071}<Z\le 0.668\\ -32.9Z^2+13Z+106\mathrm{,0.668}<Z\le 1.27\\ 93Z^{-1.22}\mathrm{,1.27}<Z\le 2.78\\ 73Z^{-0.978}\mathrm{,2.78}<Z\le 37.6\end{array}\right.\text{(kPa)}$$ (7)

$$\frac{i_d^{-}(Z)}{W_{TNT}^{1/3}}=\left\{\begin{array}{l}-724Z^2+445Z+553\mathrm{,0.071}<Z\le 0.580\\ 11.4Z^2-315Z+752\mathrm{,0.580}<Z\le 1.19\\ 462Z^{-0.88}\mathrm{,1.19}<Z\le 5.25\\ 434Z^{-0.848}\mathrm{,5.25}<Z\le 37.6\end{array}\right.\left(\text{kPa}\frac{\text{ms}}{{\text{kg}}^{\text{1/3}}}\right)$$ (8)

3 DYNAMIC ANALISYS OF PLATES

3.1 SDOF Model

Because of its simplicity and computational effortlessness, the Single Degree of Freedom (SDOF) method is a well-established approach in the dynamic analysis of two-dimensional structures (Morison, 2006). Based on some approximations, the SDOF model can be used to obtain the typical response of a structure, offering reliable results with much less computational effort than a full Finite Element analysis would require, especially when blast loads are considered. The SDOF is a simplified model that describes the structural behavior based on the response of a certain point of the structure, typically its midpoint, considering its dynamical properties (Stolz et al., 2014). Thus, observing that the SDOF method is a reliable and inexpensive mathematical procedure, the study of the nonlinear dynamic plate equation considering a single degree of freedom is the approach of this work. The US Department of Defense (2008) recommends the use of SDOF models for structures with simple geometry and elastic or perfect plastic behavior of its materials.

Even though the explosion phenomenon produces high frequency content, most of the blast’s energy excites the fundamental mode of the structure. The fundamental period is usually one or two times longer than the positive phase duration. Moreover, the American Society of Civil Engineers (ASCE, 2021) states that, in the general case, SDOF models can be used when the explosive source is far from the bulkhead and the wave pressure is uniformly distributed in the surface of the structure. All the examples analyzed in the paper are modeled with such conditions and the simplicity of SDOF method allows the extension of the analysis to other combinations of scaled distances and masses of explosive, thus enabling a comprehensive parametric study.

For the following analyses, a rectangular plate is considered, with dimensions a and b, parallel to the x and y axes, respectively, as shown in Fig. 3. The thickness of the plate has dimension h, much smaller than a and b, which allows the Kirchhoff simplification, therefore neglecting shear strains. Moreover, all the equations are written considering the plate’s average plane, located at z = 0. The loading p_z is always normal to the surface of the plate and is considered as uniformly distributed.

Figure 3
Plate dimensions and loading feature (Adapted from Feldgun et al., 2016)

Since explosive loading causes high impact on the structure, it should be considered that large displacements can be generated, requiring the inclusion of higher order terms in the strain tensor. In this work, von Karman's equations are used, according to:

$$L\left(u_z,\varphi \right)\Rightarrow D{\nabla }^4{u}_z+\rho h\frac{{\partial }^2{u}_z}{\partial t^2}-h\left(\frac{{\partial }^2\varphi }{\partial y^2}\frac{{\partial }^2{u}_z}{\partial x^2}+\frac{{\partial }^2\varphi }{\partial x^2}\frac{{\partial }^2{u}_z}{\partial y^2}-2\frac{{\partial }^2\varphi }{\partial x\partial y}\frac{{\partial }^2{u}_z}{\partial x\partial y}\right)=P_z\left(x,y,t\right)$$,(9)

$${\nabla }^4\varphi =E\left[{\left(\frac{{\partial }^2{u}_z}{\partial x\partial y}\right)}^2-\frac{{\partial }^2{u}_z}{\partial x^2}\frac{{\partial }^2{u}_z}{\partial y^2}\right]$$,(10)

where L is the differential operator, D is the plate’s bending stiffness, u_z is the displacement in the z direction, ρ is the material specific mass, h is the plate thickness, $\varphi$ is the Airy function and E is the Young modulus.

Considering that the plate is simply supported a Navier series is used as an approximation for displacement u_z:

$$u_z\left(x,y,t\right)=A\left(t\right)\cos \left(\frac{\pi x}{a}\right)\cos \left(\frac{\pi y}{b}\right)$$,(11)

where A(t) corresponds to the amplitude that arises from the separation of time and space variables.

3.2 Membrane Boundary conditions

Substituting Eq. (11) in Eq. (10) the following is obtained:

$${\nabla }^4\varphi =-\frac{E{\left[A\left(t\right)\right]}^2{\pi }^4{h}^2}{2{\left(ab\right)}^2}\left[\cos \left(\frac{2\pi x}{a}\right)+\cos \left(\frac{2\pi y}{b}\right)\right]$$ (12)

Based on Yamaki (1961), the function ϕ can be generalized as a series, considering the sum of a homogeneous and a particular solution, which can be rewritten as:

$$\varphi ={\left[\frac{A(t)}{h}\right]}^2\left(\frac{1}{2}{p}_x^{\text{'}}{y}^2+\frac{1}{2}{p}_y^{\text{'}}{x}^2+E{h}^2{\displaystyle \sum _{p=0}^{\infty }{\displaystyle \sum _{q=0}^{\infty }{\Phi }_{pq}{X}_{2p}{Y}_{2q}}}\right)$$ (13)

Where:

$${\Phi }_{p,q}={\phi }_{p,q}+\phi {\text{'}}_{p,q}$$

$${\phi }_{p,q}=\left\{p\wedge q\in \mathbb{N}/{\phi }_{\mathrm{0,1}}=-\frac{1}{32{\beta }^2},{\phi }_{\mathrm{1,0}}=-\frac{{\beta }^2}{32}\right\}$$ (14)

$$\begin{array}{l}{\phi }_{p,q}^{\text{'}}=\frac{4\beta }{\pi {\left(p^2+{\left(q\beta \right)}^2\right)}^2}\left[\frac{p{\left(-1\right)}^q{\epsilon }_q{\sinh }^2\left(\frac{p}{\beta }\pi \right)A_p}{\sinh \left(\frac{p}{\beta }\pi \right)\cosh \left(\frac{p}{\beta }\pi \right)+\frac{p}{\beta }\pi }+\frac{q{\left(-1\right)}^p{\epsilon }_p{\sinh }^2\left(q\beta \pi \right)B_q}{\sinh \left(q\beta \pi \right)\cosh \left(q\beta \pi \right)+q\beta \pi }\right]\\ p,q=\left\{p\wedge q\in \mathbb{N}/\left(p\vee q=0\Rightarrow {\epsilon }_p\vee {\epsilon }_q=1/2\right)\wedge \left(p\vee q\ne 0\Rightarrow {\epsilon }_p\vee {\epsilon }_q=1\right)\right\}.\end{array}$$ (15)

$$X_{2p}=\cos \left(\frac{2\pi px}{a}\right),Y_{2q}=\cos \left(\frac{2\pi qy}{b}\right)$$ (16)

In Eq. (13), p’_x and p’_y are the nondimensional coefficients of the homogeneous solution of the Airy function, φ_p,q and φ’_p,q are the coefficients of the Fourier series for the and homogeneous solution of the Airy function, respectively. At this point, it is important to mention that there are three boundary conditions corresponding to the membrane effect, which affect the Airy function: stress free, immovable, and movable. The physical meaning of these conditions in terms of displacements, forces and/or stresses are given in Table 1 (Yamaki, 1961). In this table, P_i is the resultant of the forces acting on the i-axis.

Table 1
Specifications of membrane stress conditions

Case	$x=\pm a/2$	$y=\pm b/2$
Stress Free	${\varphi }_{,yy}={\varphi }_{,xy}=0$	${\varphi }_{,xx}={\varphi }_{,xy}=0$
Immovable	$u_x={\varphi }_{,xy}=0$	$u_y={\varphi }_{,xy}=0$
Movable	$P_x={\varphi }_{,xy}=0$	$P_y={\varphi }_{,xy}=0$

The general solution presented in Equation (13) is interesting, since it simplifies the imposition of membrane boundary conditions (stress free, immovable, and movable), by only requiring the calculation of parameters A_p, B_q, p’_x and p’_y.

Applying the conditions in Table 1 based on Airy's general expression (Equation (13)) , p’_x and p’_y expressions are determined for each support condition as in Table 2.

Table 2
Specifications of membrane stress conditions

Case	$p_x^{\text{'}}$	$p_y^{\text{'}}$
Stress Free	0	0
Immovable	$\frac{{\pi }^{2}E{h}^2\left({\beta }^2\nu +1\right)}{8a^2\left(1-{\nu }^2\right)}$	$\frac{{\pi }^{2}E{h}^2\left({\beta }^2+\nu \right)}{8a^2\left(1-{\nu }^2\right)}$
Movable	0	0

For immovable and movable conditions, parameters A and B are null. For the stress-free case, A_p and B_q can be obtained by solving the following linear system of equations:

$$\begin{array}{l}\left\{\begin{array}{l}{B}_n+{\displaystyle \sum _{m=1}^{\infty }{\left(-1\right)}^m{A}_m\eta \left(\frac{m}{\beta },n\right)}=-{\beta }^2{n}^2{\displaystyle \sum _{p=0}^{\infty }{\left(-1\right)}^p{\phi }_{p,n}}\\ A_n+{\displaystyle \sum _{m=1}^{\infty }{\left(-1\right)}^m{B}_m\eta \left(m\beta ,n\right)}=-n^2{\displaystyle \sum _{q=0}^1{\left(-1\right)}^q{\phi }_{n,q}}\end{array}\right.\\ \eta \left(\gamma ,\xi \right)=\frac{4{\left(-1\right)}^m\gamma {\xi }^2}{\pi {\left({\gamma }^2+{\xi }^2\right)}^2}\frac{{\sinh }^2\left(\gamma \pi \right)}{\sinh \left(\gamma \pi \right)\cosh \left(\gamma \pi \right)+\gamma \pi }.\end{array}$$ (17)

3.3. Nonlinear Dynamic differential equation

Applying the Galerkin Method, the problem can be written in variational form as the minimization of the potential energy:

$${\displaystyle \underset{-a/2}{\overset{a/2}{\int }}{\displaystyle \underset{-b/2}{\overset{b/2}{\int }}L\left(u_z,\varphi \right)u_{z}dxdy=0}}$$ (18)

Another simplification is the consideration of the load as uniformly distributed, thus disregarding its dependence on spatial variables x and y. From now on, p(x,y,t) will be simply written as p(t). After solving Equation (18) using the approximation for u_z and ϕ given in previous sections, the following equation is obtained:

$$\frac{{\partial }^{2}A(t)}{\partial t^2}+\frac{{\pi }^{4}D{\left({\beta }^2+1\right)}^2}{a^{4}h\rho }A(t)+\left[\frac{{\pi }^2}{a^2\rho }\left(p_x^{\text{'}}+{\beta }^2{p}_y^{\text{'}}\right)-\frac{2{\pi }^{4}E{h}^2{\beta }^2}{a^4\rho }\left({\Phi }_{\mathrm{0,1}}+{\Phi }_{\mathrm{1,0}}\right)\right]{\left[A(t)\right]}^3=\frac{16}{{\pi }^2{h}^2\rho }p(t)$$,(19)

where β is the relation between the largest and smallest dimensions of the plate (considered in this work as x and y dimensions, respectively), i.e., β = a/b.

According to Feldgun et al. (2016), Eq. (19) can be looked at as a SDOF and the coefficients related to the linear and third-degree terms can be summarized, obtaining a much simpler equation:

$$\frac{{\text{d}}^{2}A(t)}{\text{d}{t}^2}+K_{1}A(t)+K_3{\left[A(t)\right]}^3=f\left(t\right)$$ (20)

It becomes evident that the dynamic differential equation of the plate is directly dependent on a linear portion K₁, which corresponds to the plate (bending) boundary conditions (simply supported in this work), and a nonlinear portion K₃ related to the membrane conditions (immovable, movable and stress-free). Equation (20) is a general formulation that is valid for every combination of boundary conditions (plate and membrane). For a simply supported plate, parameters K₁ and K₃ are given in Table 3.

Table 3
Parameters K₁ and K₃

Case	$K_3$	$K_1$
Stress Free	$\frac{2{\pi }^{4}E{h}^2{\beta }^2\left({\Phi }_{\mathrm{0,1}}+{\Phi }_{\mathrm{1,0}}\right)}{a^6\rho }$	$\frac{{\pi }^{4}D{\left({\beta }^2+1\right)}^2}{a^{4}h\rho }$
Immovable	$\frac{{\pi }^{4}E{h}^2\left[\left({\beta }^4+1\right)\left({\nu }^2-3\right){\nu }^2-4{\beta }^2\nu \right]}{16a^4\rho \left(1-{\nu }^2\right)}$
Movable	$\frac{{\pi }^{4}E{h}^2\left({\beta }^4+1\right)}{16a^4\rho }$

Thus, based on the data presented in Table 3, a nondimensional analysis is performed comparing the bending and membrane coefficients (K₁ and K₃, respectively) with respect to β, as shown in Figure 4. It is verified that the immovable system presents a greater influence of the nonlinear regime in comparison to the bending behavior of the structure. This is justified by the fact that the immovable system is characterized by having greater restrictions imposed on the edges of the plate. Thus, the nonlinearity becomes more evident for this condition. The stress-free case is characterized by presenting the opposite behavior of the immovable system, since there is greater freedom for lateral displacement, and less influence of the membrane effect in the system and greater influence of bending stiffness. The movable case is characterized by an intermediate behavior between immovable and stress-free conditions.

Figure 4
Comparison between parameters K₁ and K₃, case simply supported

3.4. Solution Strategy for the Nonlinear Dynamic Differential Equation

Noting that Equation (20) presents a high degree of nonlinearity, analytical methods are not able to directly solve such differential equation. The application of numerical methods becomes necessary and the fourth order Runge-Kutta numerical method was chosen. Since the Runge-Kutta method requires the differential equations to be of first order, Equation (20) should be rewritten as a system of two first order equations, as presented in Equation (21). An auxiliary function H(t) is used for this purpose.

$$\begin{array}{l}\frac{dA(t)}{dt}=H(t)\\ \frac{dH(t)}{dt}=F(t)-K_{1}A(t)-K_3{\left[A(t)\right]}^3\end{array}$$ (21)

It is important to mention a special case of Equation (20), when the external force f(t) is null. This case is known as the free vibration behavior, and it is useful to obtain the natural periods (or frequencies) and modes of the structure. These periods are directly influenced by the bending and membrane behavior, and, for that reason, they are referred to as linear (bending) and nonlinear (bending plus membrane) natural periods, respectively. Thus, Eq. (21) becomes a non-linear homogeneous equation known as the Duffing equation, a structural model that includes a non-linear force restoration which is proportional to the cubic power of the displacement (Rand, 2005). This equation has an analytical solution based on Jacob's concept of elliptical functions, as demonstrated by Soudack (1964), according to equations (22) to (24), in which Cn(x,y) is the elliptic cosine and A(t) is the solution of the free vibration problem.

$$A(t)=A_{0}Cn\left(\lambda t,k\right)$$ (22)

$$k=\sqrt{\frac{K_3{A}_0^2}{2\left(K_1+K_3{A}_0^2\right)}}$$ (23)

$$\lambda =\sqrt{K_1+K_3{A}_0^2}$$ (24)

Based on equations (22) to (24), it is possible to determine the expressions referring to the non-linear and linear natural periods:

$$T_{NL}=\frac{4{\displaystyle \underset{0}{\overset{\pi /2}{\int }}\frac{dx}{\sqrt{1-k^2{\sin }^{2}x}}}}{\lambda }$$ (25)

$$T_L=\frac{2\pi }{\sqrt{K_1}}$$ (26)

Finally, another important analysis to be considered is the Dynamic Amplification Factor (DAF), which is the ratio between the maximum dynamic displacement (u_z,max) and the static displacement (u_z,st), obtained by using the peak overpressure as the static load.

$$DAF=\frac{u_{z,\max }}{u_{z,st}}$$ (27)

4 NUMERICAL EXAMPLES

Two examples were formulated based on results found in the literature. The first one is a thin steel plate studied by Houlston et al. (1985). In this research, the results of the experimental data are used as loading parameters. Experimental results are compared with numerical ones (Feldgun et al., 2016) and then a parametric study of the plate is presented.

The second example is based on the one presented by Krauthammer and Altenberg (2000). The structure is a laminated glass plate subjected to a blast load, in which the TNT mass and the scaled distance Z (m/kg^1/3) are known. However, the loads considered are linear and bilinear approximations for the positive and negative phases, respectively.

4.1 Houlston et al. (1985)

The first example is a study based on the experimental data obtained by Houlston et al. (1985), also used as a numerical comparison by Feldgun et al. (2016). Geometric and material properties of this example are given in Table 4. As the blast load parameters in the original work are results from the experimental tests, a solver was used to determine ideal values for Z (m/kg^1/3) and W (kg) that can be then fed to the numerical model. The parametric study was performed in a model developed in Matlab.

Table 4
Example 1 parameters (Houlston et al., 1985)

Parameters	Valor
Dimension (a x b)	0.508 m x 0.508 m
Thickness (h)	0.0034 m
Young’s modulus (E)	207 GPa
Poisson’s coefficient (ν)	0.3
Mass density (ρ)	7770 kg/m3

The main parameters of experimental data resulting from the blast load are presented according to Table 5. Since the distance between the plate and TNT's mass used in the experiment is unknown, a solver was used to determine the ideal value of Z (m/kg^1/3) and W (kg) based on the results of the experimental data. The least squares method was used in the resulting values were Z = 5.64 m/kg^1/3 and W = 0.24 kg.

Table 5
Experimental blast parameters (Houlston et al., 1985)

Parameters	Valor
Pmax	57086.0019 Pa
Pmin	15420.247 Pa
td	0.002026061 s
id	37.3129 Pa.s
im	40.536 Pa.s
tm	0.004673336 s

Figure 5 compares the curves of the blast load with experimental data. Curve [1] represents the results obtained directly from the experimental test. The second curve [2] is obtained by the direct application of the blast parameters (see Table 5), according to Equations (1) and (2). Finally, the last curve [3] is obtained by using values for Z (m/kg^1/3) and W (kg) given by the least squares solver directly in the equations developed by Rigby (2014) and the abacus of the US Department of Defense (2008).

Figure 5
Comparison between experimental data (Houlston et al., 1985), analytical data and solver data

Figure 6 shows the curves plotted for comparison between the experimental test and numerical analyses. Numerical solutions are based in Equation (20) which was solved using a Runge-Kutta scheme as shown in section 3.4. Curve [1] is based on experimental data obtained by Houlston et al. (1985), curves [2] to [4] are obtained considering simply supported boundary condition. It is important to emphasize that in the models that consider the negative phase it is approximated by cubic polynomial. It is evident that when membrane effect is neglected (curve [4]) results are not in good agreement with experimental ones. For comparison purposes, curves [5] and [6] are the ones obtained by Feldgun et al. (2016) considering or not the negative phase, respectively. Finally, curve [7] is obtained by applying a load whose parameters are calculated using optimization techniques (least squares solver). According to Figure 6 and the studies by Feldgun et al. (2016), the simply supported boundary condition is the one that best represents the behavior of the structure used by Houlston et al. (1985). This case presents natural frequency and maximum displacements closest to the experimental results. For the experimental test, the largest displacement is 5.839 mm. In the case of the analysis with simply supported plate, u_z,max = 5.571 mm, an error of 4.58%; for the clamped boundary condition, u_z,max = 4.63 mm, an error of 20.70%. The importance of the negative phase consideration can also be assessed in Figure 6, which shows that its non-consideration affects the plate’s behavior in terms of displacement and frequency.

Figure 6
Midpoint displacement using experimental data, analytical and solver solution for Houlston et al. (1985)

The clamped (CC) boundary condition is assessed in Figure 7 as the plate was originally modeled using this boundary condition in the experiment by Houlston et al. (1985). It is possible to verify that results for simply supported boundary condition are better than the clamped case, which was also shown by Feldgun et al. (2016). All the expressions needed to implement the CC boundary conditions are presented in Reis (2019).

Figure 7
– Comparison between SS and CC conditions with negative phase considered

4.1.1. Parametric Study

After the validation of the model in previous sections, a parametric study was conducted to check the influence of some variable changes in the overall behavior of the plate. The first study used different blast loads, obtained from constant Z (m/kg^1/3) and different TNT masses. Using the parameters given in Houlston et al. (1985) the solver solution was found as Z = 5.64 m/kg^1/3 and this value was kept constant. For this purpose, the ratio between the linear and nonlinear maximum displacement amplitudes (u_L/u_z) was plotted for each membrane boundary condition as shown in Figure 8. In addition, all results were compared to the case of a blast load pressure without the negative phase.

Figure 8
Ratio uL/uZ versus WTNT for various values of Z when is used the Solver data

As expected, the plate with immovable boundary condition presented larger values of u_L/u_z, suggesting that the membrane effect is dominant. Not only the membrane effect is essential, but it is also possible to observe a peak behavior located between W = 0 kg and W = 25 kg in all curves in Figure 8, which is characteristic of the negative phase in the analysis of the blast load. As will be shown later, the inclusion of the negative phase increases the amount of energy transferred to the system, changing the natural frequency of the plate, and possibly leading to resonance.

The inclusion of the negative phase into the blast pressure load showed to be relevant for lower values of scaled distance and equivalent mass of TNT. The study was extended by using different values of scaled distance, resulting in Figure 9. Only the immovable case was investigated since it is the worst case in terms of displacements. A limiting curve beyond which the negative phase may be disregarded could be established. It is worth mentioning that this limiting curve was built by evaluating displacement results for several values of scaled distances, even though only three curves were plotted for the sake of clarity. It should also be noted that the negative phase of the blast load is dismissed beyond the limiting curve for the evaluation of maximum displacement amplitudes, but not necessarily for frequency values.

Figure 9
Ratio uL/uZ versus WTNT for various values of Z when is used the Solver data

Figure 10 is obtained with the variation of Z (m/kg^1/3) and its influence in the relative displacement u_z/h. TNT mass calculated by the solver based on the experimental results from Houlston et al. (1985) was W = 0.24 kg. It was used as a start value and compared with various other masses to analyze its influence in the relative displacement of the structure. It is verified, naturally, that larger masses cause larger displacements. As the scaled distance Z increases, the displacement decreases, which is also a natural conclusion since the energy of the blast wave tends to decrease with Z.

Figure 10
Ratio uz/h versus Z (m/kg^1/3) varying values of W_TNT for the Solver data for Houlston et al. (1985)

As shown in Figure 10 and 11, there is no evident limit for Z values as to disregard the negative phase. In this way, it is not possible to obtain a limit value for Z below which the negative phase can be neglected. This fact is also corroborated by Granström (1956), which shows that a value of Z = 5 m/kg^1/3 is a lower bound when a cubic approximation is used for the negative phase. In other words, for values below Z = 5 m/kg^1/3 the negative phase is insignificant when compared to the impulse provided by the positive phase.

Figure 11
Ratio uz/h versus Z (m/kg^1/3) comparing analyses with and without negative phase

Considering the structure studied by Houlston et al. (1985), a parametric analysis can be performed in terms of the Dynamic Amplification Factor (DAF), using four case combinations, as shown in Figure 12: the consideration or not of the negative phase, as well as the activation or not of the membrane effect. Curve [1] represents the real conditions of the test, that is, nonlinearity and presence of the negative phase. In curves [1] and [3], which are those that consider the negative phase, it is possible to verify the same "peak" behavior presented in Figure 8 It is worth mentioning that in these cases the nonlinear period may have an important role in the resonance behavior.

Figure 12
DAF versus td / TL for different values of W

Curves [2] and [4] are the ones obtained without the application of the negative phase in blast load. As mentioned before, the contribution of the negative phase is evident in this type of structure. The dots in Figure 12 represent the points in each curve using the parameters obtained from the experimental test. It is observed that the actual test conditions (nonlinear analysis and consideration of negative phase) generate a ratio of t_d/TL which falls within the region where “peaks” occur. In this way, it is worth noting that performing a linear analysis in the structure or the non-consideration of the negative phase in the loading can lead to the errors up to 73%.

The existence of peaks in the DAF behavior suggests the occurrence of resonance, due to proximity between the natural periods of the load and the plate. This hypothesis can be tested by performing a Fast Fourier Transform (FFT) in the load data and comparing the results with the natural periods of the plate. There are two aspects to consider in this analysis: 1) since the load is non-periodic its FFT will not generate clearly defined natural frequencies/periods; 2) the natural periods of the plate should be the nonlinear ones, since the impulsive characteristics of the blast loading will alter the fundamental periods, as shown previously. Figure 13 shows this analysis considering (Figure 13-a) or not (Figure 13-b) the negative phase. There is also the possibility of a great influence of the loading characteristics and the nonlinear period which can also affect this kind of analysis. This is still a work in progress and further studies should be performed to support that statement.

Figure 13
FFT analysis of the load, immovable condition, (a) with negative phase; (b) without negative phase

4.2 Krauthammer and Altenberg (2000)

Krauthammer and Altenberg (2000) presented a numerical parametric study of laminated glass plates, considered as homogeneous structures. As the plate is made of a fragile material, different from steel, an elastic linear analysis is performed, with simple support as boundary conditions and without damping. The blast load also includes the negative phase using a bilinear approximation. Table 6 shows the physical characteristics of the plate. In their parametric analysis, Krauthammer and Altenberg (2000) studied the behavior of four plates, two of which are addressed in this work. Both are square plates, but with different blast loads, that is, W = 10 kg and W = 100 kg. Data on the blast load are shown in Table 7.

Table 6
Example 2 parameters (Krauthammer and Altenberg, 2000)

Parameter	Value
Dimension (a x b)	1.397 m x 1.448 m
Thickness (h)	0.00963 m
Young’s modulus (E)	72 GPa
Poisson’s coefficient (ν)	0.25
Mass density (ρ)	2500 kg/m3

Table 7
Loading parameters (Krauthammer and Altenberg, 2000)

Parameters	Valor
Z	0 – 39.64 m / kg1/3
W	10 kg and 100 kg

The first parametric study presents the relationship between the ratio u_z/h and Z (m/kg^1/3). In this case, the reference result is compared with 4 types of analysis: the combination of considering or not the nonlinearity in geometry, together with the consideration or not of the negative phase. Krauthammer and Altenberg (2000) presented a study considering a probabilistic case of collapse of the structure. As expected, good results are obtained for larger values of scaled distances and small masses of explosive, conditions that guarantee that the structure will remain within linear material behavior. Nevertheless, results are acceptable even for smaller scaled distances, as shown in Figure 14. Moreover, this analysis also validates the consideration of the negative phase, especially for the case with W = 10kg.

Figure 14
uz/h analysis of the load with Solver data solution, immovable condition, (a) W = 10 kg; (b) W = 100 kg

By analyzing the structure subjected to a TNT equivalent mass of W = 10 kg, Figure 14-a, it is confirmed that the linear analysis with consideration of the negative phase, curve [4], is the one that best fits the reference.

For the case with W = 100 kg, Figure 14-b, it is verified that the results of Krauthammer and Altenberg (2000) are different from those given by the curve [4], especially as Z grows. This case is justified by the equation of the negative phase, which is represented by a bilinear approximation in the original paper, while this work uses the cubic formulation. This difference becomes more evident for bigger masses.

When analyzing the influence of the load time t_d in the Dynamic Amplification Factor (DAF), the consideration or not of nonlinearity and the negative phase plays an important role in the behavior of the plate, as shown in Figure 15.

Figure 15
DAF versus td / TL

In Figure 15, curves [1] and [2] represent the behavior of the structure when geometric nonlinearity is activated. Curves [3] and [4] represent linear analysis, with and without the negative loading phase, respectively. Curve [5] shows the results for Z = 7.928 m/kg^1/3 and W = 10 kg, and curve [6] is obtained for the same scaled distance and W = 100 kg. The value for the parameter Z (m/kg^1/3) was chosen since it is the lowest value among the set of scaled distances Z (m/kg^1/3) used by Krauthammer and Altenberg (2000). These combinations for Z and W give t_d/T_L = 0.2014 and t_d / T_L = 0.4341, as shown.

Finally, FFT was performed considering linear analysis (since it was the best fit) considering or not the negative phase (Figure 16-a and Figure 16-b, respectively). The FFT analysis, as previously discussed, is justified since "peaks" are generated in the DAF plot, suggesting resonance. Nevertheless, since a linear analysis is performed, there is no change in the plate’s natural period.

Figure 16
FFT of the load based on Linear Analysis (a) with Negative Phase (b) without Negative Phase

As stated before, the best agreement between the numerical model and experimental results occurs for low masses and high scaled distances. For those cases, the negative phase has the greatest influence on displacements and its consideration is paramount. Results presented in this work and the one by Krauthammer and Altenberg (2000) are close especially for those conditions that guarantee that the structure will remain within linear material behavior, which is true for glass. Besides, Krauthammer and Altenberg (2000) also explored the use of the negative phase in their study, although here its expression is obtained using a cubic polynomial equation instead of piecewise linear.

4 FINAL REMARKS

This work aimed to evaluate the nonlinear dynamic behavior of rectangular plates considering the membrane effect, when subjected to blast loads. Von Karman's plate theory for large displacements was used to obtain the differential equations, Galerkin Method was employed to determine stiffness parameters and Runge-Kutta's numerical routine was finally used to solve the nonlinear equations. The explosive load imposed on the structure was divided into positive and negative phases, represented by the Friedlander equation and the cubic polynomial equation, respectively.

Considering the application of the negative phase in the analysis of plates subjected to blast loads, this work corroborates its importance by showing good accordance of the formulation with results from literature. It is also worth noticing that the negative phase approximation by a cubic polynomial equation presents good results when compared with other formulations (piecewise linear for instance).

The membrane effect was considered in the plate behavior, and three types of lateral conditions were analyzed: immovable, movable and stress free. It was shown that for greater lateral restrictions (immovable is more restrictive than stress free) influence of the membrane effect on the time history is more evident; displacements in this case are smaller, generating bigger ratios between nonlinear and linear displacements. This conclusion is in accordance with previous studies and corroborates the need for the consideration of membrane effect.

Results from the literature, both numerical and experimental, were used to validate the formulation. After validation, parametric studies were performed. Considering the new load formulation, it is possible to optimize the parameters to obtain the best fit for experimental blast data. This approach made the parametric study possible since different values of explosive mass and scaled distance can be easily applied to the structure.

Based on the presented results, it was found that the consideration of the negative phase in the dynamic behavior of the plate is relevant, since in some cases the largest displacements in the time history are observed during the negative phase. There is also influence of the negative phase consideration especially for low scaled distances and large masses of explosive, in which the energy transferred to the system can induce changes in natural frequencies (periods).

The study of the Dynamic Amplification Factor for different ratios between load positive time and the plate’s natural period shows that there are some peaks which suggest resonance when it is considering the negative phase in analysis. The influence of the nonlinear period which arises from the membrane effect is still a subject for further studies.

Acknowledgements

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Finance Code 001 and FAPERJ − Grant Number E-26/010.002150/2019.

• https://doi.org/10.1590/1679-78256706

References

• ASCE, American Society of Civil Engineers. Structural Design for Physical Security. ASCE Manuals and Reports on Engineering Practice N° 142. Structural Engineering Institute, 2021.
• Chandrasekharappa, G., Srirangarajan, H.R. (1987). Nonlinear Response of Elastic Plates to Pulse Excitations. Computers & Structures 27(3):373-378.
• Dharaneedpathy, M.V., Sudhesh, K.G. (1990). Optimal Stiffening of Square Plates Subjected to Air-Blast Loading. Computers & Structures 36(5):891-899.
• Dharani, L.R., Wei, J. (2004). Dynamic response of laminated glass under blast loading: effect of negative phase. Proceedings of Structures under shock and impact VIII, 181-190, Crete, Greece.
• Feldgun, V.R., Yankelevsky, Y.S., Karinski, Y.S. (2016). A nonlinear SDOF model for blast response simulation of elastic thin rectangular plates. International Journal of Impact Engineering 88:172–188.
• Friedlander, F.G. (1940). The diffraction of sound pulses I. Diffraction by a semi-infinite plane. Communicated by G. I. Taylor, F.R.S., 1940.
• Gantes, C., Pnevmatikos, N. (2004). Elastic-plastic response spectra for exponential blast loading. International Journal of Impact Engineering 30(3):323-343.
• Granström, S.A. (1956). Loading characteristics of fair blasts from detonating charges. Technical Report 100, Transactions of the Royal Institute of Technology, Stockholm.
• Gupta, A.D., Gregory, F.H., Bitting, R.L., Bhattachary, S. (1987). Dynamic Analysis of an Explosive Loaded Hinged Rectangular Plate. Pergamon Journal Ltd, Computers & Structures 26339-344.
• Houlston, R., Slater, J. E., Pegg, N., Desrochers, C. G. (1985). On Analysis of Structural Response of Ship Panels Subjected to Air Blast Loading. Computers & Structures 21:273-289.
• Kang, K.Y, Choi, K.H, Choi, J.W, Ryu, Y.H., Lee, J.M. (2016). An Influence of Gas Explosions on Dynamic Responses of a Single Degree of Freedom Model. Shock and Vibration 2016, ID 9582702.
• Khaledy, N., Habibi, A., Memarzadeh, P.A. (2018). Comparison between different techniques for optimum design of steel frames subjected to blast. Latin American Journal of Solids and Structures 15(9):e106.
• Kinney, G.F., Graham, K.J. (1985). Explosive Shocks in Air. Springer Science & Business Media (New York).
• Krauthammer, T., Altenberg, A. (2000). Negative phase blast effects on glass panels. International Journal of Impact Engineering 24(1):1-17.
• Librescu, L., Oh, S. Y., Hohe, J. (2004). Linear and non-linear dynamic response of sandwich panels to blast loading. Composites: Part B 35:673-683.
• Morison, C.M. (2006). Dynamic response of walls and slabs by single-degree-freedom analysis – a critical review and revision. International Journal of Impact Engineering 32(8):1214-1247.
• Naval Facilities Engineering Command. (1986). Blast resistant structures. Alexandria, Va., USA NFE DM 2-08.
• Needham, C.E., (2010). Blast Waves, Shock Wave and High-Pressure Phenomena. Springer-Verlag (Berlin).
• Rand, R. H. (2005). Lecture Notes on Nonlinear Vibrations. Department Theoretical & Applied Mechanics, Cornell University (New York).
• Reis, A.W.Q.R. (2019). Dynamic analysis of plates subjected to blast load. M. Sc. Dissertation (in Portuguese), Rio de Janeiro State University, Brazil.
• Reis, A.W.Q.R., Burgos, R.B., Oliveira, M.F.F. (2019). Blast Wave Analysis in Plate Structures Considering Membrane Effect, Proceedings of the Ibero-Latin-American Congress on Computational Methods in Engineering (ABMEC), Natal/RN, Brazil.
• Rigby, S. E., Andrew, T., Bennett, T., Clarke, S. D., Fay, S. D. (2013). The Negative Phase of the Blast Load. International Journal of Protective Structures 5(1):1-19.
• Rigby, S.E. (2014). Blast Wave Clearing Effects on Finite - Sized Targets Subjected to Explosive Loads. Ph. D. Thesis. Department of Civil and Structural Engineering, University of Sheffield, UK.
• Soudack, A.C. (1964). Nonlinear Differential Equations Satisfied by the Jacobian Elliptic Functions. Mathematics Magazine 37(3):138-147.
• Stolz, A., Fischer, K., Roller, C., Hauser., S. (2014). Dynamic bearing capacity of ductile concrete plates under blast loading. International Journal of Impact Engineering 69:25-38.
• Teich, M., Gebbeken, N. (2010). The influence of the underpressure phase on the dynamic response of structures subjected to blast loads. International Journal of Protective Structures 1(2):219-233.
• US Department of Defense. (2008). Structures to resist the effects of accidental explosions, US DoD, Washington DC, USA, UFC-3-340-02.
• Von Karman, T. (1910). Festigkeitsprobleme im Maschinenbau. Encyklopädie der Mathematischen Wissenschaften IV:311-385.
• Wei, J. Dharani, L. R. (2006). Response of laminated architectural glazing subjected to blast loading. International Journal of Impact Engineering 32(12):2032-2047.
• Yamaki, N. (1961). Influence of large amplitudes on flexural vibrations of elastic plates. ZAMM 41(12):501-510.

Edited by

Editor: Rogério José Marczak

Publication Dates

• Publication in this collection
  05 Jan 2022
• Date of issue
  2022

History

• Received
  20 Aug 2021
• Accepted
  16 Nov 2021