Structural response and optimization of airtight blast door under gas explosion load

2.1. Basic principle for SPH algorithm

SPH algorithm is proposed by Lucy and GINGOLD [28-29], mainly used for Astrophysics and Universe. Up till now, it has become one of the oldest meshfree particle methods and is applied to deal with the dynamic response of solid structure for the first time in the year 1991 [30]. Nowadays, great advantages have been made in dealing with large deformation problems prompted by a long period of theoretical improvement [31]. SPH is a Lagrangian meshless particle method that can overcome the defects based on grid method. The free surfaces can be traced naturally in the process of simulation by proper distribution of the particles, therefore, it is suitable to follow large surface deformation involving free surface and interfacial [32]. A part simulated with SPH particles and can be defined a contact with FEM mesh in a model is called FEM-SPH coupling algorithm [33]. Such an algorithm can deal with the destruction of the material and large deformation in simulation process effectively, and it has aroused great concern in the application of explosion and shock dynamics.

The basic idea of SPH algorithm is to describe a domain with a series of particles and approximate the field function with integral representation. Prominently, such a method can elucidate the field function and its derivative at each time step based on the current support domain distribution of particles, as well as applying the particle approximation method to the related items of all the field functions.

2.2. SPH algorithm

The core idea of SPH algorithm is to establish a kernel function for kernel estimation of discrete particles so as to calculate various kinetic parameters, transform the conservation equations into integral equations in differential forms, and compute the kernel estimation of all the variables at an arbitrary point. The kernel estimation for function $f\left(x\right)$ at a certain point $x$ in the space can be obtained by integral in the domain $\mathrm{\Omega }$ of $f\left(x\right)$:

where, $x$ is spatial coordinate of the points to be estimated, $x\text{'}$ is spatial point contributing to $x$, $w(x-x\text{'},h)$ is named as the kernel function or smooth function, $h$ is defined as a smoothing length of affected areas of smooth function $w$, $\mathrm{\Omega }$ is integral or support domain. In SPH method, a smooth function satisfies the following conditions:

