Dynamic analysis of shell structures with application to blast resistant doors

This paper concerns the dynamic analysis of shell structures, with emphasis on application to steel and steel-concrete composite blast resistant doors. In view of the short duration and impulsive nature of the blast loading, an explicit integration method is adopted. This approach avoids time-consuming computations of structural stiffness matrix and solving of simultaneous nonlinear equations. Single-point quadrature shell elements are used, with numerical control to suppress spurious hourglass modes. Composite shells are handled by an appropriate integration rule across the thickness. Both material and geometric nonlinearities are accounted for in the formulation. Contact and gap problems are considered using bilinear spring elements in the finite element analysis. Numerical examples are presented for some benchmark problems and application study to blast resistant doors. Good correlation is generally obtained between the numerical results based on the software developed and the results obtained by other means including field blast tests.


Introduction
Blast resistant doors are commonly used in defence shelters, ammunition storage dumps and blast-resistant structures alike.These doors have to be designed to withstand blast loading and shock waves, so as to protect occupants and objects inside, and in the case of ammunition storage, to prevent chain reaction of explosions from one compartment to another.Inadequate design of the blast resistant door would undermine the operational performance and survivability of the structure and users.On the other hand, an overly conservative design would result in incurring unnecessarily high costs.
Recent studies on blast resistant structures mainly dealt with reinforced concrete structures, see [1][2][3][4][5] for example.Lok and Xiao [6] extended the study to steel-fibre-reinforced concrete panels subjected to blast loading.Krauthammer carried out numerical studies on the blast response of structural concrete and structural steel connections [7].The present study focuses on the dynamic analysis of steel blast doors and composite blast doors in the form of steel-concrete-steel sandwich panels.
Geometrically, blast resistant doors may be flat or curved, though the doors are usually rectangular in front elevation.The structural resistance of blast doors is primarily provided by both membrane and bending actions.Correct modelling of the nonlinear behaviour of blast resistant doors under severe loading is therefore essential in their safe and cost-effective design.Dynamic analysis of shell structures, among all structural forms, is perhaps the most challenging both mathematically and physically.Mathematically, it is relatively difficult to formulate and obtain solutions for shell structures, particularly when nonlinear dynamic response is involved.Physically, the response of shell structures involves the interaction between the mem-ISSN 1070-9622/03/$8.00 2003 -IOS Press.All rights reserved brane mode and bending mode which, in general, are of different orders of magnitude in stiffness effects.
Analytical approaches do not apply to nonlinear problems in general, and numerical approaches have to be used.With rapid advances in computer hardware, the finite element method has proven to be a versatile and powerful tool for nonlinear dynamic analysis of shell structures.Many shell theories and solution methods have been developed in the last few decades.The finite element formulation of shell elements can largely be categorised into two approaches: (a) direct approach based on a selected shell theory and (b) approach based on degenerated continuum element first proposed by Ahmad et al. [8].A good literature review of shell elements may be found in [9].
As this study involves short-duration loading, explicit dynamic analysis, which avoids the time consuming simultaneous equation solver, is preferred.For explicit analysis, 4-node shell element derived by the degenerated continuum approach is an attractive choice.In contrast, higher-order elements are sensitive to mesh distortion and less stable in an explicit dynamic analysis [10].Reduced integration with single quadrature point is preferred for reason of computational efficiency, but numerical means have to be introduced to suppress hourglass modes.
In blast resistant door designs, steel and concrete materials are commonly used.A combination of steel and concrete would take advantage of high strength of steel and added inertia of concrete in-fill.The sandwich-like structure would require composite elements in the analysis.In this regard, the explicit dynamic approach developed by Koh et al. [11] for 4-node composite shells is adopted herein for the dynamic analysis of blast resistant doors.The constitutive models adopted for steel and concrete are, respectively, Von-Mises yield criterion with hardening and a scalar-damage model (established in the framework of continuum mechanics).Bilinear spring elements are used to handle contact and gap between the door and door frame.Nonlinear effects due to geometry updates are accounted for by means of the updated Lagrangian formulation.

