Fluidization Prediction of Liquefiable Cargo: Based on FINN and UBCSAND Constitutive Models
https://doi.org/10.1007/s1180402200266x

Abstract
Fluidization mechanism of liquefiable cargo is similar to that of sand liquefaction. But the cargo is subjected to more complex external loads during sea transportation. Based on two different constitutive models, numerical predictions were conducted in this study. And the effects of ship motions including motion acceleration, frequency and relative density of cargo on fluidization was investigated. By comparing with available the experimental data, validation is carried out. Results show that the roll motion is the most important one related to cargo fluidization. When the motion acceleration increases, the possibility of cargo fluidization increases. The higher of the cargo density the lower of the possibility of cargo fluidization. The effect of frequency on cargo fluidization is not unique, and it exists a critical value. The fluidization behavior of cargo could be described both by UBCSAND and FINN models. And the comparisons were discussed and summed up.
Keywords:
 Liquefiable cargo ·
 Fluidization ·
 FINN ·
 UBCSAND ·
 Ship motion ·
 Relative density
Article Highlights• A more comprehensive numerical model was built and the effects of cargo density and different ship motions on the fluidization of the cargo were investigated.• Performance of the two widely used liquefaction constitutive models FINN and UBCSAND on the prediction of the cargo fluidization were presented and compared. 
1 Introduction
The liquefiable cargo composed of a certain proportion of fine particles and moisture. During sea transportation, when its moisture content exceeds a certain limit and under external loads such as ship motions, it may fluidize and form a fluid state. The stability of the ship would be weakened due to the fluidized cargo (Zhang et al. (2017; 2019); Abbas et al. 2019). From 1986 to 2019, a total of 25 cargo fluidization accidents occurred, resulting in the sinking of 19 ships (Sakar et al. 2020). The common causes of the accidents are cargo fluidization and the cargo movement in cargo holds. To ensure the safety of sea transportation of liquefiable cargo, it is necessary to study the factors influencing the fluidization of cargo.
At present, experimental tests and numerical methods have been adopted to carry out the research. An indoor small shaking table was used by Jian et al. (2015) to study the microscopic behavior during bulk iron ore fluidization. Under action of vibration, the iron ore concentrate compacts and the moisture between the particles would squeeze out. And moisture migration is the key inducement of cargo fluidization. The fluidization characteristics of iron ore concentrate were studied by Zhou et al. (2013) based on model experiments. Results show that acceleration and moisture content are the two key factors affecting cargo fluidization. A scaled cargo hold model was used by Munro and Mohajerani (2016) to conduct vibration tests of cargo samples with different moisture contents. They found that the cargo with a moisture content lower than transportable moisture limit (TML) may also fluidize. In the above studies, due to limitations of scaled model experiments, the actual stress distribution of cargo could not be reproduced.
In terms of numerical simulation, the UBC3DPLM constitutive model in PLAXIS software was used by Ju et al. (2016) to assess cargo fluidization potential. They reported that as the saturation degree, motion frequency and amplitude increase, the cargo is more likely to fluidize. To evaluate the potential of cargo fluidization under rolling motion, the UBCSAND constitutive model in FLAC software was used by Daoud et al. (2018). Results show that the cargo with a lower relative density is more prone to fluidization and the cargo slope is more prone to have plastic deformtion than other areas. The FINN model was used by Hao (2019) to simulate the fluidization of saturated iron ore concentrates at various frequencies. It is reported that there is critical frequency for the excess pore water pressure of iron ore concentrates. Ding (2020) adopted the FINN model to conduct a numerical simulation on the influence of physical parameters on the fluidization of cargoes. To ensure the safe transportation of liquefiable cargo, the internal friction angle should be greater than 32°, Poisson's ratio should be greater than 0.41, initial pore ratio should be less than 1.48, and permeability coefficient should be less than 6×10^{−6} m/s. The PFC software was used by Zhou et al. (2016) to simulate cargo fluidization. They found that the main reason for cargo fluidization is the rise of the water level. Of the above research, the external loads were simplified and a single sway or roll motion was considered. For the problem of sand liquefaction, vertical vibration is also important. For example, Vasiliki et al. (2016) studied the influence of vertical ground motion on liquefaction of saturated sand. Results show that when the deviator stress generated by resonance is strong enough, sand liquefaction can also happen.
In addition, the acceleration, frequency of the motion and relative density of the sample are also the key factors for liquefaction. For example, it is reported by Vargehese and Iatha (2014) that the increase of acceleration and frequency can cause sand to liquefy. And the sand with higher relative density is less likely to liquefy. Based on the 1 g shaking table test, Wang et al. (2020) suggested that when the relative density of the sample continues to rise, its liquefaction resistance improved gradually. Similarly, acceleration and frequency of ship motion are also the key factors influencing cargo fluidization. Considering the variety of ship motion acceleration and frequency in different navigation areas, the safety load thresholds were proposed by Wu et al. (2020).
To sum up, numerical simulation research on the influence of ship motion forms on cargo fluidization is still needed, especially in the actual sizes. In addition, the selection of appropriate constitutive model is also the key problem to complete the work. In this study, firstly, to predict the fluidization of the cargo in a real sized cargo hold, a more comprehensive numerical model was built. Compared with previous studies, the wall zones were established and the friction conditions between cargo and wall surface were taken into account and permeable boundaries were added for the top surface of the cargo. And this is closer to reality. Furthermore, the effects of cargo density, motion acceleration and frequency on the fluidization of the cargo were investigated. Secondly, the performance of the two widely used liquefaction constitutive models FINN and UBCSAND on the prediction of the fluidization of the liquefiable cargo were presented and compared. The research could provide some reference for safe sea transportation of liquefiable cargo.
2 Theoretical background and constitutive models
2.1 Theoretical background
The finite difference method is used to divide the solution area into quadrilateral grids (Wang and Xing 1995), as shown in Figure 1. The grid is composed of nodes and zones. Usually, 'point' variables such as displacement, velocity, and acceleration are stored in nodes, while 'domain' variables such as density, porosity, stress and strain are stored in zones.
In the calculation process, the node will be affected by the resultant force formed in its surrounding area at each time step. The node will move when the resultant force is not equal to 0. The strain of the zone is generated by the movement of the node, which further generates a new stress increment. The resultant force acting on the node will be updated by the new stress increment, so that the calculation enters the next cycle. It will continue until the resultant force reaches the convergence condition. The calculation cycle is shown in Figure 2.
Two equations and four variables are mainly involved in a calculation cycle. First, the velocity of the node under the action of the unbalanced force is solved by the following equation:
$$ \frac{{\partial v_i^l}}{{\partial t}} = \frac{{F_i^l(t)}}{{{m^l}}} $$ (1) where F_{i}^{l}(t) is the unbalanced force component of the l node in the
direction at time t, and m^{l} is centralized quality of the l node. Secondly, the zone strain increment at a certain time step is solved by velocity, the formula is as follows:
$$ \Delta {e_{ij}} = \frac{1}{2}\left({{v_{ij}} + {v_{ji}}} \right) \cdot \Delta t $$ (2) where Δe_{ij} is the tensor of strain increments, i, j = 1, 2, v_{ij} is node velocity, and Δt is time step.
Then, the zone stress increment is solved by strain increment, and the formula is as follows:
$$ \Delta {\sigma _{ij}} = f\left({\Delta {e_{ij}}, {\sigma _{ij}}, \cdots } \right) $$ (3) where f (…) is constitutive function, Δσ_{ij} is stress increment, and Δe_{ij} is strain increment.
Finally, the unbalanced force of the node in the next cycle is calculated by the principle of virtual work. For details, see the help document (Itasca 2011).
From the numerical calculation process, it can be seen that the constitutive equation plays a key role for numerical results. Two widely recognized liquefaction constitutive models were adopted in this study, the FINN and the UBCSAND models. At present, they are widely used in the numerical simulation of freefield seismic liquefaction and model box centrifuge liquefaction. With the aid of the FINN model, the expected simulation results of material liquefaction were obtained by Wang (2019), Banerjee et al. (2017) and Noui et al (2019). In terms of UBCSAND model, the sample liquefaction behavior is well reproduced in FLAC software by Ecemis (2013), Beaty (2018), Dashti and Bray (2013). For the prediction of soil liquefaction, both the FINN model and the UBCSAND model were used by Chou et al. (2021) to simulate a cyclic undrained direct simple shear (DSS) test by a singleelement. They reported that, different from the UBCSAND model, the FINN model cannot consider the bananashaped stress–strain path and the butterflyshaped stress path observed in the laboratory test. But both models can still provide reasonable simulations of the excess pore pressure accumulation during cyclic loads. Therefore, for the prediction of cargo fluidization, both the FINN and the UBCSAND models were selected and compared.
2.2 FINN constitutive model
The FINN model is a strain pore pressure model, which is essentially a dynamic pore pressure rise model. And it was added to the MohrCoulomb constitutive model. The MohrCoulomb constitutive model is composed of MohrCoulomb failure criterion, tensile failure criterion and fluid plastic dynamic incremental equation. The stressstrain relationship is reflected by the flowplastic dynamic incremental equation and updated under different states.
According to the incremental equation, the elastic incremental law of the principal stress space can be expressed as:
$$ \begin{array}{l} \Delta {\sigma _1} = {S_1}\left({\Delta \varepsilon _1^e, \Delta \varepsilon _2^e, \Delta \varepsilon _3^e} \right) = {\alpha _1}\Delta _{{\varepsilon _1}}^e + {\alpha _2}\left({\Delta \varepsilon _2^e + \Delta \varepsilon _3^e} \right)\\ \Delta {\sigma _2} = {S_2}\left({\Delta \varepsilon _1^e, \Delta \varepsilon _2^e, \Delta \varepsilon _3^e} \right) = {\alpha _1}\Delta _{{\varepsilon _2}}^e + {\alpha _2}\left({\Delta \varepsilon _2^e + \Delta \varepsilon _3^e} \right)\\ \Delta {\sigma _3} = {S_3}\left({\Delta \varepsilon _1^e, \Delta \varepsilon _2^e, \Delta \varepsilon _3^e} \right) = {\alpha _1}\Delta _{{\varepsilon _3}}^e + {\alpha _2}\left({\Delta \varepsilon _2^e + \Delta \varepsilon _1^e} \right) \end{array} $$ (4) where S_{i} is the elastic strain increment equation in the i direction of the principal stress. Δε_{i}^{e} is the elastic strain increment in the i direction of principal stress. α_{1} and α_{2} is material constant, related to bulk modulus K and shear modulus G, it can be expressed as:
$$ {\alpha _1} = K + \frac{4}{3}G $$ (5) $$ {\alpha _2} = K  \frac{2}{3}G $$ (6) When the material undergoes shear or tensile failure, the stress increment law needs to be plastically modified. The stress after shearing plastic flow can be expressed as:
$$ \begin{array}{c} \sigma _1^P = \sigma _1^e  \lambda \left({{\alpha _1}  {\alpha _2}{N_\emptyset }} \right)\\ \sigma _2^P = \sigma _2^e  \lambda {\alpha _2}\left({1  {N_\emptyset }} \right)\\ \sigma _3^P = \sigma _3^e  \lambda \left({{\alpha _2}  {\alpha _1}{N_\emptyset }} \right) \end{array} $$ (7) $$ \begin{array}{c} \lambda = \frac{{{f_s}\left({\sigma _1^e, \sigma _3^e} \right)}}{{\left({{\alpha _1}  {\alpha _2}N\phi } \right)  \left({{\alpha _2}  {\alpha _1}N\phi } \right)N\phi }}\\ {N_\phi } = \frac{{1 + \sin \phi }}{{1  \sin \phi }} \end{array} $$ (8) In the same way, the stress after tensile failure can be expressed as:
$$ \begin{array}{c} \sigma _1^P = \sigma _1^e  \left({\sigma _3^e  {\sigma _t}} \right)\frac{{{\alpha _2}}}{{{\alpha _1}}}\\ \sigma _2^P = \sigma _2^e  \left({\sigma _3^e  {\sigma _t}} \right)\frac{{{\alpha _2}}}{{{\alpha _1}}}\\ \sigma _3^P = {\sigma _t} \end{array} $$ (9) where σ_{i}^{e} is the principal stress value when failure is not considered. σ_{i}^{P} is the stress considering the plastic flow after entering the failure zone. σ_{t} is the strength of extension. ϕ is frictional angle.
In addition, the relationship between the increase of dynamic pore pressure and plastic volumetric strain is as follows:
$$ \Delta u = {\bar E_r} \cdot \Delta {\varepsilon _{vd}} $$ (10) where Δu is the increase of pore water pressure, E_{r} is rebound modulus and Δε_{vd} is the plastic volumetric strain increment.
The Byrne model is adopted in this study to calculate the plastic volumetric strain increment.
$$ \frac{{\Delta {\varepsilon _{vd}}}}{\gamma } = {C_1}\exp \left({  {C_2}\frac{{{\varepsilon _{vd}}}}{\gamma }} \right) $$ (11) $$ {C_1} = 7600{\left({{D_r}} \right)^{  2.5}} $$ (12) $$ {C_2} = \frac{{0.4}}{{{C_1}}} $$ (13) where γ is the shear strain. C_{1} and C_{2} are model constants. D_{r} is relative density of sample.
2.3 UBCSAND constitutive model
UBCSAND is an effective stress plasticity model and it was used to directly evaluate the response of the soil particle skeleton under external loads. The elastoplastic stressstrain relationship of the soil skeleton is expressed as follows:
$$ \left\{ {\Delta {\sigma ^\prime }} \right\} = \left[ {{D^{{e_P}}}} \right]\{ \Delta \varepsilon \} $$ (14) where {Δσ'} is the soil skeleton stress increment vector.
and {Δε^{P}} are the strain increment vector, the elastic strain increment vector and the plastic strain increment vector of the soil skeleton, respectively. is the elastoplastic stiffness matrix of the soil skeleton. Strain in UBCSAND model is divided into elastic strain and plastic strain, and plastic strain is divided into plastic shear strain and plastic volume strain. The plastic shear strain increment is calculated according to the hyperbolic relationship of stress ratio. And the plastic volume strain increment is calculated according to its flow law of the plastic shear strain increment. UBCSAND is mainly divided into two parts: elastic response and plastic response. See the description document of UBCSAND for detains (Beaty and Byrne 2011).
In the elastic response, assuming that the elastic component is isotropic and represented by the shear modulus
and bulk modulus . $$ {G^e} = K_G^e \cdot {P_a} \cdot {\left({\frac{{{\sigma ^\prime }}}{{{P_a}}}} \right)^{ne}} $$ (15) $$ {B^e} = \alpha \cdot {G^e} $$ (16) where K_{G}^{e} is shear modulus. P_{a} is atmospheric pressure. σ' is mean effective stress. ne is the elastic shear modulus index. α dependents on the Poisson's ratio.
In the plastic response, the plastic shear strain increment dγ^{P} is related to the increment of stress ratio dη.
$$ {\rm{d}}{\gamma ^P} = \frac{1}{{{G^P}/{\sigma ^\prime }}} \cdot {\rm{d}}\eta $$ (17) $$ {G^P} = G_i^P \cdot {\left({1  \frac{\eta }{{{\eta _f}}} \cdot {R_f}} \right)^2} $$ (18) where G^{P} is the plastic shear modulus. G_{i}^{P} is the plastic shear modulus under low stress. η_{f} is the stress ratio at failure. R_{f} is the stress ratio at failure.
According to the flow law, the calculation formula of plastic volumetric strain increment is as follows:
$$ \Delta \varepsilon _v^p = \left({\sin {\mathit{\Phi }_{cv}}  \sin {\mathit{\Phi }_{mob}}} \right) \cdot \Delta {\gamma ^P} $$ (19) where Φ_{cv} is the constant volume friction angle. The definition of Φ_{mob} is defined by the maximum stress ratio η before the current half loading cycle, where Φ_{mob} = arcsin (η_{max}).
In addition, the material liquefaction response is predicted by coupling of soil skeleton and pore fluid. Under undrained conditions, assuming that the volumetric strain of the equivalent fluid is equal to the volumetric strain of the soil skeleton (Wang et al. 2002).
$$ \Delta {\varepsilon _v} = \Delta {\varepsilon _x} + \Delta {\varepsilon _y} = \Delta \varepsilon _v^f $$ (20) where Δε_{v} is the skeleton volume strain. Δε_{v}^{f} is the equivalent fluid volume strain.
Therefore, the pore fluid pressure increment Δu is expressed as:
$$ \Delta u = \left({{B_f}/n} \right)\Delta {\varepsilon _v} $$ (21) where B_{f} is the fluid bulk modulus. n is the porosity of the soil skeleton.
2.4 Verification
Take sand liquefaction as verification (Banerjee et al. 2017), and the excess pore pressure is calculated and compared. The calculation grid is shown in Figure 3, the twodimensional geometric model has a height of 650 mm and a width of 900 mm. Consistent with the instrument layout in the shaking table test, pore water pressure is monitored for zones at a depth of 0.3 m (monitoring point 2) and 0.5 m (monitoring point 1).
The numerical prediction of liquefaction is divided into two stages: static and dynamic analysis. In the static analysis, by assigning material properties to the mesh (see Table 1), the initial stress field of the sample induced by gravity alone was established.
Material Parameter Value Sand Density, ρ (kg/m^{3}) 1 600 Elasticity modulus, E (Pa) 2.938e6 Poisson's ratio, v 0.3 Angle of friction, φ (°) 32 Porosity, n 0.4 Permeability, k (m^{2}/(Pa^{*}s)) 1e8 Relative density, D_{r} 50% Density(kg/m^{3}), ρ 2 000 Model box Bulk modulus, K (Pa) 2e8 Shear modulus, G (Pa) 1e8 Water Density, ρ (kg/m^{3}) 1 000 Bulk modulus, K (Pa) 2e8 In the dynamic analysis, the FINN and UBCSAND constitutive models were added separately. The constitutive model parameters are shown in Table 2. The water permeability condition of the sample surface was set by a fixed pore pressure and saturation value. A sinusoidal velocity history was applied to the boundaries of the hold to simulate the horizontal loading process. The velocity equations are as follows:
$$ v = \left\{ {\begin{array}{*{20}{l}} {t/2 \cdot 2{\rm{ \mathsf{ π} }} \cdot A \cdot f \cdot \sin (2{\rm{ \mathsf{ π} }} \cdot f \cdot t), }&{0 < t < 2}\\ {2{\rm{ \mathsf{ π} }} \cdot A \cdot f \cdot \sin (2{\rm{ \mathsf{ π} }} \cdot f \cdot t), }&{\ \ \ \ t \geqslant 2} \end{array}} \right. $$ (22) Model Parameter Value UBCSAND Relative density, D_{r} (%) 50 N_{1, 60} SPT 11 Elastic shear modulus number, K_{G}^{e} 964.4 Stress exponent for elastic bulk modulus, m_{e} 0.5 Elastic bulk modulus number, K_{B} 675.1 Stress exponent for elastic shear modulus, n_{e} 0.5 Plastic shear modulus number, K_{G}^{P} 450.1 Stress exponent for plastic shear modulus, n_{P} 0.4 Critical state friction angle, Φ_{cs}(°) 33 Peak friction angle, ϕ_{peak}(°) 35.2 Failure ratio, R_{f} 0.89 Model parameter, m_{}hfac1 1.0 Model parameter, m_{}hfac2 1.0 Relative density, D_{r} (%) 50 ρ Density(kg/m3), 1 600 Elasticity modulus, E (Pa) 2.938e6 Poisson's ratio, v 0.3 Angle of friction, φ (°) 32 FINN porosity, n 0.4 Permeability, k (m^{2}/(Pa^{*}s)) 1e8 Model constant, C_{1} 0.429 9 Model constant, C_{2} 0.930 4 where v is velocity. A is amplitude. f is frequency. t is time.
To reproduce the shaking table motion acceleration of 0.35 g at a frequency of 2 Hz, A is set as 0.022 m and f is 2 Hz. To use the two constitutive models, the specific values of the model parameters (see Table 2) should be calculated firstly. For the FINN model, the parameter values are acquired based on Table 1, and the additional parameters C_{1} and C_{2} are calculated from the relative density of the sample though equation (12) and (13). For the UBCSAND model, the parameter values are calculated by the empirical formulas through the relative density of the sample. The calculation formulas are shown in Table 3 (Kutter et al. 2020).
Parameter Calculation formula N_{1, 60} (D_{r}/15)^{2} K_{G}^{e} K_{G}^{e} = 21.7 × 20 ×(N_{1, 60})^{0.333} K_{G}^{P} K_{G}^{P} = K_{G}^{e}×(N_{1, 60})^{2} × 0.003 + 100 n_{P} 0.4 m_{e} 0.5 Φ_{cs}(°) 33 K_{B} K_{B} = 0.7 × K_{G}^{e} ϕ_{peak}(°) ϕ_{Cs} + N_{1, 60}/5 n_{e} 0.5 R_{f} 1 − N_{1, 60}/100 Comparison of experimental and numerical results is shown in Figure 4. At the two monitoring points, results of shaking table test (Banerjee et al. 2017) show that the excess pore water pressure has a rapid rise and then go into a stable phase. And the excess pore pressure of the point 2 tends to stabilize earlier than that of the bottom. These features can be well captured by the two different constitutive models. Still, there are some differences between the experimental and simulation results, mainly are the difference of the time for the pore water pressure begins to stabilize and the averaged pore pressure value. The possible reasons can be divided into the following two parts:
(1) In the simulation, the averaged pore water pressure of an area of 0.1 m × 0.1 m is monitored, while in the experiment the data is recorded at a specific point.
(2) At present, the liquefaction constitutive parameters are calculated by empirical formula, additional test corrections are required for more accurate values.
In the test, initial mean effective stress σ_{m}' is used as the threshold of sand liquefaction, which is calculated by the initial vertical effective stress σ_{v}' and coefficient of static lateral pressure K_{0}.
$$ \sigma _m^\prime = \frac{{\left({1 + 2{K_0}} \right)\sigma _v^\prime }}{3} $$ (23) The experimental results show that both the excess pore water pressures of middle and bottom areas are slightly greater than the mean effective stresses. Thus, liquefaction occurred in the middle and bottom of the sand.
Similarly, liquefaction threshold is adopted to judge the liquefaction happens or not for the simulation zones. And its value in the simulation is close to that in the experiment. The vertical effective stress of the element at the depth of 0.3 m and 0.5 m is 2.7 kPa and 4.8 kPa. And under the same K_{0} of 0.47, the mean effective stress is 1.746 kPa and 3.104 kPa, which are close to the experiment results.
For the two constitutive models, results of pore water pressure have deviations under the same conditions. Different from the FINN model, UBCSAND model is the new comprehensive liquefaction constitutive model. And the FINN model is a simple liquefaction constitutive model based on MohrCoulomb model. Same as the experimental results, liquefaction occurred in the middle and bottom of the sand when using the UBCSAND model. The FINN model cannot capture the phenomenon at the bottom of the cargo. Thus, UBCSAND model can capture the liquefaction state of samples better than FINN model.
3 Simulation cases and result analysis
3.1 Model setup
With reference to Daoud et al. (2018), the numerical calculation grid of actual cargo hold in this study is shown in Figure 5. The load size of the liquefiable cargo is having a width of 25 m and a height of 10 m, and three monitoring zones (53, 16), (53, 11) and (53, 5) are set in the interior area from top to bottom. The pore water pressure values of the three different areas of cargo are recorded and analysed.
In static analysis, the initial stress field induced by gravity alone was established. The model parameters are shown in Table 4. In dynamic analysis, two liquefaction constitutive models were added, and the parameters are shown in Table 5.
Cargo Cabin Parameter Value Parameter Value D_{r} (%) 65 80 90 ρ (kg/m^{3}) 7 900 ρ (kg/m^{3}) 1 650 1 700 1 720 K (Pa) 2e11 E (Pa) 1.7e7 1.9e7 2e7 G (Pa) 7e10 v 0.3 0.3 0.3 φ (°) 36.9 38.9 40.5 n 0.585 0.542 0.515 k (m^{2}/(pa^{*}s)) 1e8 1e8 1e8 FINN UBCSAND Parameter Value Parameter Value D_{r} (%) 65 80 90 D_{r} (%) 65 80 90 C_{1} 0.223 1 0.132 8 0.098 9 1, 60 19 29 37 C_{2} 1.792 8 3.012 1 4.044 5 Secondly, the velocity history was applied to the boundary nodes of the cargo hold to realize the three different motions. Since the twodimensional cargo hold was considered in this present study, so the motions are sway, heave and roll, separately. The implementation of the motion is shown in Figure 6.
Considering different motion parameters and cargo relative density, the specific simulation cases are shown in Table 6. The acceleration and frequency values are changed by the amplitude and period of the velocity function. Considering constant acceleration, the relation between amplitude and frequency is as follows:
$$a = 4{{\rm{ \mathsf{ π}}} ^2}{f^2}A$$ (24) Model Motion Acceleration (g) Frequency (Hz) Density (%) FINN sway 0.1  0.5 5 65 0.1 1  5 65 0.1 5 65, 80, 90 heave 0.1  0.6 1 65 0.1 1  9 65 0.1 1 65, 80, 90 roll 0.1  0.4 5 65 0.1 1  5 65 0.1 5 65, 80, 90 sway 0.1  0.15 1 65 0.1 1  5 65 0.1g 1 65, 80, 90 UBCSAND heave 0.1  0.4 1 65 0.1 1  9 65 0.1 1 65, 80, 90 roll 0.1  0.15 1 65 0.1 1  5 65 0.1 1 65, 80, 90 where a is the acceleration. f is the frequency. A is the amplitude.
Furthermore, the excess pore water pressure ratio is used as the basis to judge whether the fluidization happens or not. When its value is approximately equal to 1, the cargo is regarded as fluidized. The formula for excess pore pressure ratio is presented as follows:
$${r_u} = \frac{{\Delta u}}{{\sigma _{{\rm{ini}}}^\prime }}$$ (25) where Δu is excess pore water pressure. σ_{ini}' is the initial vertical effective stress.
The formula for shear strength is presented as follows:
$${\tau _f} = (\sigma  \mu)\tan \phi = {\sigma ^\prime }\tan \phi $$ (26) where τ_{f} is shear strength. σ is total stress. μ is pore water pressure. σ' is effective stress. ϕ is Internal friction angle.
Shear strength is used to describe the ultimate ability of a skeleton to resist shear failure. With the increase of the excess pore pressure ratio, the shear strength of the cargo decreases. Therefore, the value of r_{u} should be paid more attention to.
3.2 FINN constitutive model
3.2.1 Sway motion
Variation of the excess pore pressure ratio and its maximum value at different accelerations are shown in Figure 7. Changes of the excess pore water pressure at the middle monitoring point are significantly affected by acceleration. For different accelerations, the excess pore pressure ratio gradually rises to its maximum first, and then continues to decrease. And the larger of the acceleration, the higher of the excess pore pressure ratio. In addition, after smoothing the excess pressure ratio curve and comparing with the maximum value, it is found that with the increase of acceleration, the maximum values of excess pressure ratio in each region are approaching to 1.0, and the possibility of cargo fluidization increases. The maximum value at the top area descends later because of a shorter steady calculation time under high acceleration conditions.
For the frequency, taking the monitoring point in the middle of the cargo as an example, the variation of the excess pore pressure ratio and its maximum value are shown in Figure 8. Results show that the excess pore pressure ratio also rises first and then drops as frequency increases. In addition, the excess pore pressure ratio reaches its maximum when the frequency is around 3 Hz. It indicates that, the possibility of cargo fluidization is bigger when the motion frequency is about 3 Hz.
For the relative density, also takes the monitoring point in the middle of the cargo as an example, the variation of the excess pore pressure ratio and its maximum value are shown in Figure 9. Results show that the excess pore pressure ratio rises rapidly and then drops slowly. When the relative density is smaller, the corresponding excess pore pressure ratio is higher. In addition, the maximum excess pore pressure ratio in each region gradually decreases as the relative density of the cargo increases. Therefore, cargo with higher relative density is not prone to fluidization.
3.2.2 Heave motion
Taking the monitoring point in the middle of the cargo as an example, the variation of the excess pore pressure ratio and its residual value at different accelerations are shown in Figure 10. As the loading direction is consistent with gravity, the ratio of excess pore pressure fluctuates greatly and the positive value is about 2 times of the negative. The increase of the heave acceleration directly expands the fluctuation range of the excess pore water pressure asymmetry. In addition, due to the maximum excess pore water pressure would be overestimated during such large fluctuations. Therefore, the buffer function is added to the motion equation, and the residual pore water pressure is used for comparison. The buffer function is as follows.
$$v = [(30  t)/10] \cdot 2{\rm{ \mathsf{ π} }} \cdot A \cdot f \cdot \sin (2{\rm{ \mathsf{ π} }} \cdot f \cdot t)$$ (27) where v is velocity. A is amplitude. f is frequency. t is time, 20 ≤ t ≤ 30.
Results show that with the increase of acceleration, the residual pore pressure of the middle and bottom increases. However, the residual pore pressure value of top monitoring point decreases under the high acceleration. Compared with the middle and bottom, the top area is less constrained and deforms larger under high acceleration. That would lead to a faster dissipation, so the residual value of the top area has a downward trend under high acceleration conditions.
Taking the monitoring point in the middle of the cargo as an example, the variation of the excess pore pressure ratio and its residual value at different frequencies are shown in Figure 11. The excess pore pressure ratio rises firstly and then falls. In addition, there are critical values for the residual excess pore pressure ratio, corresponding to 3 Hz and 8 Hz, respectively.
Taking the monitoring point in the middle of the cargo as an example, the variation of the excess pore pressure ratio and its residual value at different relative densities are shown in Figure 12. Results show that under different relative densities, the excess pore pressure ratio also increases firstly and then decreases. In addition, the possibility of cargo fluidization decreases as the relative density of cargo increases.
3.2.3 Roll motion
For different accelerations, taking the monitoring point in the middle of the cargo as an example, the variation of the excess pore pressure ratio and its maximum value are shown in Figure 13. The excess pore pressure ratio at the middle gradually increases firstly and then becomes stable. When the acceleration is larger, the excess pore pressure ratio reaches its maximum value earlier. In addition, when the acceleration increases a little, the maximum excess pore pressure ratio of each area of the cargo would get close to 1.0, and fluidization occurs in all areas of the cargo.
For different frequencies, taking the monitoring point in the middle of the cargo as an example, the variation of the excess pore pressure ratio and its maximum value are shown in Figure 14. The excess pore pressure ratio increases gradually and then decreases. It decreases faster at higher frequencies. In addition, the maximum excess pore pressure ratio increases as the frequency increases, and the possibility of cargo fluidization increases accordingly.
For different relative densities, taking the monitoring point in the middle of the cargo as an example, variation of the excess pore pressure ratio and its maximum value are shown in Figure 15. The excess pore pressure ratio increases firstly and then decreases. In addition, the maximum excess pore pressure ratio decreases with the increase of the relative density of the cargo. Therefore, cargo with higher relative density is not easy to fluidization.
3.2.4 Analysis
The above results show that acceleration and relative density have the same effect on liquefaction under different motion forms. The increase in motion acceleration would increase the possibility of cargo fluidization, as shown in Figures 7, 10 and 13. The increase of the cargo relative density would reduce the possibility of cargo fluidization, as shown in Figures 9, 12 and 15. However, the effect of frequency is not the same at different motion forms. In sway and heave motions, there is a critical frequency, as shown in Figures 8 and 11. In rolling motion, the increase of frequency would increase the possibility of cargo fluidization, as shown in Figure 14. Essentially, the accumulation of pore water pressure is related to the shrinkage of the pores, caused by the external load. The shrinking pores squeeze the pore water, resulting in an increase of the pore water pressure. The higher of the motion amplitude, the more of the settlement. Therefore, with the increase of acceleration, the pore water pressure can accumulate higher values. For the frequency, only when the frequency closes to the resonance of the system, can it significantly promote the accumulation of pore water pressure. For the double peaks phenomenon of the pore water pressure in heave motion, it is related to multistage resonance. That is, the resonance state of the system is corresponding to more than one single frequency. For the relative density, the denser of the cargo, the smaller of the pore volume. There is less space for a further compression under load. So, under higher load conditions, it is difficult for highdensity cargo to accumulate pore water pressure. And these were also found in the study of soil liquefaction. For example, in the work of Vargehese and Iatha (2014), Wang et al. (2020) and Chen et al. (2018).
To further study the effect of ship motion forms on cargo fluidization, numerical simulations of the cargo with a relative density of 65% under the heave, sway and roll motion at 0.1 g and 5 Hz were performed respectively. Due to the inconsistency of the distance between the boundary nodes and the center of rotation, the acceleration of each node is different in the roll motion. The farthest node having the maximum acceleration, and the nearest node having the minimum acceleration. To compare more clearly, two sets of simulation schemes were adopted in the roll motion, setting 0.1 g to the nearest node or the farthest node acceleration.
Figure 16 shows the excess pore pressure ratio and its residual value under different motion forms. It can be seen that, for the three different motions, the excess pore pressure ratio at the top rises rapidly first and then decreases gently. Compared with the residual excess pore pressure ratio, it shows that the possibility of cargo fluidization of the heave motion is smaller than the others two. Secondly, under the condition of minimum acceleration of 0.1 g, the residual excess pore pressure ratio of roll motion is obviously higher than that in sway motion. And under the condition of maximum acceleration of 0.1 g, the value in roll motion is close to that in sway motion. Therefore, the risk of cargo liquefaction is greater under rolling motion. What is more, referring to Figure 13, not only the top and middle part of the cargo are extremely easy to fluidize, but also the bottom of the cargo can reach the fluidization threshold after a small increase of motion acceleration. Thus, the potential risk of cargo fluidization under the rolling movement is the highest. Therefore, the roll motion should be paid more attention to for preventing cargo fluidization.
3.3 UBCSAND constitutive model
3.3.1 Sway motion
The variation of the excess pore pressure ratio and its maximum value are shown in Figure 17. The excess pore pressure ratio rises gradually first and then becomes gentle. When the acceleration is bigger, the corresponding excess pore pressure ratio is also higher. Furthermore, the top and middle parts of the cargo both are likely to liquefy, and their maximum excess pore pressure ratios are close to 1.0. For the bottom part of the cargo, the maximum excess pore pressure value increases as the acceleration increases, and the possibility of cargo fluidization becomes higher. This is the same as the result of the FINN constitutive model, but the accumulation of excess pore water pressure can be significantly affected by a small increment of acceleration in the UBCSAND constitutive model.
For the frequency, selecting the middle monitoring point as a representative, the variation of the excess pore pressure ratio and its maximum value are shown in Figure 18. When the frequency is relatively high, the excess pore pressure ratio is smaller. And it continues to rise and cannot reach a stable state in the present simulation. Regard to the maximum excess pore pressure ratio, its value also decreases as frequency increases, which means that the possibility of cargo fluidization becomes smaller. For the top and the middle parts of the cargo, the excess pore pressure ratio still can rise to 1.0 when the frequency is 2 Hz. And that means the middle and upper parts of the cargo are prone to liquefy under same motion frequencies. This is different from the conclusion of FINN constitutive model. There are two trends in the influence of frequency in FINN, the maximum pore pressure increases/decreases with the increase of frequency, and there is only one downward trend in UBCSAND.
For the relative density, selecting the middle monitoring point as a representative, the variation of the excess pore pressure ratio and its maximum value are shown in Figure 19. At high relative density, although the excess pore pressure ratio continues to increase, the corresponding excess pore pressure ratio is relatively small. For different parts of the cargo, the maximum excess pore pressure ratio decreases as the relative density increases. This is similar to the conclusion of the FINN constitutive model. However, the excess pore pressure ratio of the top and middle areas can still rise to around 1.0 when the relative density of the cargo is 80%. Thus, in the UBCSAND model, the top and middle areas of the cargo are likely to liquefy.
3.3.2 Heave motion
Considering the effect of acceleration on cargo fluidization, selecting the middle monitoring point as a representative, the variation of the excess pore pressure ratio and its residual value are shown in Figure 20. The excess pore pressure ratio rises wavelike. The larger of the acceleration, the higher of the excess pore pressure. For different parts of the cargo, the residual excess pore pressure ratio increases as the acceleration increases. In addition, with the increase of acceleration, the increment of residual pore pressure of the top area decreases. Because of that the pore pressure of the top area would dissipate under the influence of larger deformation at high acceleration. This phenomenon also exists in the FINN model.
Secondly, considering the effect of frequency, selecting the middle monitoring point as a representative, the variation of the excess pore pressure ratio and its residual value are shown in Figure 21. The excess pore pressure ratio first rises and then reaches a stable state. In addition, the critical values for the influence of frequency for cargo fluidization are corresponding to 3 Hz and 8 Hz, respectively. This is the same as the conclusion in the FINN model, and the two peak points in the FINN model both appear at 3 Hz and 8 Hz.
Finally, considering the effect of relative density, selecting the middle monitoring point as a representative, the variation of the excess pore pressure ratio and its residual value are shown in Figure 22. Under different relative densities, the excess pore pressure ratio increases firstly and then stabilizes. The corresponding excess pore pressure ratio is higher when the relative density is lower. In addition, as the relative density of the cargo increases, the ratio of the residual excess pore pressures in each area decreases, and the possibility of cargo fluidization becomes smaller. This is similar to the conclusion in the FINN constitutive model, but the difference is that the residual pore pressure values in FINN are all negative.
3.3.3 Roll motion
For different acceleration, taking the monitoring point at the middle as an example, the variation of the excess pore pressure ratio and its maximum value are shown in Figure 23. Under high acceleration, the excess pore pressure ratio continues to rise and stabilize to 1.0. And the time for the excess pore pressure ratio reaches stable at higher acceleration is shorter. Furthermore, the top and middle parts of the cargo both are likely to liquefy, and their maximum excess pore pressure ratios are close to 1.1. For the bottom part of the cargo, the larger of the acceleration value, the higher of the maximum excess pore pressure ratio. This is similar to the conclusion in the FINN constitutive model, and both of them can reach the maximum pore pressure near 1.0 under low acceleration.
For different frequencies, the variation of the excess pore pressure ratio at middle and its maximum value at different frequencies are shown in Figure 24. Under the different frequencies, the excess pore pressure ratio of the top and middle area can still rise to around 1.0. Differently, the excess pore pressure ratio of the bottom area is affected by the frequency obviously. And there is a critical value for the cargo fluidization, corresponding to 3 Hz. This is different from the conclusion in the FINN constitutive model. In the UBCSAND model, the maximum pore water pressure at the bottom has two trends with the increase of frequency, increasing first and then decreasing, while in the FINN model, there is only an increasing trend.
For different relative densities, taking the monitoring point at the middle as an example, the variation of the excess pore pressure ratio and its maximum value are shown in Figure 25. At lower relative density, the excess pore pressure ratio is higher. As the relative density increases, the maximum excess pore pressure ratio of the whole cargo decreases, and the possibility of cargo fluidization becomes smaller. This is the same as the conclusion in the FINN constitutive model.
3.3.4 Analysis
Based on the above results, acceleration and relative density have the same effect on the possibility of cargo fluidization under different motion forms. With the increase of acceleration, the possibility of cargo fluidization increases, as shown in Figures 17, 20 and 23. With the increase of relative density, the possibility of cargo fluidization reduces, as shown in Figures 19, 22, and 25. For frequency, its effect is different at different motions. For sway motion, the increase of frequency would reduce the possibility of cargo fluidization, as shown in Figure 18. For heave and roll motion, it exists a critical frequency on the possibility of cargo fluidization, as shown in Figures 21 and 24. Due to different liquefaction constitutive models, the specific influences of acceleration, frequency and relative density on cargo fluidization are different. But the influence trends of each factor on cargo fluidization are basically the same.
To further study the influence of different ship motions on cargo fluidization, simulations under same conditions were carried out. The curves of excess pore pressure ratio and its residual value are shown in Figure 26. The trend of the excess pore pressure ratio at each part under different motion forms is similar. And the possibility of cargo fluidization is the least under heave motion. In addition, the residual excess pore pressure ratio in the two rolling motions is higher than that in the swaying motion. And as show in Figure 23, the whole cargo is more likely to fluidize under rolling motion. Therefore, the rolling motion is the most critical ship motion that affects cargo fluidization.
3.4 Comparison of FINN and UBCSAND
As a FLAC embedded constitutive model, the FINN model is a simple liquefaction constitutive, which is used to solve liquefaction problems. Likewise, UBCSAND is also widely used in numerical simulation of liquefaction. It is a comprehensive liquefaction constitutive model, which has a fully coupled effective stress procedure (Itasca, 2011).
In view of computational efficiency and convenience, the FINN constitutive model is more suitable for the present problem. The parameters of FINN model are directly inherited for the static analysis. However, UBCSAND parameters need to be obtained additionally. And although the parameters of UBCSAND can be acquired by the relative density of the cargo from the empirical formula, its value is not accurate and calibration experiments are required when necessary.
For the accuracy of results, UBCSAND constitutive model may have more advantages. The FINN model is able to simulate changes in pore water pressure during material liquefaction and can better capture the response brought by changes of the above factors in prediction simulations. However, it does not capture the liquefied state at the bottom of the actual sample in the verification.
4 Conclusions
The prediction model of cargo fluidization considering ship motion forms, motion acceleration, frequency and relative density of cargo was established in this study. And two different liquefaction models were used and compared. Conclusions are summarized as follows:
• Generally, both the UBCSAND and the FINN constitutive model can simulate the cargo behavior better during cargo fluidization. The FINN model is relatively simple and convenient. By comparing the experiment results, UBCSAND model is more accurate for the present problems.
• Compared with sway and heave motion of the ship, the rolling motion is the key motion that influences the fluidization of the liquefiable cargo during sea transportation.
• The possibility of cargo fluidization is increased with the increase of acceleration. And it decreases with the increase of the relative density of cargo. However, the influence of frequency is not unique. In heave motion, there exists a critical frequency for cargo fluidization under the two constitutive models. In roll motion, the frequency to trigger the fluidization of the whole cargo in the UBCSAND model is smaller than that of FINN. In sway motion, for UBCSAND, the possibility of fluidization continuously reduced with the increase of the frequency. Different from that, there exists a critical value in the FINN.

Table 1 Material properties
Material Parameter Value Sand Density, ρ (kg/m^{3}) 1 600 Elasticity modulus, E (Pa) 2.938e6 Poisson's ratio, v 0.3 Angle of friction, φ (°) 32 Porosity, n 0.4 Permeability, k (m^{2}/(Pa^{*}s)) 1e8 Relative density, D_{r} 50% Density(kg/m^{3}), ρ 2 000 Model box Bulk modulus, K (Pa) 2e8 Shear modulus, G (Pa) 1e8 Water Density, ρ (kg/m^{3}) 1 000 Bulk modulus, K (Pa) 2e8 Table 2 Liquefaction constitutive model parameters
Model Parameter Value UBCSAND Relative density, D_{r} (%) 50 N_{1, 60} SPT 11 Elastic shear modulus number, K_{G}^{e} 964.4 Stress exponent for elastic bulk modulus, m_{e} 0.5 Elastic bulk modulus number, K_{B} 675.1 Stress exponent for elastic shear modulus, n_{e} 0.5 Plastic shear modulus number, K_{G}^{P} 450.1 Stress exponent for plastic shear modulus, n_{P} 0.4 Critical state friction angle, Φ_{cs}(°) 33 Peak friction angle, ϕ_{peak}(°) 35.2 Failure ratio, R_{f} 0.89 Model parameter, m_{}hfac1 1.0 Model parameter, m_{}hfac2 1.0 Relative density, D_{r} (%) 50 ρ Density(kg/m3), 1 600 Elasticity modulus, E (Pa) 2.938e6 Poisson's ratio, v 0.3 Angle of friction, φ (°) 32 FINN porosity, n 0.4 Permeability, k (m^{2}/(Pa^{*}s)) 1e8 Model constant, C_{1} 0.429 9 Model constant, C_{2} 0.930 4 Table 3 UBCSAND parameter calculation formula
Parameter Calculation formula N_{1, 60} (D_{r}/15)^{2} K_{G}^{e} K_{G}^{e} = 21.7 × 20 ×(N_{1, 60})^{0.333} K_{G}^{P} K_{G}^{P} = K_{G}^{e}×(N_{1, 60})^{2} × 0.003 + 100 n_{P} 0.4 m_{e} 0.5 Φ_{cs}(°) 33 K_{B} K_{B} = 0.7 × K_{G}^{e} ϕ_{peak}(°) ϕ_{Cs} + N_{1, 60}/5 n_{e} 0.5 R_{f} 1 − N_{1, 60}/100 Table 4 Dynamic model parameters
Cargo Cabin Parameter Value Parameter Value D_{r} (%) 65 80 90 ρ (kg/m^{3}) 7 900 ρ (kg/m^{3}) 1 650 1 700 1 720 K (Pa) 2e11 E (Pa) 1.7e7 1.9e7 2e7 G (Pa) 7e10 v 0.3 0.3 0.3 φ (°) 36.9 38.9 40.5 n 0.585 0.542 0.515 k (m^{2}/(pa^{*}s)) 1e8 1e8 1e8 Table 5 Parameters of liquefaction constitutive model
FINN UBCSAND Parameter Value Parameter Value D_{r} (%) 65 80 90 D_{r} (%) 65 80 90 C_{1} 0.223 1 0.132 8 0.098 9 1, 60 19 29 37 C_{2} 1.792 8 3.012 1 4.044 5 Table 6 Numerical simulation cases
Model Motion Acceleration (g) Frequency (Hz) Density (%) FINN sway 0.1  0.5 5 65 0.1 1  5 65 0.1 5 65, 80, 90 heave 0.1  0.6 1 65 0.1 1  9 65 0.1 1 65, 80, 90 roll 0.1  0.4 5 65 0.1 1  5 65 0.1 5 65, 80, 90 sway 0.1  0.15 1 65 0.1 1  5 65 0.1g 1 65, 80, 90 UBCSAND heave 0.1  0.4 1 65 0.1 1  9 65 0.1 1 65, 80, 90 roll 0.1  0.15 1 65 0.1 1  5 65 0.1 1 65, 80, 90 
Abbas M, Joshua D, Michael MC (2019) An Overview of the Behaviour of Iron Ore Fines Cargoes, and Some Recommended Solutions for the Reduction of Shifting Incidents During Marine Transportation. Ocean Engineering, 182: 451474. https://doi.org/10.1016/j.oceaneng.2019.04.073 Banerjee R, Konai S, Sengupta A, Deb K (2017) Shake Table Tests and Numerical Modeling of Liquefaction of Kasai River Sand. Geotechnical and Geological Engineering, 35(4): 13271340. https://doi.org/10.1007/s107060170178z Beaty MH (2018) Application of Ubcsand to the Leap Centrifuge Experiments. Soil Dynamics and Earthquake Engineering, 104: 143153. https://doi.org/10.1016/j.soildyn.2017.10.006 Beaty MH, Byrne PM (2011) Ubcsand Constitutive Model Version 904ar. Available from https://www.itasca.ca/software/udmlibrary/ubcsand. Chen WY, Wang ZH, Chen GX, Jeng DS, Wu M, Zhao HY (2018) Effect of vertical seismic motion on the dynamic response and instantaneous liquefaction in a twolayer porous seabed. Computers and Geotechnics, 99: 165176. https://doi.org/10.1016/j.compgeo.2018.03.005 Chou JC, Yang HT, Lin DG (2021) Calibration of Finn Model and UBCSAND Model for Simplified Liquefaction Analysis Procedures. Applied Sciences, 11(11): 52835283. https://doi.org/10.3390/APP11115283 Daoud S, Said I, Ennour S, Bouassida M (2018) Numerical Analysis of Cargo Liquefaction Mechanism Under the Swell Motion. Marine Structures, 57: 5271. https://doi.org/10.1016/j.marstruc.2017.09.003 Dashti S, Bray JD (2013) Numerical Simulation of Building Response on Liquefiable Sand. Journal of Geotechnical and Geoenvironmental Engineering, 139(8): 12351249. https://doi.org/10.1061/(ASCE)GT.19435606.0000853 Ding YH (2020) Study on Classification of A/C Group of Marine Solid Bulk Cargo. Master thesis, Dalian Maritime University, Dalian. (in Chinese) Ecemis N (2013) Simulation of Seismic Liquefaction: 1g Model Testing System and Shaking Table Tests. European Journal of Environmental and Civil Engineering, 17(10): 899919. https://doi.org/10.1080/19648189.2013.833140 Hao DZ (2019) Study on the Influence of Load Form on the Fluidization Process of Liquefiable Cargo. Master thesis, Zhejiang Ocean University, Zhoushan. (in Chinese) Itasca (2011) Online Manual of FLAC. Available from https://www.itascacg.com/software/flac. Jian QW, Li N, Zhou J, Li C (2015) Mesoscopic Mechanism of Fluidization for Bulk Iron Ore Concentrates Based on Model Test. Journal of Tongji University (Natural Science), 43(7): 10191024 (in Chinese) Ju L, Vassalos D, Boulougouris E (2016) Numerical Assessment of Cargo Liquefaction Potential. Ocean Engineering, 120: 383388. https://doi.org/10.1016/j.oceaneng.2016.01.024 Kutter BL, Manzari MT, Zeghal M (2020) Model Tests and Numerical Simulations of Liquefaction and Lateral Spreading. Springer, Cham, Switzerland. Munro MC, Mohajerani A (2016) Variation of the Geotechnical Properties of Iron Ore Fines Under Cyclic Loading. Ocean Engineering, 126: 411431. https://doi.org/10.1016/j.oceaneng.2016.09.006 Noui A, Karech T, Bouzid T (2019) A Numerical Investigation of Dynamic Behavior of a Unit Cell of a Loose Sand Reinforced By Stone Column Under the Effect of Gravity Using Finn Model. Indian Geotechnical Journal, 49(3): 255264. https://doi.org/10.1007/s4009801803262 Sakar C, Koseoglu B, Toz AC, Buber M (2020) Analysing the Effects of Liquefaction on Capsizing Through Integrating Interpretive Structural Modelling (ism) and Fuzzy Bayesian Networks (fbn). Ocean Engineering, 215. https://doi.org/10.1016/j.oceaneng.2020.107917 Vargehese RM, Iatha GM (2014) Shaking Table Tests to Investigate the Influence of Various Factors on the Liquefaction Resistance of Sands. Natural Hazards, 73(3): 13371351. https://doi.org/10.1007/s1106901411423 Vasiliki T, Stavroula K, Taborda DMG, Potts DM (2016) Vertical Ground Motion and Its Effects on Liquefaction Resistance of Fully Saturated Sand Deposits. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472(2192): 20160434. https://doi.org/10.1098/rspa.2016.0434 Wang H (2019) Research on Several Crucial Problems of Geotechnical Centrifuge Modeling Techniques. PhD thesis, Institute of Engineering Mechanics, China Earthquake Administration, Heilongjiang (in Chinese) Wang JT, Salam S, Xiao M (2020) Evaluation of the Effects of Shaking History on Liquefaction and Cone Penetration Resistance Using Shake Table Tests. Soil Dynamics and Earthquake Enginee ring, 131: 106025. https://doi.org/10.1016/j.soildyn.2019.106025 Wang MY, Zhao YT, Qian QH (2002) Study on Dynamic Behavior and Numerical Method for Saturated Sand. Chinese Journal of Geotechnical Engineering, 2002(6): 737742. (in Chinese) Wang YJ, Xing JB (1995) Discrete Element Method and Lagrangian Element Method and Their Applications in Geomechanics. Rock and Soil Mechanics, 1995(2): 114. https://doi.org/10.16285/j.rsm.1995.02.001(inChinese) Wu WQ, Ding YH, Xu HD, Zhang CC, Zhang QG, Zhang CL, Han JS, Chen YC (2020) Quantitative Safety Research of Water Transport of Liquefiable Solid Bulk Cargo Based on Wave Load. Journal of Dalian Maritime University, 46(02): 19. https://doi.org/10.16411/j.cnki.issn10067736.2020.02.001 (inChinese) Zhang JW, Wu WQ, Hu JQ (2017) Parametric Studies on Nickel Ore Slurry Sloshing in a Cargo Hold By Numerical Simulations. Ships and Offshore Structures, 12(2): 209218. https://doi.org/10.1080/17445302.2015.1131004 Zhang JW, Wu WQ, Zhao ZH, Chen YL (2019) Numerical Study on Coupled Effect of a Vessel Loaded with Liquefied Nickel Ore. Journal of Marine Science and Technology, 25: 116. https://doi.org/10.1007/s00773019006589 Zhou J, Du Q, Jian QW, Li C, Zhang ZQ (2016) Discrete Element Numerical Simulation of Model Experiments of Iron Concentrate Ore Fluidization During Shipping. Journal of Hohai University (Natural Sciences), 44(3): 219225 (in Chinese) Zhou J, Jian QW, Wu XH, Li N, Zhu YM (2013) Model Experimental Study of Fluidization of Iron Concentrate Ore In Bulk. Chinese Journal of Rock Mechanics and Engineering, 32 (12): 25362543 (in Chinese)