Large Liquefied Natural Gas (LNG) tanks are prone to damage during strong earthquakes, and accurate seismic analysis must be performed during the design phase to prevent secondary disasters. However, the seismic analysis of large LNG tanks is associated with high computational requirements, which cannot be satisfied by the calculation efficiency of traditional analytical techniques such as the Coupled Eulerian–Lagrangian (CEL) method. Thus, this paper aims to employ a less computationally demanding algorithm, the Smoothed Particle Hydrodynamics-Finite Element Method (SPH-FEM) algorithm, to simulate large LNG tanks. The seismic response of a 160,000 m3 LNG prestressed storage tank is evaluated with different liquid depths using the SPH-FEM algorithm, and simulation results are obtained with excellent efficiency and accuracy. In addition, large von Mises stress at the base of the tank indicates that strong earthquakes can severely jeopardize the structural integrity of large LNG tanks. Therefore, the SPH-FEM algorithm provides a feasible approach for the analysis of large liquid tanks in seismic engineering applications.

1. Introduction

Natural gas is a reliable source of energy that is used globally to meet growing energy demands. With the increasing consumption of natural gas over the past few years, Liquefied Natural Gas (LNG) tanks have become a major component of urban infrastructure. As a result, the scale of tank construction has also increased, but this change has been associated with various safety risks. LNG storage tanks have high seismic risks compared to traditional buildings because they can lead to secondary disasters, such as explosions and environmental pollution, that result in significant property damage or loss of life. In addition, LNG tanks are more vulnerable to earthquakes owing to their low redundancy, low ductility, and low energy-dissipating capacity compared to conventional structures [1]. For example, the destruction of an LNG tank during the 1964 Japan earthquake caused fires and explosions that resulted in serious societal losses and pollution. Another example was the 1976 Tangshan earthquake in which the bottom ring of a storage tank buckled and led to liquid leakage [2]. Thus, it is of paramount importance to investigate the seismic performance of LNG tanks to ensure their structural safety.

To analyze the dynamic performance of large LNG tanks, it is essential to understand fluid–structure interaction (FSI) and dynamic performance of the structure. Different from conventional structures, seismic analysis of liquid storage tanks must account for fluid–structure interaction as a result of inertial earthquake loading and hydrodynamic pressure. The seismic design of storage tanks was initially established on the rigid wall model proposed by Housner [3] in 1957, and many subsequent models have improved on this theory. However, many tanks designed according to this model have suffered damage during earthquakes [2]. Consequently, Veletsos and Yang [4, 5] proposed a single-degree-of-freedom model for a flexible tank wall that assumed that the tank vibrated under a given bending mode. In 1981, Haroun and Housner [6] developed a more realistic model that considered the dynamic interaction of liquid and the tank wall. This method is known as the coupled vibration effect model, and it divides the liquid into three components: a convection component, a flexible impulsive component, and a rigid impulsive component. In addition, Hu et al. [7] presented a method for mass addition using user subroutines to study the seismic response of a large cylindrical thin-walled storage tank based on the theories of Housner and Veletsos. Liquid mass was added to the solid wall without considering liquid sloshing. Hariri-Ardebili et al. [8] compared the endurance time analysis (ETA) method to the time-history analysis (THA) and incremental dynamic analysis (IDA) method to investigate the seismic performance of a high-rise telecommunication tower. Hadj-Djelloul and Djermane [9] used the three-dimensional finite element technique to study the seismic response of perfect and imperfect elevated water tanks. They also investigated the effect of geometric imperfection on the dynamic performance of elevated water tanks.