Composite shell elements and dynamic analysis
In this study, a bilinear four-node quadrilateral shell element with single quadrature point [12] is adopted in the finite element modelling of blast resistant doors.In essence, a co-rotational coordinate system is established for each element.The nodal coordinates and ve-locities are transformed to the co-rotational coordinate system.The velocity strains for membrane stretching, bending and shear are computed.This involves straightforward strain-displacement matrices evaluated at the origin of the co-rotational coordinate system.As the element is under-integrated (in plane), a stabilization approach called perturbation hourglass control as proposed in [12] is adopted to suppress the spurious singular modes.
Since blast resistant doors may be composite in material, it is necessary to split the element in several plies that are not necessarily made of the same material.Different plies may have different thicknesses.In this regard, the approach for the explicit dynamic analysis of laminated composite shells [11] is appropriate.It is assumed herein that the bonding between plies remains intact with no delamination even though deformation may be large.The focus of the analysis is the global behaviour (e.g.deflection) of the structure, rather than the local behaviour (e.g.strain variation through the thickness).Therefore, the assumption that the plane section across the thickness remains plane is made.To account for possible large thickness strain, the thickness is updated using log-strain rather than engineering strain.Geometric nonlinear effects due to large deformations are included by means of the updated Lagrangian formulation.
In terms of dynamic analysis, there are two possible approaches: the implicit and explicit integration methods.In the implicit methods, element stiffness matrices have to be computed, assembled and stored in the global stiffness matrix.The equations of motion are converted into a set of algebraic equations by some numerical integration schemes such as the Newmark-Wilson family [13] and then solved by a simultaneous equation solver.Since the system is nonlinear, an iterative scheme such as the Newton-Raphson method is required in the equation solving until certain convergence criterion is met.As the system is nonlinear and coupled, the equation solving process normally forms the bulk of the computational effort.Thus, though the implicit methods with appropriate parameters can be unconditionally stable, the cost per time step is often high.
Explicit methods, on the other hand, do not require the formation of stiffness matrices.The equations of motion in the global co-ordinate system may be written as where U , P , F , M and C are the displacement vector, internal force and external force vectors, mass and damping matrices, respectively.Note that no stiffness matrix is required.With the assumption of lumped mass matrix and mass-proportional damping matrix, the equations of motion are uncoupled.Expressing the acceleration and velocity in terms of displacement, the displacement vector can be solved for the current step without a simultaneous equation solver.The drawback is that explicit methods are conditionally stable; thus a small time step (smaller than a critical value) has to be used.This, however, is not a hindrance in this study because the time step used for impulsive (blast) loading is normally small in order to accurately capture the loading function.In such applications, the overall computational cost for explicit methods is generally less than implicit methods.
In the present study, a simple yet effective explicit method known as the central difference method is chosen [13].Let subscripts n and n − 1 denote the current and previous time steps, respectively.The velocity for the i-th degree of freedom is replaced by the following central difference approximation: where ∆t n denotes time step used at the n-th step (current step) and subscript n + 0.5 denotes the specific time midway between the n-th step and (n + 1)-th step.Similarly, Note that the formulation allows variable time steps.For certain situations where long-term response is needed, the computational time would be reduced substantially with the application of a variable time step scheme.The velocity and acceleration at the current step are given by: where Making use of the above equations, it can be readily shown that Eq. (1) yields To account for gap As the explicit scheme is conditionally stable, the time step has to be less than a critical value to ensure numerical stability.The critical time step is computed based on a safe estimate as follows [14] where ρ is the density, E Young's modulus, v Possion's ratio, h element thickness, and L s the shortest distance between adjacent nodes.The stresses are computed at all the integration points through the thickness.In general, composite layer thicknesses are arbitrary and need not necessarily be symmetrical with respect to the mid-surface of the element.Each layer can be further divided into sub-layers (of the same material), depending on the relative thickness and desired numerical accuracy.In this instance, the commonly used rule of Gauss quadrature is not convenient for integration through the thickness.Furthermore, this rule is not necessarily the most accurate strategy for integration of nonlinear stress through the thickness [15] and it is thus justifiable to use a simpler integration rule.In the present study, the mid-point rule is used, i.e. integration points are the midpoints of sub-layers and quadrature weights are proportional to the sub-layer thicknesses, though not equal in general.For ease of reference, this is termed as the composite integration rule.This is different from the centroidal method as presented in [15], which divides the element thickness into sub-layers of equal thickness.

Material models
In the design of blast resistant doors, steel and concrete materials are commonly used for reason of costeffectiveness though more advanced materials are avail- (A) Steel-air-steel (SAS) door -The door comprises two "skin" plates, stiffened by side plates along the four edges and some stiffener beams (e.g.I or C-channel cross section) at intermediate locations.All the plates and beams are made of steel.(B) Steel-concrete-steel (SCS) door -The door is similar to the SAS door, except that the void is filled with concrete.

Elastoplastic model with hardening for steel
Steel usually exhibits a work-hardening behaviour when stressed into the elastoplastic zone, and the isotropic Von Mises material model is appropriate.The yield function is generally in the form of f y (σ, κ, . ..),where κ is the hardening parameter which depends on the accumulated plastic strain in some manner.The Von Mises criterion for the assumed plane stress state can be written as where σ Y is the yield stress in pure shear, and the yield function (f ) is given by For stress update, the key issue is often the integration of the flow rule in each time step.A commonly adopted approach is the backward Euler method which is stable and usually efficient [16].This involves iterative solutions to achieve convergence by means of, say, the Newton-Raphson scheme which is adopted in this study.Non-iterative schemes are also available if the advantage of high concurrency for the explicit dynamic analysis is to be exploited towards vectorization and parallel computing [11].