1) Positive and semi-definited, satisfies $w(x-x\text{'},h)\ge 0$ in the support domain.

2) Compactness, satisfies $w(x-x\text{'})=0$ out of the support domain.

3) Normalization, $\underset{\mathrm{\Omega }}{\int }w(x-x\text{'},h)dx\text{'}=1$.

4) $w(x-x\text{'},h)$ decreases monotonously with the addition of $d$, where $d=‖(x-x\text{'})‖$.

5) If the scale of support domain $h$ approaches to zero, then the function $w(x-x\text{'},h)$ converge towards $\delta (x-x\text{'})$.

6). The function $w(x-x\text{'},h)$ is even and needs to be smoothed sufficiently.

The whole system of SPH method can be expressed by several finite particles with independent quality and independent space, the function of an arbitrary particle can be approximated by a smooth function of all the particles in a support domain that corresponding to the function value of a weighted average, as shown in Fig. 1.

Fig. 1Particle interactions in SPH within the influence domain governed by the kernel function

If the volume $\mathrm{\Delta }{V}_j$ of a particle is applied to replace the infinite body element $dx\text{'}$ at the front integral of the particle $j$, the mass $m_j$ of the particle can be stated as follows:

where, ${\rho }_j$ is density of particle $j$, $N$ is numbers of particles in the support domain.

The continuous SPH integral of $f\left(x\right)$ can be expressed in the following discrete form by particle approximations:

SPH method can be extended to the corresponding field of material under dynamic load. Given the elastic plastic effect of the material, differential statements of SPH in the whole stress-tensor space are deduced as below:

Mass conservation equation:

Momentum conservation equation:

Energy conservation equation:

Position equation of particle motion:

where, $\rho$, $m$, $\sigma$, $v$, $e$ and $w$ are the density, mass, stress, velocity, internal energy and kernel function of a particle, respectively.

Smooth length is a crucial parameter in SPH method for which can affect the calculating efficiency and accuracy directly. If $h$ is too small in value, a lack of enough particles in the support domain will lead to a low accuracy. The calculating efficiency and accuracy will be affected greatly if $h$ is too large. Then, the changeable smooth length can be described as follows:

where, $v$ is velocity of a particle, $h\left(t\right)$ is smoothing length vary with time in the support domain.

The smooth length will increase with the separation of particles in their expansion process, and will decrease with the vicinity of particles in their compression process. The number of particles in the support domain would be stable owing to a changeable smoothing length.

In this paper, we choose Gauss kernel function as the smooth function, which can be expressed as follows:

where, $q=\left|{\overline{x}}_i-{\overline{x}}_j\right|/h$, ${\overline{x}}_i$ and ${\overline{x}}_j$ are the position vectors of the particles $i$ and $j$, respectively.

2.3. FEM-SPH contact algorithm

Fig. 2 shows the contact force added between SPH and finite elements, small solid circles represent the SPH particles, large dashed circle depicts the support domain of SPH particle $i$, small dotted circles means background particles set on the finite element nodes. The background particles have their own SPH properties and can be passively searched by SPH particles. The mass, position, velocity and stress of a particle are updated according to the corresponding finite element nodes. Contact force generates as the distance between a finite element node and a SPH particle is two times of the smooth length. Calculation of the contact force refers to the contact idea of SPH algorithm, with regard to finite element nodes as particles, however, any SPH particle located at finite element node within the support domain will exert contact force on the nodes. Conversely, any finite element node in the support domain will impose contact force on a SPH particle. The calculation process is shown in Fig. 3.

Fig. 2Contact between SPH particles and finite elements

Fig. 3Solution procedure for SPH-FEM contact algorithm

A contact potential is defined in order to calculate the contact force:

where, $N_c$ is the numbers of particles that adjacent to particle $i$ in support domain but belongs to different bodies (finite element nodes for each particle), for particle $i$ (Fig. 1), $N_c=\text{3}$, nodes $n_6$, $n_7$ and $n_8$ will work on the particle $i$. $m_j$ and ${\rho }_j$ are mass and density of particle $j$, respectively. $W$ is the SPH smoothing kernel function, if both $x_A$ and $x_B$ belong to a same body, then $W(x_A-x_B)=0$. $r_{ij}$ is the distance between two particles, $\mathrm{\Delta }{p}_{av}$ is the average distance between different particles. $K$ and $n$ are user-defined parameters, and $K$ is similar to the contact stiffness in the finite element contact algorithm, affected by material properties and impact velocity.

The contact force between a SPH particle and a finite element node can be calculated as follows:

The direction of contact force is determined by the gradient of SPH smoothing kernel function, which can be considered as an external force of the SPH momentum equation and the finite element dynamic formulation, such a force can be calculated as follows:

where, ${\sigma }_i$ is stress tensor of particle $i$, $f\left(x_i\right)$ is contact force exerted on particle $i$. $u$, $M$, and $D$ are the displacement matrix, lumped mass matrix and stiffness matrix of FEM algorithm, respectively.

The integration of SPH and FEM is required to be synchronized in the Lagrange framework, and the data must be transmitted at the same time point. It is demanded that time step used in each step remains the same, therefore, selection for smaller time steps between the SPH and FEM is a good choice, that is:

where, $\mathrm{\Delta }{t}_{SPH}=\alpha h/c$ and $\mathrm{\Delta }{t}_{FEM}\le L/c$ are time step of SPH and FEM algorithms, respectively. $h$ and $L$ are smooth length and the smallest unit size, respectively. $c$ is sound velocity of the material, and $\alpha$ is scale coefficient of time step.

The FEM-SPH contact algorithm is consistent with the solution format of SPH. Accordingly, the contact force can be applied on the relevant particles (unnecessary to define) and nodes as long as the distance between the particle and the node is two times of smooth length. Different from the finite element contact algorithm, the SPH-FEM algorithm can be easily extended to 3D owing to the unnecessity of defining the contact interface and the normal direction at each time step.

2.4. Topology optimization algorithm

Topology optimization, a mathematical approach to optimize the material distribution in a certain region for the given load, constraint conditions and performance indexes and can be used for structure optimization, Structural optimization includes size optimization, shape optimization, topology optimization and topography optimization. The main research field of topology optimization consists of continuum topology optimization and discrete topology optimization. The continuum topology optimization method includes homogenization method, variable density method, evolutionary structure optimization (ESO) method and level set method. The basic idea of topology optimization is to find an optimal topology of the structure that can be transformed into a problem of finding the optimal material distribution in a certain designing area. Topology optimization design mainly includes three basic elements: design variables, constraint conditions and objective functions.

The mathematical model of optimal design can be expressed as follows:

Minimization:

Constraint conditions:

where, $X=(x_1,x_2,...,x_n)$ and $f\left(X\right)$ are design variable and objective function, respectively. $g_j\left(x\right)$ and $h_k\left(x\right)$ are constraint conditions, superscripts $L$ and $U$ are upper limit and lower limit, respectively.

3.1. Geometry and mesh

The structure of an airtight blast door must be strength enough so that can protect a structure against explosion shock wave and instantaneous high temperature once a gas explosion accident occurs. According to the design and development of an airtight blast door of a certain model, an airtight blast door was composed of 4 parts: explosion resistance structure, sealing structure, adiabatic structure and hand wheel. The original anti-knock structure of the original model is a rectangular knock back with a massive stiffened plate. Fig. 4(a) showed the shape and geometric model of the airtight blast door. Solid element is applied to express the entity structure and observation window. The observation window is welded on the top of the door. Table 1 lists the material type and thickness of each component on the structure.

Taking the actual size of the model established above into account, the door hinge is described with rotatable elements mesh and the door shell is meshed with the mapping method and simulated using 163 element, each unit size was 0.01 m. Grid of the door is shown in Fig. 4(b).

Fig. 4Geometry and mesh of the original model

a) Geometry of the original model

b) Mesh of the original model

3.2. Material models and state equation

The process of impact and dynamic under gas explosion is a typical high-strain-rate phenomenon. The mechanical properties of steel and hardened glass will change noticeably when affected by the rate of explosive load. The yield strength of the material will increase substantially under high strain rate. In this paper, the nonlinear plastic material model PLASTIC-KINEMATIC proposed by Symonds and Cowper [34] is employed to simulate the strain rate effect of steel and hardened glass, constitutive curve of the model is described in Fig. 5.

where $\overline{\sigma }$ and $\overline{\epsilon }$ are stress and strain of the material, respectively. $C$ and $P$ are Cowper-Symonds constants, the suggested values for steel are $C=$ 40.4 and $P=$ 5, respectively. $\beta$ is hardening parameter, which means that kinematic hardening if $\beta =0$, but isotropic hardening if $\beta =$ 1. ${\sigma }_0$ is the initial yield stress, $\epsilon \text{'}$ is the strain rate, ${\epsilon }_P^{eff}$ is the effective plastic strain, ${\sigma }_y$ is the yield strength, $E_p$ is the plastic hardening modulus and determined by the following calculation:

where, $E_{tan}$ is the Yield hardening modulus, $E_c$ is the Elastic modulus. The nonlinear plastic material model PLASTIC-KINEMATIC is applicable for metal materials under blast loading. The calculated results are in consistent with the experimental data. The mechanical properties of material parameters are shown in Table 2.

Fig. 5Constitutive curve for steel and hardened glass. Where σ¯ and ε¯ are stress and strain, respectively

3.3. Simulation assumptions

Some further hypotheses need to be made in order to eliminate the impact of various unfavorable factors:

a) Regardless of the changes of material properties and the influence of welding angle on the structure.

b) Material is homogeneous and isotropic.