Based on applicable theories regarding fluid–structure interaction, since the 1950s, there have been several studies on the dynamic analysis of LNG storage tanks. Hwang [10] investigated the dynamic response of an LNG tank using a combination of Boundary Element Method (BEM) and Finite Element Method (FEM) programs. Wu [11] performed time-history analysis of LNG storage tanks with FEM software and found that deformation of the inner tank was effectively prevented by adding a prestressed concrete wall and an insulating layer. Christovasilis and Whittaker [12] investigated the seismic response of a conventional and isolated vertical cylindrical LNG tank by applying finite element analysis to mechanical models. The results obtained from two numerical models were in good agreement and demonstrated that a mechanical model could be used with confidence for the preliminary analysis and design of conventional and isolated LNG tanks. Sun et al. [13] derived and calculated the dynamic response of base-isolated LNG storage tanks considering soil–structure interaction, the effect of which was proven to be insignificant. Chen et al. [14] investigated the dynamic response of a typical 160,000 m3 LNG prestressed concrete outer tank under impact loading, and types of impact damage were determined based on dynamic response results including stress, displacement, energy, and critical impact velocity. Li et al. [15] conducted numerical dynamic analysis of LNG storage tanks and concluded that filling the space between the outer and inner tank could improve seismic performance with little effect on the sloshing height. Du [16] adopted the Coupled Eulerian–Lagrangian (CEL) analysis technique in ABAQUS to simulate fluid–structure interaction [17] in LNG tanks, but the CEL method exhibited liquid leakage and low simulation efficiency.

The aforementioned numerical simulation methods are very computationally intensive due to the sheer size of large LNG tanks. Given the recent rise in large LNG tank construction, it is essential to employ a fast analysis method that can simulate the seismic response of LNG tanks in practical engineering applications. Thus, an efficient smoothed particle hydrodynamics-finite element method (SPH-FEM) algorithm [18], which has not been previously applied to large land-based LNG tanks, is proposed in this paper. SPH-FEM methods retain the advantages of FEM in simulating FSI problems while adopting the advantages of SPH in modeling fluid behaviour.

SPH is suitable for modeling fluid behaviour because it is a mesh-free method that discretizes a liquid field into particles [19]. A kernel function is used to approximate differential equations, and the function value at any point is then expressed by neighboring nodes based on a local approximation. The SPH algorithm can overcome the weakness in mesh-based methods regarding fluid simulation by effectively addressing free surfaces and moving interfaces in addition to large deformation of materials. As a result, after its development by Lucy [20] and Gingold [21] in 1977, the SPH method has been successfully adopted for studies in continuous solid mechanics and fluid mechanics to solve problems related to structural failure, large deformation, and liquid sloshing [22]. Huang et al. [23] verified the feasibility of liquid sloshing analysis using the SPH method by obtaining results that corresponded to theoretical analysis. Liu et al. [24] performed a three-dimensional numerical simulation of free surface liquid sloshing in a prismatic tank subjected to a sinusoidal excitation using the SPH algorithm, and the numerical results exhibited good agreement with experimental data. Shao et al. [25] also presented an improved SPH algorithm to model liquid sloshing dynamics, and the numerical results agreed well with the experimental observations. The SPH algorithm has also been widely applied to simulate embankment flow [26], wave breaking [27], and explosions [28]. Therefore, the SPH algorithm is a versatile analysis method that can be applied to a variety of structural problems, including the fluid component of large LNG tanks.

Coupled with a FEM method for the solid component, an SPH-FEM algorithm greatly improves the calculation accuracy and efficiency for fluid–structure impact problems. The coupled SPH-FEM method was first proposed by Attaway et al. [29] using a master–slave algorithm to account for the contact between FE elements and SPH particles. Kalateh and Koosheh [30] used the SPH-FEM method to simulate the interaction between a convergent-divergent nozzle and cavitating flow, and the results showed that the behaviour of liquid and vapor at the interface matched other numerical and experimental methods. Liang and Chen [31] applied the SPH-FEM method to soil–structure interaction problems during the seismic analysis of a rectangular underground structure, and results indicated that the distribution and magnitude of seismic earth pressure were influenced by the magnitude of soil deformation. Fragassa et al. [32] analyzed the effect of air in FSI problems by coupling FEM and SPH. They concluded that air had little effect on stress prediction, but it had a greater effect on the behaviour of the structure after the initial fluid impact. In this study, the SPH-FEM algorithm provides an accurate and efficient method for the analysis of a large LNG tank under dynamic earthquake loading.