Scalar damage model for concrete
The concrete model adopted here is a scalar damage model developed in the framework of the continuum damage mechanics [17].It is relatively simple in formulation and easy to implement in the finite element analysis.More importantly, the strain rate effect on concrete behaviour is accounted for.The computational model is partially verified by an experimental investigation program involving impact tests of concrete specimens.This model is appropriate for steelconcrete composite doors where concrete in-fill serves mainly to increase the inertia (rather than strength) and severe concrete damage is not expected.
In general, the mechanical properties of concrete are enhanced with increasing strain rate.The domain of stress and strain are divided into two parts, one for the low (quasi-static) to intermediate rate ( ε 30s −1 , and the other, from intermediate to high rates.The CEB-1988 recommendations [18]   the rate-dependent concrete material properties under compressive and tensile actions. The rate dependent damage model is formulated using the concept of internal variable to represent the extent of concrete damage.The stress-strain relationship of the concrete model is given as (11) where σ ij and ε ij are the components of the stress and strain tensors, respectively, E ijkl denote the initial moduli, and D is a scalar damage parameter that ranges from 0 (virgin material) to 1 (asymptotic failure).In this concrete damage model, concrete is assumed to remain isotropic up to failure.The positive strains determine the growth of damage that is mainly associated with opening of micro-cracks.The damage consists of both tension and compression components, defined as follows D = α t D t + α c D c (12) where subscripts t and c denote tension and compression, respectively, and α t and α c are non-dimensional functions of the principal strains.More details on the damage model adopted can be found in [17] which covers experimental verification of the model.

Contact/Gap modelling
There are two types of blast resistant doors in practice, namely the air-tight door and blast-tight door.In the former case, rubber gaskets are used between the door and door frames for airtightness, giving rise to a structural gap between the door and the door frame.In the latter case, the door is in contact with the door frame (no gap) when closed.In either case, the contact and gap behaviour is modelled by using nonlinear spring elements in the finite element analysis.The stress-strain relationship of a typical contact/gap spring is illustrated in Fig. 1.The tension part represents the situation when the door and the frame are not in contact, causing no stress in the element.When the door comes into contact with the door frame, the element becomes very stiff under compressive action.

Numerical study and discussion
Instead of using a general purpose software, a specific PC-based software is developed as a research tool to implement the explicit dynamic analysis of blast resistant doors with the main features as illustrated in Fig. 2. In particular, the composite shell element as mentioned earlier is incorporated.The flowchart of the software called BASS (Blast Analysis of Shell Structures) is illustrated in Fig. 3.The software has the option of using Gauss integration rule for non-composite shells or the composite integration rule otherwise.The following numerical examples are presented to verify the formulation and implementation (more numerical examples can be found in [19]).Example 1: Cantilever beam Consider a cantilever beam with rectangular cross section and parameters as follows [12]: length = 0.254 m, width = thickness = 2.54 cm, Young's modulus = 82.8MPa, density = 12.87 kg/m 3 and Possion's ratio = 0.2.The beam is subjected to a suddenly applied pressure loading is analysed using the software BASS.The applied pressure loading is 69 Pa (0.01 psi), and the material is assumed to be linearly elastic.
The present solution agrees well with the published results [12] as presented in Table 1.It is noted that the present solution is closer to the beam element solution as compared to other solutions as summarised in [12].

Example 2: Spherical cap
Clamped all around the edge, a spherical cap is subjected to a suddenly applied pressure of 4.14 MPa (600 psi).The geometric and material properties are as fol- lows [20]: internal radius = 0.5657 m, thickness = 1.041 cm, semi-angle (θ) = 26.27• , Young's modulus = 72.4GPa, density = 2778 kg/m 3 , Poisson's ratio = 0.3, yield stress = 169 MPa, and hardening parameter = 0.By virtue of symmetry, only a quarter of the sphere is considered and is discretised into12 elements (Fig. 4).Both linear and nonlinear analyses are carried out.The problem was analysed by Liu [20] using 8-node degenerated thick-shell elements and an implicit dynamic analysis approach.The same material model is adopted by Liu [20] and in the present study, and no rate effect is considered in this example.Figure 5 shows the vertical displacement at the centre of the sphere for linearly elastic material, in comparison with Liu's result.For the nonlinear case, the Von Mises yield criterion as described earlier is considered but without hardening.The corresponding dynamic response is presented in Fig. 6 and, due to softening effect (yielding), is generally larger than the response for linearly elastic material as shown in Fig. 5.Note that in [20], the results were computed using 6 layers (accurate up to 5th order polynomial in the thickness direction), whereas in the present study, 5 Gauss points (accurate up to 9th order polynomial) are used.Therefore, the present solution should be more accurate than the results given in [20].

Example 3: Numerical study of blast resistant door
In this example, a single-panel rectangular SAS door subjected to blast loading is analyzed.The door considered is made of two 9-mm mild steel skin plates, stiffened by four side plates along the edges and three Cchannels at intermediate locations, as shown schematically in Fig. 7 with the blast pressure) and inner skin plate are modelled by the shell elements.Stiffeners are modelled as beam elements in this study as they deform predominantly in flexural bending, though shell elements may also be used.The door is surrounded by a rigid door frame around all four edges, with two hinges on the left edge.There is an initial gap of 1 mm between door and door frame.There is no initial gap between latch and door frame.An idealized pressure loading as shown in Fig. 8 is used.The displacements at critical points on the outer and inner skin plates are presented in Figs 9 and 10, respectively.The displacements at points A, C and E at the stiffeners are smaller than the other two points (B and D) away from the stiffeners, as expected.It is also observed that the outer skin plate undergoes larger deformation than the inner skin plate.
Numerical results for up to 8 ms are obtained by the software BASS.The maximum displacements of the skin plates are presented in Table 2.These results are found to be generally in good agreement with the results obtained by an implicit finite element analysis software called LUSAS [21], particularly for the two most critical points B and D where the deformations are large.In terms of computational time, BASS based on the explicit approach is found to be about twenty times more efficient than the implicit approach to achieve roughly the same accuracy.