c) Deformation caused by production and installation can be ignored, and relative movement between different components can be without consideration.

d) The bolt connection is reliable, and pre-stress has no effect on the whole structure.

3.4. Topology optimization and dynamic response for a structure under gas explosion load

In this paper, ANSYS is applied to conduct shape and structure optimization of a structure. There is no doubt that the load (static load and dynamic load) on a structure should be determined firstly, and the next is to establish finite element model of the structure, exert loading, displacement constraints based on the boundary conditions, and define the optimization variables, objective function, constraint functions and optimal parameters and finally submit to the solution. The flow chart of topology optimization and mechanical response analysis for a structure is shown in Fig. 6.

Fig. 6Flow chart for topology optimization and mechanical response of a structure based on topology and FEM-SPH

4.1. Models and conditions

Fluid-solid coupling algorithm can be employed to simulate the explosion of the initial point, shock wave propagates and the interaction between shock wave and structures. However, it is inevitable to generate a large number of units and lead to variation for a long time using such an algorithm [3]. In this work, according to the requirements specification of numerical simulation, the gas explosion load can be simplified as the shock wave induced by a volume of 200 m³ that 100 m away from the structure in a tunnel.

Fig. 7Schematic diagram of impact load

However, explosion shock wave interacts with the structure is a typical transient process, in accordance with equivalent momentum, blast shock wave can be transformed into a triangular shock wave which is evenly applied to the surface of airtight blast door [35]. Here, the triangular shock wave loading curve is determined by three parameters: the peak overpressure, time of pressure increasing and actuation duration of overpressure. The maximum pressure of the triangle shock wave is 0.6 MPa, and the pulse width is 300 ms. The impact load curve is shown in Fig. 7.