Section 2 of this paper presents an overview of the SPH-FEM methodology. Section 3 verifies the SPH-FEM algorithm by comparing it to the traditional Coupled Eulerian–Lagrangian (CEL) method on the analysis of a cubic water tank under sinusoidal excitation. Subsequently, Section 4 analyzes the dynamic response of a 160,000 m3 LNG storage tank under three earthquakes with the same site conditions. Section 5 concludes the study and reiterates the advantages of the SPH-FEM algorithm as an innovative, accurate, and efficient method to simulate large LNG tanks.

2. Methodology

2.1. SPH Theory

The formulation of the SPH method is often divided into two parts: the integral representation of the field function and the particle approximation [33]. The integral representation of the function in the SPH algorithm begins with the following identity [18]:where is a function of the three-dimensional position vector x, and the Dirac delta function is given by the following equation:

If a smoothing function replaces the Delta function kernel , the integral representation of can be written as follows:where is the smoothing kernel function and is the smoothing length that defines the area of influence of the smoothing kernel function . The smoothing kernel function should satisfy the following conditions (see equations (4)–(6)):(1)The normalization condition states that(2)The Delta function property is observed when the smoothing length approaches zero(3)The compact condition implies thatwhere , which is a constant related to the smoothing function of the point at , defines the effective (nonzero) area of the smoothing function.

The basic governing fluid dynamics equations are based on three fundamental laws of conservation: conservation of mass, momentum, and energy. The governing equations for dynamic fluid flows can be written as a set of partial differential equations using the Lagrangian method; this system of partial differential equations is the famous Navier–Stokes equations, the governing equations of which are given as follows (see equations (7)–(9)):(1)The continuity equation:(2)The momentum equation:(3)The energy equation:

A continuous density particle approximation method is employed to approximate the density. This approximation method is obtained by applying the SPH approximation concept to transform the continuity equation. Different forms of density approximation equations can be obtained by applying different conversions and operations to (7). One such approach is to apply the SPH approximation to the velocity divergence; then, the particle density in (7) can be evaluated in the following form of a gradient:

The particle approximation method for the momentum equation is similar to the abovementioned continuous density method, and some transformations are required. The momentum equation approximation can be derived in differential form based on different transformations. Equation (11) can be obtained using the SPH particle approximation method to transform the gradient term of the momentum equation (equation (8)):

There are several approximation expressions of the work performed by pressure; thus, the internal energy calculation for the work performed by pressure has many alternative forms. The following form (12) is the most commonly used:

Time stepping in ABAQUS is explicit and is limited by the Courant Condition as shown below (equation (13)):where is the local speed of sound. The time step used in all simulations is 1.0 × 10−6 s, and the convergence of the simulation is guaranteed [34].

2.2. Coupled SPH-FEM

The main concept in SPH-FEM coupling is centered around the relationship between the two algorithms. In each time step, the velocity and displacement of FEM nodes are transferred to particles. Simultaneously, the coordinate and velocity of particles are transferred to FEM nodes. While particles supply boundary conditions for FEM, FEM maintains the continuity of the particles by preventing the boundary effect [35]. In this study regarding fluid–structure interaction, velocity and displacement fields on the solid–fluid interface are exchanged in order to couple the interaction forces between SPH and FEM. A schematic computation flowchart of the SPH-FEM coupling method is displayed in Figure 1.

3. Example Verification

Compared to the SPH-FEM method, the coupled Eulerian–Lagrangian (CEL) method developed by Du [16] is another analysis technique that is suitable for most fluid–structure interaction problems. Qiu et al. [36] applied the CEL method to large deformation problems in geotechnical engineering such as pile sinking and ship grounding to demonstrate its reliability. Meng et al. [37] used the CEL method to simulate the influence of pile-sinking construction on the soil around the pile under undrained conditions. In this section, a simple example is provided to compare the accuracy and efficiency between the SPH-FEM and the CEL method using ABAQUS software [34].

3.1. Description of the Problem