Example 4: Comparison with field test results for blast resistant doors
Field tests of blast resistant doors of both SAS and SCS types were carried out, each under 100 kg bare charge of TNT with a stand-off distance of 5 m.The peak pressure loading is computed by a computer program called CONWEP [22] and its shape is idealized as linear, as shown in Fig. 11.The dimensions of all test specimens are 2.2 m × 1.2 m, as shown in Fig. 12.Each door specimen comprises two skin plates, three C-channel stiffeners, and side plates all around the four edges.A total of eight specimens were blast tested in the field: four SAS doors and four SCS doors.In the case of SCS doors, plain concrete of grade 30 (designed to have compressive strength of 30 MPa at 28 days conforming to British Standard BS 5328 [23]) was used as in-fill.Table 3 shows the geometrical dimensions of the specimens.The specimens were vertically held by clamping their two shorter edges on a specially constructed concrete supporting frame.Displacement gauges were installed behind the door specimen to measure the maximum and permanent deflection of the specimen.The field test results are presented in Table 4 and compared with the numerical results obtained by BASS.
Response measurement in the field is generally a difficult task and particularly so for blast tests where the peak response occurs in a very short time.Some field test results may be inaccurate due to measurement errors or even unavailable.In fact, some results for permanent displacements were not available due to dam- age of instruments during the blast test.Furthermore, it should be noted that the specified blast time history (peak value generated by CONWEP and idealized linearly) used in the numerical simulation studies differs from the real blast in the field tests.In this light, the dynamic analysis by means of software BASS gives reasonably good correlation with the blast test results, as shown in Table 4.The other possible source of discrepancy is the modeling of actual boundary conditions, which play an important role in the displacement response of blast resistant doors.Both the numerical and field test results confirm the effectiveness of the concrete in-fill in that the displacement response of SCS blast doors is considerably smaller than that of SAS blast doors.

Conclusions
This paper illustrates the use of explicit dynamic analysis for the analysis and design of actual blast resistant doors.The main characteristic of the quadrilateral shell element used is that only single quadrature point is needed in the plane of each element.This results in a simplified formulation and low computational cost per element-time step.The shell element allows the extension to treat composite materials using the composite integration rule.
In terms of nonlinearity, three aspects have been included.(A) For material nonlinearity, Von Mises yield criterion with isotropic hardening is adopted for steel, whereas a scalar damage model developed for concrete is used.(B) Geometric nonlinearity due to possible large displacement motion is accounted for by means of Lagrangian formulation.(C) The contact/gap behavior is modeled using bilinear spring elements.
The computational procedure has been coded leading to the software called BASS.The numerical study presented includes comparison with published solutions for some benchmark problems of beams and shells, with another finite element (implicit) software and with field measurements for actual blast resistant doors.In general, the reasonably good correlation validates the numerical procedure and software implementation.The software developed is useful in the parametric studies to optimize the design of blast resistant doors.
are adopted with regards to

Fig. 4 .
Fig.4.Geometry of the spherical shell and its finite element mesh.

Fig. 12 .
Fig. 12. Schematic diagram of blast door specimens for field tests (plan view).

Table 1
Dynamic response of cantilever beam

Table 3
Details of blast resistant doors used in field tests Note: N.A. = Not available due to instrumentation problem.