Airtight blast door is the main part of a coal mine rescue capsule or refuge chamber, therefore, the damage modes and criterion of airtight blast wall fit in well with that of the coal mine rescue devices, as given in Table 3 [36]. Based on the relative deformation of the main body structure, the damage modes and failure of the key components of the door can be evaluated according to the criteria of the numerical simulation.

As it can be seen from Table 3, the deformation failure condition for the airtight blast door is that the maximum deflection of plate is no less than 2 % or the deformation is larger than 20 mm. Furthermore, the failure condition is that the maximum stress exceeds yield strength of the material, the failure condition for sealing is that the relative displacement of required connection requirements is larger than 1 mm or the deformation is larger than 10 mm.

4.2. Stress field

The stress contours of the door at different moments are shown in Fig. 8, which indicates that the pressure is gradually increased from blue to red.

As it can be seen from Fig. 8, higher stress appears at the door while lower stress on the stiffened plate. Influenced by the bolt connection and rigid connection, it is easy for the plate to generate stress concentration. Therefore, the overall strength of the door would be affected by the reinforcement method.

Fig. 8Stress contours for the original model

The bolt should be rigid enough so as to meet the antiknock performance under gas explosion load in actual working conditions. The safety factor of the door can be defined as $n={\sigma }_s/{\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ , demonstrating that the structural strength of the door does meet the requirements but such a reinforcement method is unreasonable. On the other hand, the original door is a site of materials and needs to be optimized.

4.3. Displacement field

Fig. 9 shows the displacement contour of the door at different moments, the transition of color from blue to red represents the displacement changes from small to large.

As it can be seen from Fig. 9, the most severe serious deformation area locates in the middle part and the stiffened plate of the door. The maximum deformation of the door is less than 20 mm, indicating that the original model door body can meet safety requirements. According to the numerical simulation and performance analysis, the maximum misses stress is about 94.5 MPa, smaller than the material yield limit of 345 MPa, the maximum deformation is 0.85 mm, which does not exceed the limitation of damage failure and sealing failure. Meanwhile, the strength of the door meets the general technical requirements but the stiffened plate needs to be optimized.

Fig. 9Displacement contours for the original model

5.1. Topology optimization

Topology optimization, a novel approach that can assist us to obtain the optimal material distribution or channel in a given design space, so as to meet the requirements of different performance or getting the lightest weight. Some commonly-used topological expression forms and interpolation methods are the homogenization method, density method, variable thickness method and topological function description method. In this study, the topology optimization method based on variable density method of continuum is applied for the original door optimization, meaning that in a given design space, it is necessary to discrete the cell density of retain part on the structure to 1 and the unwanted part to 0 after several topology optimization iterations [37]. The optimal material distribution scheme can be searched in an evenly distributed material designed by topology optimization.

As for the original door, topology optimization and shape optimization based on the variable density method were applied on purpose of decreasing the cost of the door as much as possible under the premise that anti-knock and sealing performance meet the safety requirements. It should be demanded that the minimum overall quality is the optimization objective and the limiting conditions are: the overall maximum stress < 345 MPa and the maximum displacement < 20 mm in the iterative calculation process. To achieve this aim, determining the design space is the first step, while the layout for strengthening ribs is decided by the original structure of the door. The reinforcing plate is the design space for topology optimization. Without considering the specific structure of the reinforcing plate on the inner edge, the design space for strengthening plate on back of the door is shown in Fig. 10, the topology optimization results based on element density distribution is shown in Fig. 11 (the density of red part is 1, and the density of green part is 0).

Fig. 10Design space for topology optimization

Fig. 11Element density distribution for topology optimization

As it can be seen from Fig. 11, according to the topology optimization results, the number of strengthening ribs on the door body should be 2. Both the longitudinal and transverse reinforcing ribs play a connection and auxiliary role for the door. Considering the stiffened plate on the back of the door, the design space for topology optimization would make a difference. Namely, the structure of reinforcing plate on the back has already been identified, therefore it is unnecessary to consider the impact of force exerted on a structure in the topology optimization.

5.2. Static response for the optimization structure

Having finished the topology optimization and shape optimization of the original door, the new model is shown in Fig. 12. The distribution characteristics of stress field and displacement field of the new model are simulated, analyzed and compared with the original model.

The stress contour of the door is shown in Fig. 13, which indicates that the pressure is gradually increased from blue to red.

Fig. 12Geometry and mesh of new airtight blast door after structure design optimization

a) Geometry of the optimization model