A cubic water tank is taken as a model and is designed according to the following geometric specifications: the side length is 1 m, the wall thickness is 10 mm, and the liquid height is 0.8 m. The tank model is composed of concrete with the following material properties: density , Young’s modulus E = 3.25 × 109 N/m2, and Poisson’s ratio . The behavior of a material experiencing liquid sloshing is extremely complex, and an equation of state is used to describe the pressure, volume, and energy of the material. The Newtonian fluid constitutive material formulation is used for water in this section, and the Mie–Grüneisen equation of state, which describes the volumetric strength and density of a material, is defined in Table 1. The isotropic pressure in the Mie–Grüneisen equation of state is defined as follows:where H is the Hugoniot curve, Em is the internal energy per unit mass, is the volumetric compressibility, and c is the sonic speed of the material. The coefficients s and c can be described by the linear relationship between shock velocity and particle velocity:where is shock wave velocity and is particle velocity. This equation of state describes the fluid motion.

The equation of state of provides a model for the hydrodynamic material, and the properties of water are listed in Table 1.

The water-filled models constructed based on the two algorithms are shown in Figures 2 and 3. The models both use 1600 C3D8R-type concrete elements, and the SPH-FEM algorithm model uses 22000 PC3D-type liquid elements. In the SPH-FEM algorithm, the ratio of liquid element size to the side length of the container is 5:270. The SPH method also uses a cubic spline as the interpolation polynomial. Smoothing length is an important parameter and has significant influence on the accuracy of the prediction [38, 39]. By default, ABAQUS calculates a smoothing length at the start of the analysis so that an element is associated with approximate 30 to 50 particles on average. Since the smoothing length remains constant throughout the analysis, depending on whether the behavior in the model is expansive or compressive, the average number of particles per element can either decrease or increase [34]. In the CEL model, a total of 80,000 EC3D8R-type liquid elements are used, including 20,000 elements in the initial position. The Eulerian elements in the CEL method must be finely meshed over the entire tank trajectory and allow liquid to move freely in the grid, because the quantity and quality of the mesh affects simulation accuracy. Gravity is applied to the entire model , and the tank is subjected to a horizontal excitation given by .

3.2. Comparison of Results between the SPH-FEM and CEL Method

The results of the two algorithms during one period are compared based on three factors: sloshing wave pattern, stress distribution on container wall, and computational efficiency. Both methods can accurately reflect the sloshing wave pattern (as shown in Figure 4), but the CEL method exhibits liquid leakage, especially at the bottom of the water tank. The CEL method is also not applicable when the liquid volume fraction in the Eulerian element is less than 0.5, because some elements are no longer contained within the container wall [34].

Using the two algorithms, the von Mises stress distribution on the container wall at 0.8 s and 1.6 s in one period is shown in Figure 5. The maximum von Mises stress and the von Mises stress distribution on the walls calculated using the two methods are in good agreement. Thus, the liquid element size in the SPH-FEM method can be scaled up based on an optimal ratio of 5 : 270 and applied to a large LNG tank while maintaining calculation accuracy.

As shown in Table 2, the simulation time for the CEL method is twice that of the SPH-FEM method, indicating that the SPH-FEM algorithm is much more efficient than the CEL method. In summary, both methods can accurately simulate stress on tank walls and the liquid sloshing pattern, but the CEL method exhibits liquid leakage. Fluid leakage is a common problem in CEL method and has not been solved so far. Therefore, the SPH-FEM method has a significant advantage over CEL in the analysis of large liquid storage tanks.

4. Seismic Analysis of a 160,000 m3 LNG Tank

4.1. Project Overview

This study employed a 160,000 m3 LNG prestressed storage tank, which was adapted from an LNG technical manual as shown in Figure 6 [40]. The inner diameter of the outer tank is 82 m, the wall height is 38.55 m, the wall thickness is 0.8 m, and the dome thickness is 0.4 m. The tank also has components to secure the moisture-proof linings and roof-bearing rings.

The inner tank has a total of 10 layers as shown in Figure 6, each of which has a height of 3.55 m. For example, Layer R1/24.9 refers to the first layer from the bottom with a thickness of 24.9 mm. These layers gradually thicken from the top to the bottom as given in Table 3.

The LNG storage tank consists of four parts: LNG, an inner tank, an insulation layer, and an outer tank. The inner tank is composed of 9% Ni steel. The tank exhibits excellent low-temperature resistance [41, 42]. The tank is also characterized by good weldability and low susceptibility to cold cracks. The insulation layer is composed of expanded perlite and resilient felt. The expanded perlite clings to the inner wall of the outer tank, and the outer wall of the inner tank is fitted with resilient felt, which prevents the perlite from settling and provides elastic properties for the perlite. This elasticity is necessary because the tank shrinks due to temperature changes. The prestressed concrete external tanks are constructed to increase the overall safety of the tank in the event of accidental LNG leakage. Adapting to the advice of Huang [43], the material of the outer tank is 60 MPa high-strength concrete. The material characteristics of the tank and LNG are listed in Tables 4 and 5, respectively.

4.2. Finite Element Model

ABAQUS software [44, 45] is used to establish the finite element model (FEM) as shown in Figure 7. The tank model is relatively large, and the associated interactions are extremely complex. Rebar rods are embedded within the tank wall using the “Embedded” option in ABAQUS, which allows the rebar elements to work in conjunction with each other. The outer tank and the insulation layer are composed of C3D8R mesh elements, and the inner tank is composed of S4R shell elements. Shell elements are selected because the thickness of the inner tank is much smaller than its radius of curvature. The LNG is discretized into particles using the aforementioned equation in Section 3. The element size is chosen based on the small tank in Section 3 for which the sloshing frequency can be accurately calculated. The ratio of liquid element size to container size in the large LNG tank is identical to the parameters of the small tank in Section 3.

19T15S steel cables with a tensile strength of 1860 MPa are embedded in the protective cover of the tank. A total of 122 steel cables are vertically anchored at the bottom and top of the concrete wall, and 220 circumferential steel cables are anchored in separate half-circles on four vertical supporting columns arranged at 90° [40].

The basic principles of applying prestress in ABAQUS [44] include various approaches, such as the model predictive control (MPC) method, falling temperature method, initial stress method, and single rebar element method. This study adopts the commonly utilized falling temperature method that applies prestress by decreasing temperature. The initial temperature is obtained with the following formula:where N is the tension control force of the prestressed rebar, which should not exceed 75% of the standard value of the strength; denotes the thermal expansion coefficient of prestressed rebar; E represents the elastic modulus of prestressed rebar; and A is the cross-sectional area of prestressed rebar.

The bottom boundary of the model is fixed, and the interaction between soil and structure is not considered. The static condition of the LNG tank is primarily calculated based on self-weight to ensure that the LNG liquid reaches steady state. Natural frequencies are important parameters in the design of LNG prestressed storage tanks. There are three primary ways to extract eigenvalues in ABAQUS: the AMS eigensolver, the subspace iteration method, and the block Lanczos algorithm. Here, the first five frequencies of the empty tank are calculated with the Lanczos algorithm as 3.753 Hz, 4.902 Hz, 5.088 Hz, 5.097 Hz, and 5.204 Hz.

4.3. Records Selection

In China, seismic codes for large LNG storage tanks are based on the seismic codes of building structures [46] with specifications for the design of large oil storage tanks. The normative theoretical model is developed from the Housner rigid two-particle theoretical model considering the influence of the elastic tank wall. In addition, Chinese seismic codes classify design objectives into three tiers based on the degree of damage. The design of large LNG tanks belongs to the second tier meaning repairable damage without collapse. Furthermore, the Chinese seismic code divides seismic sites into four classes based on soil stability, shear wave velocity, and thickness of overburden. Soil strength ranges from strongest to weakest from Class I to Class IV with Class I having the strongest soil. In this paper, the large LNG tank is located in a Class II site. Based on guidelines for seismic wave selection [46], three waves are selected and scaled to ensure compatibility with the design spectrum at the chosen site. The details for three ground motions are shown in Table 6, and the acceleration response spectra for ground motions as well as the code acceleration response spectra are shown in Figure 8.

4.4. Results