b) Mesh of the optimization model

Fig. 13Stress contours for the optimization airtight blast door

As it can be seen from Fig. 13, higher stress appears at the door while lower stress on the strengthening ribs of the door. The maximum stress is 186.8 MPa, influenced by the bolt connection and rigid connection, it is easy to generate stress concentration for the strengthening ribs. Therefore, the overall strength of the door would be affected by the strength of strengthening ribs.

The bolt should be rigid enough so as to meet the antiknock performance under gas explosion load in actual working conditions. Here, the safety factor defined as: $n={\sigma }_s/{\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$, if such a factor of the door is less than 1, demonstrating that the structural strength of the door meet the requirements.

Fig. 14 shows the displacement contour of the door at different moments, the transition of color from blue to red represents the displacement changes from small to large.

Fig. 14Displacement contours for the optimization airtight blast door

As it can be seen from Fig. 14, the deformation of the door reaches a maximum when the time arrives at a certain value. The most serious deformation part of the door locates in the bottom and the upper surface of the door. The maximum deformation of the door is greater than 20 mm, showing that the original model door body satisfies the corresponding safety requirements.

Results indicates that the maximum stress and maximum displacement of optimization model has increase by 97 % and 118 % compared with the original model under the same load, respectively. But the weight of the optimized door is less than the original one and its antiknock performance meets safety requirement for coal mine under the same load, as shown in Table 4.

However, under the real circumstances, the load work on the airtight blast door is a dynamic shock wave, the interaction between gas explosion shock wave and the door is a complicated process, the dynamic strain and stress need to be simulated in order to reveal the dynamic response of the door, hence, we will conduct a simulation using SPH particles to replace the gas explosion shock wave.

6.1. Models and conditions

The SPH-FEM model of airtight blast door and the tunnel are established by ANSYS/LS-DYNA. The tunnel to be simulated is assumed as a semi-closed pipe model, and the corresponding size is described in Fig. 15(a). As it can be seen from the Fig. 15(b), SPH particles are established to simplify the explosion shock wave caused by gas explosion. The model is divided into two parts, part 1 (from A to B) consists of 4500 SPH particles to substitute the same gas explosion load induced by a length 28 m of methane/ air mixture and full of 200 m³ methane/air mixture with 9.5 % concentration [3]. However, in this work, we develop a new part 2 (from B to C), which includes the airtight blast door and shock wave transmission section.

Fig. 15Length and the cross-sectional area of the tunnel

a) Cross-section of the tunnel

b) SPH equivalent model of tunnel and airtight blast door

Fig. 16Finite element mesh of the door and tunnel

a) Cross-section of the model

b) Mesh of FEM-SPH model

Taking the actual size of the model established above into account, the door hinge is described with rotatable elements mesh while the door shell is meshed with the mapping method and simulated using 163 element, each unit size is 0.01 m. The tunnel is meshed with hexahedral grid and simulated with 164 element. Grid of the door and the tunnel is shown in Fig. 16. Then related projects are submitted to LS-DYNA and the stress filed as well as displacement field are analyzed as below.

6.2. Stress field

The stress contours of the door at different moments are shown in Fig. 17, which reveals that the pressure increases gradually from blue to red. The equivalent stress of the cabin door increases gradually as the shock wave propagates forward. Higher stress appears at the door while lower stress on the cabin door and emergency door during the whole process of SPH particles interact with the door.

Fig. 17Stress contour for the airtight blast door

a)$t=$ 1.89 ms

b)$t=$ 8.36 ms

c)$t=$ 15.2 ms

d)$t=$ 20.3 ms

As it can be seen form Fig. 17, stress concentration appears at the impact face of the door. The maximum stress concentration appears at the bottom surface of the door. At the moment of 12.54 ms, the maximum stress on element 30402828 is 284 MPa, as shown in Fig. 18(a), but the maximum stress of the door does not surpass the material’s yield strength. The maximum deformation of the unit for the door is about 7.8 mm, as shown in Fig. 18(b). However, no obvious plastic deformation appears and the stress of elements at the other portions of the door does not exceed the material yield strength. When the overpressure of gas explosion reaches about 0.6 MPa, the strength of the overall structure and main parts are in the elastic range, illustrating that it does not reaches the yield stress and the deformation can recover. Because of the short time of explosive shock, stress concentration and plastic zone do not lead to the failure of the whole door.

Fig. 18Maximum stress result for the airtight blast door

a) Stress contour for the maximum stress element

b) The maximum stress element time history

6.3. Displacement field

Fig. 19 shows the displacement contour of the door at different moments, the transition of color from blue to red represents the displacement changes from small to large. Deformation on the impact face of door is the most obvious owing to be subjected to the explosion shock wave and the dynamic pressure directly. As the shock wave gradually propagates forward, the tunnel wall reflects shock wave and the displacement of the door increases owing to superimposed deformation and continuous force exerted by explosion blast wave.

Fig. 19Displacement contour for the airtight blast door

a)$t=$ 1.89 ms

b)$t=$ 8.36 ms

c)$t=$ 15.2 ms

d)$t=$ 20.3 ms

As it can be seen from Fig. 19, deformation of the node changes smoothly and the largest deformation emerges on the impact surface of the door. Finite element analysis results show that the displacement of node 30499751 reaches the maximum. The maximum absolute displacement value reaches about 7.89 mm at the moment of 20.14 ms, as shown in Fig. 20(a). The maximum displacement of element’s stress-time curve is extracted and the maximum stress value of the door is about 297 MPa, as shown in Fig. 20(b), illustrating that it does not reaches the yield stress and the deformation can recover. However, the absolute maximum displacement value of the door exceeds 20 mm and could guarantee the safety of the rescue according to code for design of mine rescue equipment.

According to the static and dynamic analysis of the optimization blast door, it could guarantee the safety requirement according to the code for design of mine rescue equipment and can be good for safety production. But the maximum displacement and maximum stress of the optimization model using FEM-SPH algorithm under dynamic load is about 3 times larger than those of the door under static load.

Fig. 20Maximum displacement result for the airtight blast door

a) Displacement contour for the maximum displacement element

b) the maximum displacement element time history