The entire simulation using the SPH-FEM algorithm is conducted using an ordinary personal computer. The processor is an Intel® Core™ i7-8700 CPU @ 4.10 GHz, and the installation memory is 16.00 GB. The von Mises stress distribution of the inner tank alone under hydrostatic conditions at different LNG liquid levels, without the effect of insulation layer and outer tank, is shown in Figure 9.

Among the three seismic waves, Taft produces the greatest von Mises stress on the outer and inner tank at different LNG liquid levels. Therefore, Taft is chosen as the critical wave for this analysis. The stress distribution on the outer tank exhibits a symmetrical trend about the x- and z-axes. Figure 10 shows the respective maximum stress of the outer tanks with LNG levels of 25%, 50%, 75%, and 100% at the instant corresponding to peak displacement. When the tank is 100% filled with LNG, the maximum stress of the outer tank is 51.35 MPa, which is less than the compressive strength of concrete; thus, the structure is safe. In addition, the maximum stress of the outer tank is located at the dome for liquid levels of 25% and 50%, but it is located at the bottom of the tank for liquid levels of 75% and 100%. This phenomenon occurs because gravity plays a more important role than pressure with LNG liquid levels of 25% and 50%. However, as liquid level increases, pressure at the bottom of the tank increases at a much faster rate relative to the dome. Furthermore, while the maximum stress of the outer tank is nearly identical at 25%, 50%, and 75%, it increases significantly from 75% to 100% liquid level. In general, close attention should be paid to the dome and the bottom of the tank during the design of large LNG tanks to ensure structural safety.

The respective maximum stress of the inner tank with LNG liquid levels of 25%, 50%, 75%, and 100% is shown in Figure 11. The fact that dynamic conditions generate much greater stress than static conditions deems dynamic analysis critical for LNG tanks. When the tank is 100% full, the maximum stress of the inner tank exceeds 500 MPa and jeopardizes the safety of the structure considering that the yield strength of 9% Ni steel is between 500 MPa and 600 MPa [47]. Due to the severity of consequences after large LNG tank failure, it is imperative to reduce the stress of the inner tank to ensure its safety. At the same time, an early warning system should be implemented to control the maximum liquid level. Different from the outer tank, the stress of the inner tank changes significantly with the liquid level. When the liquid level is at 25% and 50%, maximum stress of the inner tank is below 300 MPa. At a liquid level of 75% and 100%, maximum stress exceeds 300 MPa and 500 MPa, respectively. With an increase in the amount of LNG, the stress at the base of the inner tank gradually increases to notable levels. Therefore, more attention should be paid to the control of liquid level in the dynamic analysis and seismic design of large LNG storage tanks.

According to postprocessing results, three representative paths, as shown in Figure 12, are selected to analyze the variations of the von Mises stress with the height of the tank. Stress distributions on the outer and inner tanks with respect to height at different liquid levels are illustrated in Figures 13 and 14, respectively. It can be observed that the stress distribution on the outer tank changes more gradually compared to that of the inner tank.

In addition, stress values along different paths are different for a given liquid depth, but the overall trends are the same. When the liquid level is less than 50%, the stress on the outer tank first increases and then decreases with an increase in height. When the liquid level exceeds 50%, the stress on the outer tank initially decreases, then increases, and decreases again. Moreover, stress in the lower portion of the tank is much larger than stress in the upper portion when the liquid level is greater than 50%. In addition, a greater LNG liquid volume leads to greater stored energy and greater stress at the base of the tank as shown in Figures 13 and 14. In reality, the bottom of an LNG tank is prone to buckling failure during a strong earthquake, which is consistent with the stress analysis results.

When comparing outer and inner tanks, it is observed that different liquid amounts can produce similar stress levels at a height of approximately 10 m. For the outer tank, the height that correspond to this stress level gradually decreases from Path A to Path C. In contrast, for the inner tank, the height that corresponds to this stress level is nearly identical across all paths.

The displacements of the highest point on the outer wall of the tank filled with 25%, 50%, 75%, and 100% LNG are shown in Figure 15. In a Class II site, the displacements for all three seismic waves are within the same range. For the same seismic wave, variations in displacement with respect to time are nearly identical for each liquid volume. In addition, the location for peak displacement nearly coincides with the peak seismic wave acceleration. The maximum displacement only reaches 4.16 mm when the tank is filled with 100% LNG, and the seismic performance of the LNG tank is satisfactory based on this displacement calculation.

Figure 16 illustrates that the base shear values of the tank are in the same range under different seismic waves for each liquid volume. For the same seismic wave, variations in base shear with respect to time as well as peak base shear forces follow the same general trend for each liquid volume. The maximum base shear reaches when the tank is filled with 100% LNG under the Taft wave; the minimum base shear is when the tank is filled with 25% LNG under the Taft wave. These values demonstrate that the base shear of a LNG storage tank is dependent on the liquid level in the tank.

4.5. Validation

To verify the SPH-FEM simulation results under earthquake loading, they are compared with CEL simulation results in terms of tip displacements and base shear. The maximum tip displacement and the maximum base shear using the SPH-FEM and CEL method under El Centro, Taft, and Northridge seismic waves are shown in Figures 17 and 18. The tip displacement of the outer wall and the maximum base shear calculated using the two algorithms are in good agreement. Therefore, the feasibility of using the SPH-FEM algorithm for the seismic analysis of large LNG tanks is verified.

5. Summary and Conclusions

In this study, the SPH-FEM method is compared with the CEL method to demonstrate the advantage of SPH-FEM in terms of efficiency and accuracy. Then, the SPH-FEM algorithm is used to analyze the dynamic response of a 160,000 m3 LNG prestressed storage tank subjected to three earthquake waves in a Class II site. The main conclusions drawn from the numerical simulations can be summarized as follows:(1)Under static conditions, the von Mises stress increases at a linear rate with an increasing liquid volume. Under dynamic conditions, the von Mises stress increases at a nonlinear rate.(2)The stresses produced on the outer and inner tanks under earthquake loading are very similar at a height of approximately 10 m for each liquid volume. The maximum tip displacement of the outer tank is 4.16 mm, and the maximum base shear is .(3)The maximum stress of the inner tank with 100% LNG liquid level exceeds 500 MPa under the chosen seismic waves and raises concerns for its structural safety. Development of an improved structural system and a warning mechanism is of paramount importance to the design of LNG tanks at high liquid levels.(4)The SPH-FEM algorithm is accurate and more efficient than the CEL method. Moreover, the SPH-FEM algorithm exhibits excellent capability for simulating large storage tanks under a reasonable time to provide a reference for the design and construction of large LNG tanks.

The largest LNG prestressed storage tank in China is currently 200,000 m3. The sheer size of the tank requires a significant amount of meshing during seismic simulation, which is almost impossible to accomplish on an ordinary computer using a mesh-based method. The SPH-FEM algorithm provides a feasible method to simulate extremely large LNG tanks on a personal computer, which can help designers obtain its seismic performance in an acceptable time range without significant reliance on hardware capabilities. Ongoing work is being performed to optimize the design of large LNG tanks under different earthquakes using the SPH-FEM algorithm and will be reported at a later time.


:Domain of integration
:Kernel function
h:Smoothing length
:Density, kg/ m3
t:Time step size
N:Number of particles
:Gravitational acceleration
m:Mass of particle
cs:Local speed of sound
Em:Internal energy per unit mass
:Volumetric compressibility
:Direction of coordinates
:Direction of coordinates
:Poisson’s ratio
:Small positive number
i:The ith component of a vector
j:The jth component of a vector.

Data Availability

The data of ground motions used to support the findings of this study have been deposited in the PEER Ground Motion Database.


Any opinions, findings, conclusions, and recommendations expressed here are those of the authors and do not necessarily reflect the views of the sponsors.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.


This research was supported by the China Scholarship Council, the Natural Science Foundation of China (Grant no. 51738007), and the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery (201606290) Program. The authors would like express their appreciation to Professor Ozden Turan and Professor Eren Semercigil for their input.