# Structural Optimization Design of Ship Lock Heads on Soft Soil considering Time-Varying Effects of the Structure and Foundation

Authors: Jiawei Bai; Chao Su; Heng Zhang; Shaopei Hu
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Over time, the uneven settlements of the structure and foundation are prominent in constructing ship lock heads on soft soil. These deformations endanger the safety of ship lock heads during construction. This research aimed to establish the ship lock head’s structural optimization procedure on soft soil, considering the time-varying effects of the structure and foundation. By comprehensively considering the linear viscoelastic creep of concrete and the elastoplastic consolidation characteristic of soft soil, a perfect time-dependent analysis method for the lock head on soft soil was proposed. Furthermore, a hybrid particle swarm optimization, enhanced whale optimization, and differential evolution (PSO-EWOA-DE) algorithm was proposed to optimize thirty-four design variables of a lock head. With the minimal volume of the lock head as the optimization objective, the finite element model was established. In the optimization process, three types of constraints were evaluated. The result showed that the optimized design could reduce 10.45% of structure volume. Through comparing and analysing the maximum principle stresses and vertical displacements of the lock head before and after optimization, some conclusions were drawn. The optimization procedure proposed in this paper provides a new perspective for the structural optimization of hydraulic structures on soft soil.


---

## Body

## 1. Introduction

As one of the most critical navigation structures, ship locks use hydraulic power to lift ships through dams built on natural or canalized rivers [1]. Ship locks are mostly built on soft-soil sites in coastal areas. Soft soil has the characteristics of low shear strength [2], high compressibility, and a large settlement. In practical engineering, building crashes caused by uneven foundation settlement are not uncommon. Therefore, the accurate calculation of foundation settlement is an issue of particular concern to designers and contractors. The lock head is a complicated structure with water-retaining and ship passing functions, and it separates the lock chamber from the approach channel. In the past few years, some progress in the lock head’s structural optimization has been achieved. Based on the material’s elastic assumption, Su et al. [3] optimized the lock head volume using the finite element method (FEM) to effectively use the concrete material. Subsequently, considering the impact of concrete creep on the structure during the construction period, Su and Bai [4] optimized the lock head on the rock foundation. Under the lock head’s gravity, the soft foundation will experience consolidation settlement. Therefore, the structural optimization of the lock head should consider the time-varying characteristics of the structure and foundation.Many algorithms have been used in structural engineering optimization tasks over the past few decades, from gradient-based algorithms to nongradient probabilistic-based algorithms [5]. The latter category contains many algorithms that imitate natural phenomena, such as ant colony algorithm [6] and particle swarm optimization (PSO) [7], among others. The gradient-based algorithms have difficulties in dealing with inequality constraints [8]. Evolutionary algorithms (EAs) are the most successful optimization techniques [9]. The PSO is a variant of EA based on bird flocking [10]. It was proposed based on birds’ cardinal rules: the simultaneous gathering of a large number of birds, often changing direction suddenly, dispersing, and rearranging [11]. The PSO has few parameters and a fast convergence rate, and it is easy to implement [12]. Nowadays, the PSO has achieved encouraging results in structural design applications. Fourie and Groenwold solved the shape [11] and topology [13] optimization problems using the PSO. Venter and Sobieszczanski-Sobieski [14] applied the PSO to the multidisciplinary optimization of transport aircraft wings. The whale optimization algorithm (WOA) is a novel population-based metaheuristic optimization algorithm, originally developed by Mirjalili and Lewis [15]. Some important parameters of the WOA are self-tuning along the evolutionary process, and the WOA also needs few parameters to set [16]. It mimics the hunting behaviour of humpback whales [17]. The main hunting steps are included in this method: searching for the prey, encircling the prey, and attacking the prey [18]. Compared with other metaheuristic algorithms, the whale algorithm can obtain very competitive results in solving various practical problems [19] such as determining the optimal distributed generator size [20], selecting features [21], and clustering the data [22]. Kaveh and Ilchi Ghazaan [23] proposed the enhanced whale optimization algorithm (EWOA) to improve the solution’s accuracy. Compared with WOA, EWOA has fewer parameters and conditional statements [24]. The difference evolution (DE) algorithm is a robust population-based global optimization algorithm [25]. It performs effectively in solving constrained parameter optimization problems [8]. The main DE steps include mutation, crossover, and selection [26]. The DE has a reliable search capability [27]. Recently, to improve optimization results, some optimization methods that hybridize multiple evolutionary techniques have been proposed [26]. Zhou et al. applied the PSO-DE algorithm to optimize the operation of heated oil pipelines [28].The creep phenomenon was discovered by Hatt in 1907 [29]. For early-age concrete, thermal stresses are greatly affected by creep deformations [30]. Concrete creep is the strain that changes with time caused by continuous stress [31]. It will influence the serviceability of the concrete structure [32]. Many experiments have been performed to collect short-term and long-term creep data of concrete [33]. Subsequently, a large number of creep prediction methods were proposed based on experimental results [34], such as the GL2000 model [35], the BP model [36], and the ACI 209R-92 model [37]. Bažant and Xi [38] believed that the Kelvin chain could accurately describe the linear viscoelastic creep model, and they used the Dirichlet series to represent the creep compliance function. Zhu [39] proposed a series of implicit recursive equations and provided a large number of calculation parameters suitable for hydraulic concrete structures. Wang and Liu [40] verified this method’s correctness by comparing the calculation results of the Gongzui gravity dam’s thermal creep stress with theoretical solutions.Consolidation is the time-settlement behaviour of clay. The clay instantaneously settles under the load’s action without drainage, and then the excess pore pressure gradually dissipates [41]. Terzaghi first proposed the one-dimensional consolidation theory in 1925 [42]. Biot further studied the coupling effect of three-dimensional deformation materials and pore fluid pressure and developed a more general consolidation theory [43]. Sandu and Wilson presented a solution for Biot’s consolidation theory based on the FEM [44]. Since local yielding occurs when the soil is loaded, Small et al. formulated the elastoplastic consolidation equations on the basis of incremental solution strategies [41]. The modified Cam-clay model (MCC) is a critical-state model most widely used to predict soft soils’ behaviour [45]. It only requires a few parameters which are convenient for laboratory determination [46]. For soft soils, Jain deduced the elastoplastic consolidation equations incorporating MCC [47]. However, few optimization studies have considered the consolidation characteristics of soft foundations.In this article, a structural optimization method for ship lock heads on soft soil was established considering the time-varying effects of the structure and foundation. Meanwhile, by comprehensively considering the linear viscoelastic creep of concrete and the elastoplastic consolidation characteristic of soft soil, a perfect time-dependent analysis method for the lock head on soft soil was proposed. With the minimal volume of the lock head as the optimization objective, a hybrid particle swarm optimization, enhanced whale optimization, and differential evolution (PSO-EWOA-DE) algorithm was applied to obtain the optimal shape. A comparison of the results obtained by the PSO-EWOA-DE algorithm with those of PSO, WOA, DE, EWOA, and PSO-DE methods revealed that the proposed algorithm was more effective in the structural optimization of a ship lock head and had a faster convergence rate. The remainder of this paper is structured as follows. Section2 deliberates governing equations of concrete thermal creep stress and soil consolidation theory. Section 3 describes the structural optimization problem in detail. The PSO, WOA, EWOA, and DE algorithms are briefly introduced in Section 4. Section 5 introduces the PSO-EWOA-DE algorithm and its main implementation details in structural optimization. Section 6 presents an engineering example of a ship lock head. Finally, some conclusions are given in Section 7.


## 2. Governing Equations

### 2.1. Thermal Creep Stress of Mass Concrete

For mass concrete structures, the compliance function is expressed as follows [48]:(1)Jt,t′=1Eτ+Ct,t′.The creep compliance function can be written in the following Dirichlet series form [39]:(2)Ct,t′=∑i=1Nψit′1−exp−αit−t′.In addition,(3)Et′=E01−exprt′λ,(4)ψit′=fi+git′pi,where t is the age of concrete, t′ is the loading age in days, N is the number of Kelvin chains, Et′ is the instantaneous elastic modulus, E0 is the ultimate elastic modulus, and αi, r, λ, fi, gi, and pi are parameters obtained by creep experiments.For viscoelastic materials, stress increment in 3D can be calculated by the following equation [39]:(5)Δσn=DnΔεn−ηn−ΔεnT,where(6)Dn=1Jt,t′Q−1,(7)ηn=∑i=1N1−exp−rΔtn′ωin,(8)ΔεnT=γΔT111000T.In these equations,Δσn is the stress increment, Δεn is the total strain increment, ΔεnT is the thermal strain increment, ωin is the implicit integral operator, Dn is the viscoelastic modulus, Q is Poisson’s ratio matrix, γ is the linear expansion coefficient, and ΔT is the temperature difference [4].


### 2.2. Finite Element Formulations for Soil Consolidation Theory

#### 2.2.1. MCC Model

According to the classical elastoplastic theory of geotechnical engineering, the stress-strain relationship is shown as follows [49]:(9)dσ=Ce+dε−dεp.In addition,(10)Ce=3K1−μ1+μ1μ1−μμ1−μ000μ1−μ1μ1−μ000μ1−μμ1−μ10000001−2μ21−μ0000001−2μ21−μ0000001−2μ21−μ,where dσ represents the stress increment, dε represents the total strain increment, dεp represents the plastic strain increment, Ce represents the elastic-stiffness matrix, K represents the volumetric modulus, and μ represents Poisson’s ratio.For an isotropic hardening material, the yield function is given as follows [49]:(11)fσij,A=0,where A is the hardening parameter for saturated soil. Using the chain rule for the differential, equation (11) can be converted as follows:(12)df=∂f∂σTdσ+∂f∂A∂f∂εpTdεp.According to the flow rule, the plastic strain increment is defined as the following equation [50]:(13)dεp=Λ′∂g∂σ,where(14)Λ′=∂f/∂σTCedε−∂f/∂A∂A/∂εpT∂g/∂σ+∂f/∂σTCe∂g/∂σ,where g is the plastic potential function and Λ′ is the plastic multiplier.In the MCC model, the volumetric modulusK is related to the stress history [51], and the hardening parameter A is the volumetric plastic strain [46]:(15)K=1+e0κ¯p′,(16)A=εvp.In these equations,e0 represents the initial void ratio, κ¯ represents the unloading-reloading index of soft soils, and p′ represents the mean effective stress.The yield function of the MCC model is given as follows [49]:(17)f=λ¯−κ¯1+e0lnp′p0+λ¯−κ¯1+e0ln1+q2M2p′2−εvp=0,where λ¯ represents the virgin compression index, p0 represents the preconsolidation pressure, q represents the deviatoric stress, and M represents the critical-state parameter.In the Cartesian coordinate system, the mean effective stressp′ and the deviatoric stress q can be obtained by the following equations, respectively [52].(18)p′=13σx′+σy′+σz′,(19)q=σx′2+σy′2+σz′2+3τxy2+τyz2+τzx2−σx′σy′+σy′σz′+σz′σx′.According to the consistency condition, the derivatives of the yield functionf and plastic potential function g with respect to stress σ are [49](20)∂f∂σ=∂g∂σ=λ¯−κ¯1+e0p′M2−η2M2+η2+2ηM2+η2,where(21)η=qp′.Based on equations (13) and (14), equation (9) can be converted into the following equation:(22)dσ=Ce−Ce∂g/∂σ∂f/∂σTCe−∂f/∂A∂A/∂εpT∂g/∂σ+∂f/∂σTCe∂g/∂σdε,where(23)∂f∂A=−1,(24)∂A∂εp=∂A∂εvp∂A∂εdp=10.According to equation (22), the relationship between stress and strain increment of the MCC model is described as follows.The thermal strain increment is calculated by equation (22):(25)dσ=Cepdε.


#### 2.2.2. Elastoplastic Consolidation

For the saturated continuous soil, the equilibrium equation can be discretized into the following form [47]:(26)Kδ=R¯,where K represents the stiffness matrix, δ represents the displacement vector at nodal points, and R¯ represents the load vector. The fluid flow equation is given as [47](27)Hq¯=R¯p,where H represents the fluid flow matrix, q¯ represents the pore water pressure vector at nodal points, and R¯p represents the applied flux.Based on the theory of seepage, the pore water pressurepw exerted on the soil skeleton can be expressed as follows [43]:(28)Xi=−∂pw∂Xi.Equation (28) can be discretized into the following form using the FEM:(29)Xi=−∂∂XiN¯pq¯,where N¯p represents the shape function for pore water pressure. According to Biot’s theory, the nodal force Fp is defined as the following equation [53]:(30)Fp=∫VN¯TXidV=−LTq¯.In addition,(31)LT=∫VN¯T∂∂XiN¯pdV,where N¯ represents the displacement shape function and L represents the coupling matrix. Therefore, equation (26) can be converted into the following equation:(32)Kδ−LTq¯=R¯.For a porous elastic medium, the continuity equation can be expressed as follows:(33)Q=εii˙=−∂∂t∂u∂x+∂v∂y+∂w∂z,where εii˙ represents the rate of volume change. The discrete form of equations (2)–(33) is given in the following equation:(34)Q=−∂∂t∂xiTN¯δ.Considering the fluid-skeleton compatibility, the fluid flow equation is converted into the following equation [47]:(35)Hq¯−L∂∂tδ=R¯p,where(36)L∂∂tδ=−∫VN¯pTQdV.In a time incrementΔt (t1⟶t2), equations (2)–(35) can be solved using the finite difference approximation [41].(37)Hq¯1+βq¯2−q¯1Δt−Lδ2−δ1=0,where β corresponds to the particular integration rule. The assembled form of equations (32) and (37) is shown as follows [41]:(38)K−LT−LβΔtHδ2q¯2=R¯2−Lδ1−1−βΔtHq¯1.In addition,(39)K=∫VBTCeBdV,where B represents the strain rate-velocity matrix. For the elastoplastic MCC model, Ce can be suited by Cep [54].


## 2.1. Thermal Creep Stress of Mass Concrete

For mass concrete structures, the compliance function is expressed as follows [48]:(1)Jt,t′=1Eτ+Ct,t′.The creep compliance function can be written in the following Dirichlet series form [39]:(2)Ct,t′=∑i=1Nψit′1−exp−αit−t′.In addition,(3)Et′=E01−exprt′λ,(4)ψit′=fi+git′pi,where t is the age of concrete, t′ is the loading age in days, N is the number of Kelvin chains, Et′ is the instantaneous elastic modulus, E0 is the ultimate elastic modulus, and αi, r, λ, fi, gi, and pi are parameters obtained by creep experiments.For viscoelastic materials, stress increment in 3D can be calculated by the following equation [39]:(5)Δσn=DnΔεn−ηn−ΔεnT,where(6)Dn=1Jt,t′Q−1,(7)ηn=∑i=1N1−exp−rΔtn′ωin,(8)ΔεnT=γΔT111000T.In these equations,Δσn is the stress increment, Δεn is the total strain increment, ΔεnT is the thermal strain increment, ωin is the implicit integral operator, Dn is the viscoelastic modulus, Q is Poisson’s ratio matrix, γ is the linear expansion coefficient, and ΔT is the temperature difference [4].


## 2.2. Finite Element Formulations for Soil Consolidation Theory

### 2.2.1. MCC Model

According to the classical elastoplastic theory of geotechnical engineering, the stress-strain relationship is shown as follows [49]:(9)dσ=Ce+dε−dεp.In addition,(10)Ce=3K1−μ1+μ1μ1−μμ1−μ000μ1−μ1μ1−μ000μ1−μμ1−μ10000001−2μ21−μ0000001−2μ21−μ0000001−2μ21−μ,where dσ represents the stress increment, dε represents the total strain increment, dεp represents the plastic strain increment, Ce represents the elastic-stiffness matrix, K represents the volumetric modulus, and μ represents Poisson’s ratio.For an isotropic hardening material, the yield function is given as follows [49]:(11)fσij,A=0,where A is the hardening parameter for saturated soil. Using the chain rule for the differential, equation (11) can be converted as follows:(12)df=∂f∂σTdσ+∂f∂A∂f∂εpTdεp.According to the flow rule, the plastic strain increment is defined as the following equation [50]:(13)dεp=Λ′∂g∂σ,where(14)Λ′=∂f/∂σTCedε−∂f/∂A∂A/∂εpT∂g/∂σ+∂f/∂σTCe∂g/∂σ,where g is the plastic potential function and Λ′ is the plastic multiplier.In the MCC model, the volumetric modulusK is related to the stress history [51], and the hardening parameter A is the volumetric plastic strain [46]:(15)K=1+e0κ¯p′,(16)A=εvp.In these equations,e0 represents the initial void ratio, κ¯ represents the unloading-reloading index of soft soils, and p′ represents the mean effective stress.The yield function of the MCC model is given as follows [49]:(17)f=λ¯−κ¯1+e0lnp′p0+λ¯−κ¯1+e0ln1+q2M2p′2−εvp=0,where λ¯ represents the virgin compression index, p0 represents the preconsolidation pressure, q represents the deviatoric stress, and M represents the critical-state parameter.In the Cartesian coordinate system, the mean effective stressp′ and the deviatoric stress q can be obtained by the following equations, respectively [52].(18)p′=13σx′+σy′+σz′,(19)q=σx′2+σy′2+σz′2+3τxy2+τyz2+τzx2−σx′σy′+σy′σz′+σz′σx′.According to the consistency condition, the derivatives of the yield functionf and plastic potential function g with respect to stress σ are [49](20)∂f∂σ=∂g∂σ=λ¯−κ¯1+e0p′M2−η2M2+η2+2ηM2+η2,where(21)η=qp′.Based on equations (13) and (14), equation (9) can be converted into the following equation:(22)dσ=Ce−Ce∂g/∂σ∂f/∂σTCe−∂f/∂A∂A/∂εpT∂g/∂σ+∂f/∂σTCe∂g/∂σdε,where(23)∂f∂A=−1,(24)∂A∂εp=∂A∂εvp∂A∂εdp=10.According to equation (22), the relationship between stress and strain increment of the MCC model is described as follows.The thermal strain increment is calculated by equation (22):(25)dσ=Cepdε.


### 2.2.2. Elastoplastic Consolidation

For the saturated continuous soil, the equilibrium equation can be discretized into the following form [47]:(26)Kδ=R¯,where K represents the stiffness matrix, δ represents the displacement vector at nodal points, and R¯ represents the load vector. The fluid flow equation is given as [47](27)Hq¯=R¯p,where H represents the fluid flow matrix, q¯ represents the pore water pressure vector at nodal points, and R¯p represents the applied flux.Based on the theory of seepage, the pore water pressurepw exerted on the soil skeleton can be expressed as follows [43]:(28)Xi=−∂pw∂Xi.Equation (28) can be discretized into the following form using the FEM:(29)Xi=−∂∂XiN¯pq¯,where N¯p represents the shape function for pore water pressure. According to Biot’s theory, the nodal force Fp is defined as the following equation [53]:(30)Fp=∫VN¯TXidV=−LTq¯.In addition,(31)LT=∫VN¯T∂∂XiN¯pdV,where N¯ represents the displacement shape function and L represents the coupling matrix. Therefore, equation (26) can be converted into the following equation:(32)Kδ−LTq¯=R¯.For a porous elastic medium, the continuity equation can be expressed as follows:(33)Q=εii˙=−∂∂t∂u∂x+∂v∂y+∂w∂z,where εii˙ represents the rate of volume change. The discrete form of equations (2)–(33) is given in the following equation:(34)Q=−∂∂t∂xiTN¯δ.Considering the fluid-skeleton compatibility, the fluid flow equation is converted into the following equation [47]:(35)Hq¯−L∂∂tδ=R¯p,where(36)L∂∂tδ=−∫VN¯pTQdV.In a time incrementΔt (t1⟶t2), equations (2)–(35) can be solved using the finite difference approximation [41].(37)Hq¯1+βq¯2−q¯1Δt−Lδ2−δ1=0,where β corresponds to the particular integration rule. The assembled form of equations (32) and (37) is shown as follows [41]:(38)K−LT−LβΔtHδ2q¯2=R¯2−Lδ1−1−βΔtHq¯1.In addition,(39)K=∫VBTCeBdV,where B represents the strain rate-velocity matrix. For the elastoplastic MCC model, Ce can be suited by Cep [54].


## 2.2.1. MCC Model

According to the classical elastoplastic theory of geotechnical engineering, the stress-strain relationship is shown as follows [49]:(9)dσ=Ce+dε−dεp.In addition,(10)Ce=3K1−μ1+μ1μ1−μμ1−μ000μ1−μ1μ1−μ000μ1−μμ1−μ10000001−2μ21−μ0000001−2μ21−μ0000001−2μ21−μ,where dσ represents the stress increment, dε represents the total strain increment, dεp represents the plastic strain increment, Ce represents the elastic-stiffness matrix, K represents the volumetric modulus, and μ represents Poisson’s ratio.For an isotropic hardening material, the yield function is given as follows [49]:(11)fσij,A=0,where A is the hardening parameter for saturated soil. Using the chain rule for the differential, equation (11) can be converted as follows:(12)df=∂f∂σTdσ+∂f∂A∂f∂εpTdεp.According to the flow rule, the plastic strain increment is defined as the following equation [50]:(13)dεp=Λ′∂g∂σ,where(14)Λ′=∂f/∂σTCedε−∂f/∂A∂A/∂εpT∂g/∂σ+∂f/∂σTCe∂g/∂σ,where g is the plastic potential function and Λ′ is the plastic multiplier.In the MCC model, the volumetric modulusK is related to the stress history [51], and the hardening parameter A is the volumetric plastic strain [46]:(15)K=1+e0κ¯p′,(16)A=εvp.In these equations,e0 represents the initial void ratio, κ¯ represents the unloading-reloading index of soft soils, and p′ represents the mean effective stress.The yield function of the MCC model is given as follows [49]:(17)f=λ¯−κ¯1+e0lnp′p0+λ¯−κ¯1+e0ln1+q2M2p′2−εvp=0,where λ¯ represents the virgin compression index, p0 represents the preconsolidation pressure, q represents the deviatoric stress, and M represents the critical-state parameter.In the Cartesian coordinate system, the mean effective stressp′ and the deviatoric stress q can be obtained by the following equations, respectively [52].(18)p′=13σx′+σy′+σz′,(19)q=σx′2+σy′2+σz′2+3τxy2+τyz2+τzx2−σx′σy′+σy′σz′+σz′σx′.According to the consistency condition, the derivatives of the yield functionf and plastic potential function g with respect to stress σ are [49](20)∂f∂σ=∂g∂σ=λ¯−κ¯1+e0p′M2−η2M2+η2+2ηM2+η2,where(21)η=qp′.Based on equations (13) and (14), equation (9) can be converted into the following equation:(22)dσ=Ce−Ce∂g/∂σ∂f/∂σTCe−∂f/∂A∂A/∂εpT∂g/∂σ+∂f/∂σTCe∂g/∂σdε,where(23)∂f∂A=−1,(24)∂A∂εp=∂A∂εvp∂A∂εdp=10.According to equation (22), the relationship between stress and strain increment of the MCC model is described as follows.The thermal strain increment is calculated by equation (22):(25)dσ=Cepdε.


## 2.2.2. Elastoplastic Consolidation

For the saturated continuous soil, the equilibrium equation can be discretized into the following form [47]:(26)Kδ=R¯,where K represents the stiffness matrix, δ represents the displacement vector at nodal points, and R¯ represents the load vector. The fluid flow equation is given as [47](27)Hq¯=R¯p,where H represents the fluid flow matrix, q¯ represents the pore water pressure vector at nodal points, and R¯p represents the applied flux.Based on the theory of seepage, the pore water pressurepw exerted on the soil skeleton can be expressed as follows [43]:(28)Xi=−∂pw∂Xi.Equation (28) can be discretized into the following form using the FEM:(29)Xi=−∂∂XiN¯pq¯,where N¯p represents the shape function for pore water pressure. According to Biot’s theory, the nodal force Fp is defined as the following equation [53]:(30)Fp=∫VN¯TXidV=−LTq¯.In addition,(31)LT=∫VN¯T∂∂XiN¯pdV,where N¯ represents the displacement shape function and L represents the coupling matrix. Therefore, equation (26) can be converted into the following equation:(32)Kδ−LTq¯=R¯.For a porous elastic medium, the continuity equation can be expressed as follows:(33)Q=εii˙=−∂∂t∂u∂x+∂v∂y+∂w∂z,where εii˙ represents the rate of volume change. The discrete form of equations (2)–(33) is given in the following equation:(34)Q=−∂∂t∂xiTN¯δ.Considering the fluid-skeleton compatibility, the fluid flow equation is converted into the following equation [47]:(35)Hq¯−L∂∂tδ=R¯p,where(36)L∂∂tδ=−∫VN¯pTQdV.In a time incrementΔt (t1⟶t2), equations (2)–(35) can be solved using the finite difference approximation [41].(37)Hq¯1+βq¯2−q¯1Δt−Lδ2−δ1=0,where β corresponds to the particular integration rule. The assembled form of equations (32) and (37) is shown as follows [41]:(38)K−LT−LβΔtHδ2q¯2=R¯2−Lδ1−1−βΔtHq¯1.In addition,(39)K=∫VBTCeBdV,where B represents the strain rate-velocity matrix. For the elastoplastic MCC model, Ce can be suited by Cep [54].


## 3. The Optimization Problem Statement

### 3.1. Geometric Description of the Ship Lock Head

Parametric modelling is the basis for structural optimization design. The ship lock head is a symmetrical structure along the water flow direction. Therefore, half of the geometric model was established. As shown in Figure1(a), a typical lock head consists of the bottom plate, the corridor layer, the empty-box layer, and the second-stage concrete. The model of the lock head-soft foundation-backfilled soil system is shown in Figure 1(b).Figure 1
Schematics of the spatial (a) ship lock head and (b) lock head-soft foundation-backfilled soil system.
(a)(b)The horizontal and vertical views of the bottom plate are shown in Figure2. The bottom plate is a regular structure. It can be established by determining Xpi (i = 1–3), Ypi (i = 1–8), and Zpi (i = 1–6). As shown in Figure 3, the design of the second-stage concrete also needs to determine the central angles αqi (i = 1–3), the arc radii Rqi (i = 1–2), and the length parameters Lqi (i = 1–7) [4]. The corridor layer section is shown in Figure 4.Figure 2
The basic dimensions of the bottom plate. (a) Plane figure of the bottom plate. (b) Sectional drawing of the bottom plate.
(a)(b)Figure 3
The second-stage concrete section.Figure 4
The corridor layer section.


### 3.2. Design Variables

Designers save concrete materials by setting up empty boxes in the lock head. According to this method, 34 design variables (x1–x15, y1–y12, and z1–z7) were distinguished in the x-, y-, and z-directions [4], as displayed in Figure 5.Figure 5
Design variables. (a) Horizontal view. (b) Elevation view.
(a)(b)


### 3.3. Objective Function and Constraints

The purpose of structural optimization is to find the optimal design that meets the stability and financial requirements. In this work, the lock head volume was taken as the objective function. The optimization problem can be expressed as follows:(40)minimize:ΑX,subjected  to:φjX≤0,Ximin≤X≤Ximax,where X = (x1, …, x15, y1, …, y12, z1, …, z7) is the updating design variables, ΑX is the volume function, j is the number of constraints, Ximin is the lower boundary of X, and Ximax is the upper boundary of X.According to the previous research [4], three types of constraints should be satisfied, which are formulated as follows:(41)stress constraints:σ1≤σ1σ3≤σ3,geometric constraints:Ximin≤X≤Ximax,stability constraints:Ka≥KaKo≥Ko,where σ1 and σ3 represent the maximum principal tensile stress and the maximum principal stress, respectively. Moreover, σ1 and σ3 represent allowable values, respectively. Ka and Ko represent the safety factors of the antisliding stability and anti-overturning stability, respectively. Ka and Ko represent allowable values, respectively.


## 3.1. Geometric Description of the Ship Lock Head

Parametric modelling is the basis for structural optimization design. The ship lock head is a symmetrical structure along the water flow direction. Therefore, half of the geometric model was established. As shown in Figure1(a), a typical lock head consists of the bottom plate, the corridor layer, the empty-box layer, and the second-stage concrete. The model of the lock head-soft foundation-backfilled soil system is shown in Figure 1(b).Figure 1
Schematics of the spatial (a) ship lock head and (b) lock head-soft foundation-backfilled soil system.
(a)(b)The horizontal and vertical views of the bottom plate are shown in Figure2. The bottom plate is a regular structure. It can be established by determining Xpi (i = 1–3), Ypi (i = 1–8), and Zpi (i = 1–6). As shown in Figure 3, the design of the second-stage concrete also needs to determine the central angles αqi (i = 1–3), the arc radii Rqi (i = 1–2), and the length parameters Lqi (i = 1–7) [4]. The corridor layer section is shown in Figure 4.Figure 2
The basic dimensions of the bottom plate. (a) Plane figure of the bottom plate. (b) Sectional drawing of the bottom plate.
(a)(b)Figure 3
The second-stage concrete section.Figure 4
The corridor layer section.


## 3.2. Design Variables

Designers save concrete materials by setting up empty boxes in the lock head. According to this method, 34 design variables (x1–x15, y1–y12, and z1–z7) were distinguished in the x-, y-, and z-directions [4], as displayed in Figure 5.Figure 5
Design variables. (a) Horizontal view. (b) Elevation view.
(a)(b)


## 3.3. Objective Function and Constraints

The purpose of structural optimization is to find the optimal design that meets the stability and financial requirements. In this work, the lock head volume was taken as the objective function. The optimization problem can be expressed as follows:(40)minimize:ΑX,subjected  to:φjX≤0,Ximin≤X≤Ximax,where X = (x1, …, x15, y1, …, y12, z1, …, z7) is the updating design variables, ΑX is the volume function, j is the number of constraints, Ximin is the lower boundary of X, and Ximax is the upper boundary of X.According to the previous research [4], three types of constraints should be satisfied, which are formulated as follows:(41)stress constraints:σ1≤σ1σ3≤σ3,geometric constraints:Ximin≤X≤Ximax,stability constraints:Ka≥KaKo≥Ko,where σ1 and σ3 represent the maximum principal tensile stress and the maximum principal stress, respectively. Moreover, σ1 and σ3 represent allowable values, respectively. Ka and Ko represent the safety factors of the antisliding stability and anti-overturning stability, respectively. Ka and Ko represent allowable values, respectively.


## 4. Concepts of Particle Swarm Optimization (PSO), Whale Optimization Algorithm (WOA), Enhanced Whale Optimization Algorithm (EWOA), and Differential Evolution Algorithm (DE)

### 4.1. Concept of PSO

The PSO is a swarm intelligence-based algorithm proposed by Kennedy and Eberhart [10]. It was derived from observing the bird predation behaviour. PSO uses particles to represent some candidate solutions to the optimization problem in the search space. Position, velocity, and the memory mechanism are three characteristics of each particle [55]. Initially, the particles are randomly generated within the ranges of design variables. The PSO evaluates each particle’s fitness value by calculating the objective function. Then, the position of each design variable is updated by the following equations:(42)Xi,jk+1=Xi,jk+Vi,jk+1.In addition,(43)Vi,jk+1=wkVi,jk+C1⋅r1Pbi,jk−Xi,jk+C2⋅r2gbi,jk−Xi,jk,where Xi,jk is the current jth design variable of particle i at the kth iteration, Vi,jk is the corresponding velocity of the jth design variable, Pbi,jk is the best position at the kth iteration, gbi,jk is the global best position, r1 and r2 are random numbers between [0, 1], C1 and C2 are user-supplied coefficients, and wk is the inertia weight. In this study, the inertia weight is defined as follows:(44)wk=wmax−wmax−wminItmax⋅k,where wmax is the maximum value of the inertia weight and wmin is the minimum value of the inertia weight. Itmax is the maximum generation, which is the termination criterion of this study.


### 4.2. Concept of WOA

The WOA mimics the bubble-net hunting behaviour of humpback whales. The unique hunting model is shown in Figure6. The WOA consists of two phases: the exploitation phase and the exploration phase [21]. The mathematical models of the two stages are discussed in the following sections.Figure 6
Unique bubble-net spiral hunting behaviour.


#### 4.2.1. Exploitation Phase

As shown in Figure7, the shrinking encircling mechanism and bubble-net spiral attacking mechanism are included in the exploitation phase. In WOA, the global best position is assumed as the target prey of the current generation, and other candidates will update their positions in a shrinking encircling manner. The mathematical model of the shrinking encircling mechanism is shown as follows:(45)Xi,jk+1=gbi,jk−ADi,jk+1,(46)Di,jk+1=C⋅gbi,jk−Xi,jk,where A and C are coefficients of WOA and can be calculated as follows:(47)A=2ar−a,C=2r,a=2−k2Itmax,where r is a random number in the interval of [0, 1]. To simulate the bubble-net attacking behaviour of humpback whales, a spiral equation is created as follows:(48)Xi,jk+1=Di,j′k+1⋅ebl⋅cos2πl+gbi,jk.Figure 7
Exploitation phase in WOA. (a) Shrinking encircling mechanism. (b) Bubble-net spiral attacking mechanism.
(a)(b)In this equation,Di,j′k+1=gbi,jk−Xi,jk, l is a random number in the interval of [−1, 1], and constant b defines the spiral shape. The humpback whales can swim around their prey in the above two paths at the same time. The mathematical model of this process is shown as follows:(49)Xi,jk+1=gbi,jk−ADi,jk+1if prob<0.5,Di,j′k+1⋅ebl⋅cos2πl+gbi,jkif prob≥0.5.


#### 4.2.2. Exploration Phase

To expand the searching area, the WOA updates each design variable’s position using a randomly selected candidate whenA≥1. The mathematical model for the second phase is represented by equations (50) and (51):(50)Xi,jk+1=Xri,jk−ADi,j″k+1,(51)Di,j″k+1=C⋅Xri,jk−Xi,jk,where Xri,jk represents a randomly selected search agent in the current population. The flowchart of WOA is presented in Figure 8(a).Figure 8
Flowcharts of the (a) WOA algorithm and (b) EWOA algorithm.
(a)(b)


### 4.3. Concept of EWOA

The WOA has an effective global search capability. In order to enhance the WOA in terms of reliability, EWOA was proposed. The EWOA maintains the simplicity of the WOA. In the exploitation phase, to maintain a balance between the diversity and intensification of search results, equation (45) is replaced as follows [23]:(52)Xi,jk+1=gbi,jk−ADi,j″′k+1,(53)Di,j″′k+1=rXi,jk.In the exploration phase, each candidate randomly changes the value of a design variable with a probabilityp″:(54)p″=0.31−kItmax.For candidatei, the design variable is changed when a random number q is less than probability p″. The mathematical model of this process is shown as follows [23]:(55)Xi,jk+1=Xjmin+random⋅Xjmax−Xjmin,q<p″,Xi,jk,q≥p″,where random and q represent random numbers between [0, 1], respectively. Xjmin and Xjmax represent the minimum value and the maximum value of the jth design variable, respectively. The flowchart of EWOA is presented in Figure 8(b).


### 4.4. Concept of DE

The DE is a simple population-based evolutionary algorithm proposed by Storn and Price [25]. Similar to PSO, it starts from randomly generating an initial population within the ranges of design variables [56]. The mutation process is to create a mutated vector Yi,jk+1 according to equation (56). The crossover process is to create a trial vector Zi,jk+1, as described in equation (57) [25]. In the selection process, whether the trial vector should be selected for the next generation can be determined according to equation (58) [8].(56)Yi,jk+1=Xran1,jk+FXran2,jk−Xran3,jk,(57)Zi,jk+1=Yi,jk+1ifrandi≤CRori=rnbri,Xi,jkifrandi>CRandi≠rnbri,(58)Xi,jk+1=Zi,jk+1ifAZi,jk+1≤AXi,jk,Xi,jk,otherwise,where ran1, ran2, and ran3 are three different numbers ∈1,2,…,N, N is the population size, and the randomly chosen numbers ran1, ran2, and ran3 are different from integer i; F is the scaling factor; CR is the crossover constant; and randi and rnbri are random numbers in the interval of [0, 1] and in 1,N, respectively.


## 4.1. Concept of PSO

The PSO is a swarm intelligence-based algorithm proposed by Kennedy and Eberhart [10]. It was derived from observing the bird predation behaviour. PSO uses particles to represent some candidate solutions to the optimization problem in the search space. Position, velocity, and the memory mechanism are three characteristics of each particle [55]. Initially, the particles are randomly generated within the ranges of design variables. The PSO evaluates each particle’s fitness value by calculating the objective function. Then, the position of each design variable is updated by the following equations:(42)Xi,jk+1=Xi,jk+Vi,jk+1.In addition,(43)Vi,jk+1=wkVi,jk+C1⋅r1Pbi,jk−Xi,jk+C2⋅r2gbi,jk−Xi,jk,where Xi,jk is the current jth design variable of particle i at the kth iteration, Vi,jk is the corresponding velocity of the jth design variable, Pbi,jk is the best position at the kth iteration, gbi,jk is the global best position, r1 and r2 are random numbers between [0, 1], C1 and C2 are user-supplied coefficients, and wk is the inertia weight. In this study, the inertia weight is defined as follows:(44)wk=wmax−wmax−wminItmax⋅k,where wmax is the maximum value of the inertia weight and wmin is the minimum value of the inertia weight. Itmax is the maximum generation, which is the termination criterion of this study.


## 4.2. Concept of WOA

The WOA mimics the bubble-net hunting behaviour of humpback whales. The unique hunting model is shown in Figure6. The WOA consists of two phases: the exploitation phase and the exploration phase [21]. The mathematical models of the two stages are discussed in the following sections.Figure 6
Unique bubble-net spiral hunting behaviour.


### 4.2.1. Exploitation Phase

As shown in Figure7, the shrinking encircling mechanism and bubble-net spiral attacking mechanism are included in the exploitation phase. In WOA, the global best position is assumed as the target prey of the current generation, and other candidates will update their positions in a shrinking encircling manner. The mathematical model of the shrinking encircling mechanism is shown as follows:(45)Xi,jk+1=gbi,jk−ADi,jk+1,(46)Di,jk+1=C⋅gbi,jk−Xi,jk,where A and C are coefficients of WOA and can be calculated as follows:(47)A=2ar−a,C=2r,a=2−k2Itmax,where r is a random number in the interval of [0, 1]. To simulate the bubble-net attacking behaviour of humpback whales, a spiral equation is created as follows:(48)Xi,jk+1=Di,j′k+1⋅ebl⋅cos2πl+gbi,jk.Figure 7
Exploitation phase in WOA. (a) Shrinking encircling mechanism. (b) Bubble-net spiral attacking mechanism.
(a)(b)In this equation,Di,j′k+1=gbi,jk−Xi,jk, l is a random number in the interval of [−1, 1], and constant b defines the spiral shape. The humpback whales can swim around their prey in the above two paths at the same time. The mathematical model of this process is shown as follows:(49)Xi,jk+1=gbi,jk−ADi,jk+1if prob<0.5,Di,j′k+1⋅ebl⋅cos2πl+gbi,jkif prob≥0.5.


### 4.2.2. Exploration Phase

To expand the searching area, the WOA updates each design variable’s position using a randomly selected candidate whenA≥1. The mathematical model for the second phase is represented by equations (50) and (51):(50)Xi,jk+1=Xri,jk−ADi,j″k+1,(51)Di,j″k+1=C⋅Xri,jk−Xi,jk,where Xri,jk represents a randomly selected search agent in the current population. The flowchart of WOA is presented in Figure 8(a).Figure 8
Flowcharts of the (a) WOA algorithm and (b) EWOA algorithm.
(a)(b)


## 4.2.1. Exploitation Phase

As shown in Figure7, the shrinking encircling mechanism and bubble-net spiral attacking mechanism are included in the exploitation phase. In WOA, the global best position is assumed as the target prey of the current generation, and other candidates will update their positions in a shrinking encircling manner. The mathematical model of the shrinking encircling mechanism is shown as follows:(45)Xi,jk+1=gbi,jk−ADi,jk+1,(46)Di,jk+1=C⋅gbi,jk−Xi,jk,where A and C are coefficients of WOA and can be calculated as follows:(47)A=2ar−a,C=2r,a=2−k2Itmax,where r is a random number in the interval of [0, 1]. To simulate the bubble-net attacking behaviour of humpback whales, a spiral equation is created as follows:(48)Xi,jk+1=Di,j′k+1⋅ebl⋅cos2πl+gbi,jk.Figure 7
Exploitation phase in WOA. (a) Shrinking encircling mechanism. (b) Bubble-net spiral attacking mechanism.
(a)(b)In this equation,Di,j′k+1=gbi,jk−Xi,jk, l is a random number in the interval of [−1, 1], and constant b defines the spiral shape. The humpback whales can swim around their prey in the above two paths at the same time. The mathematical model of this process is shown as follows:(49)Xi,jk+1=gbi,jk−ADi,jk+1if prob<0.5,Di,j′k+1⋅ebl⋅cos2πl+gbi,jkif prob≥0.5.


## 4.2.2. Exploration Phase

To expand the searching area, the WOA updates each design variable’s position using a randomly selected candidate whenA≥1. The mathematical model for the second phase is represented by equations (50) and (51):(50)Xi,jk+1=Xri,jk−ADi,j″k+1,(51)Di,j″k+1=C⋅Xri,jk−Xi,jk,where Xri,jk represents a randomly selected search agent in the current population. The flowchart of WOA is presented in Figure 8(a).Figure 8
Flowcharts of the (a) WOA algorithm and (b) EWOA algorithm.
(a)(b)


## 4.3. Concept of EWOA

The WOA has an effective global search capability. In order to enhance the WOA in terms of reliability, EWOA was proposed. The EWOA maintains the simplicity of the WOA. In the exploitation phase, to maintain a balance between the diversity and intensification of search results, equation (45) is replaced as follows [23]:(52)Xi,jk+1=gbi,jk−ADi,j″′k+1,(53)Di,j″′k+1=rXi,jk.In the exploration phase, each candidate randomly changes the value of a design variable with a probabilityp″:(54)p″=0.31−kItmax.For candidatei, the design variable is changed when a random number q is less than probability p″. The mathematical model of this process is shown as follows [23]:(55)Xi,jk+1=Xjmin+random⋅Xjmax−Xjmin,q<p″,Xi,jk,q≥p″,where random and q represent random numbers between [0, 1], respectively. Xjmin and Xjmax represent the minimum value and the maximum value of the jth design variable, respectively. The flowchart of EWOA is presented in Figure 8(b).


## 4.4. Concept of DE

The DE is a simple population-based evolutionary algorithm proposed by Storn and Price [25]. Similar to PSO, it starts from randomly generating an initial population within the ranges of design variables [56]. The mutation process is to create a mutated vector Yi,jk+1 according to equation (56). The crossover process is to create a trial vector Zi,jk+1, as described in equation (57) [25]. In the selection process, whether the trial vector should be selected for the next generation can be determined according to equation (58) [8].(56)Yi,jk+1=Xran1,jk+FXran2,jk−Xran3,jk,(57)Zi,jk+1=Yi,jk+1ifrandi≤CRori=rnbri,Xi,jkifrandi>CRandi≠rnbri,(58)Xi,jk+1=Zi,jk+1ifAZi,jk+1≤AXi,jk,Xi,jk,otherwise,where ran1, ran2, and ran3 are three different numbers ∈1,2,…,N, N is the population size, and the randomly chosen numbers ran1, ran2, and ran3 are different from integer i; F is the scaling factor; CR is the crossover constant; and randi and rnbri are random numbers in the interval of [0, 1] and in 1,N, respectively.


## 5. The Proposed Algorithm

### 5.1. PSO-EWOA-DE Algorithm Description

The PSO-EWOA-DE algorithm is a population-based evolutionary algorithm. Initially, the PSO-EWOA-DE algorithm explores the search space globally by using PSO and EWOA and then utilizes DE to fine-tune the selected solution locally. This process is repeated until the stopping criterion is met. This algorithm can ensure diversity and maintain convergence accuracy. The PSO-EWOA-DE algorithm steps are summarized as follows:(1)
Set the population sizeN and the value range of the design variables. Then, generate the initialized population for design variables.(2)
Calculate the fitness value for each candidate.(3)
Update the parametersa, r, k, l, p″, and q in the EWOA process and the parameters r1, r2, and wk in the PSO process.(4)
Generate the trial populationPOP_EWOA according to the EWOA process shown in Section 4.3, and calculate the fitness values APOP_EWOA of population POP_EWOA.(5)
Generate the trial populationPOP_PSO according to the PSO process shown in Section 4.1, and calculate the fitness values APOP_PSO of population POP_PSO.(6)
Select the trial populationPOP_trial according to equation (59) and the global best individual Gbest according to equation (60).(59)POP_triali=POP_EWOAiifAPOP_EWOAi≤APOP_PSOi,POP_PSOi,otherwise,(60)Gbest=Gbest_EWOAifAGbest_EWOA≤AGbest_PSO,Gbest_PSO,otherwise,wherePOP_triali is the ith candidate of the trial population, POP_EWOAi is the ith candidate of POP_EWOA, POP_PSOi is the ith candidate of POP_PSO, Gbest_EWOA is the global best individual of POP_EWOA, and Gbest_PSO is the global best individual of POP_PSO.(7)
Generate a new population for the next generation according to the DE process shown in Section4.4.(8)
Terminate the PSO-EWOA-DE algorithm by satisfying the termination condition. If so, go to Step 9; otherwise, go to Step 2.(9)
Output the optimal solution obtained in Step 6.The flowchart of the PSO-EWOA-DE algorithm is presented in Figure9(a).Figure 9
Flowcharts of the (a) PSO-EWOA-DE algorithm and (b) structural optimization procedure.
(a)(b)


### 5.2. Optimization Procedure

This section presents a structural optimization procedure of the ship lock head structure based on the PSO-EWOA-DE algorithm. The optimization procedure was implemented by using a self-developed Python script, which was worked in Abaqus 6.14. The heat of hydration was calculated by using a user subroutine (HETVAL). Furthermore, a user subroutine (UMAT) was developed to implement the MCC model and calculate the thermal creep stress.The structural optimization process of a ship lock head on soft soil, considering the time-varying effect of the structure and foundation, can be summarized as follows:(1)
The population sizeN and the value range of the design variables were determined. Then, the initialized population for 34 design variables was generated.(2)
The FEM model for a ship lock head was generated. Furthermore, the initial temperature, convection conditions, construction steps, and the heat of hydration were applied to calculate the temperature field.(3)
The element type was automatically modified to calculate the initial geostress field.(4)
The pore water pressure boundaries, displacement boundaries, and design loads were applied to calculate the stress field by calling Abaqus.(5)
Three types of constraints shown inSection 3.3 were evaluated, and the fitness value was calculated.(6)
A new population was generated using the PSO-EWOA-DE algorithm shown inSection 5.1.(7)
The stopping criterion was checked, and Steps 2 to 6 were repeated whenk<Itmax.(8)
The optimum solution with the highest fitness obtained in Step 6 was output.The flowchart of the structural optimization procedure is shown in Figure9(b).


## 5.1. PSO-EWOA-DE Algorithm Description

The PSO-EWOA-DE algorithm is a population-based evolutionary algorithm. Initially, the PSO-EWOA-DE algorithm explores the search space globally by using PSO and EWOA and then utilizes DE to fine-tune the selected solution locally. This process is repeated until the stopping criterion is met. This algorithm can ensure diversity and maintain convergence accuracy. The PSO-EWOA-DE algorithm steps are summarized as follows:(1)
Set the population sizeN and the value range of the design variables. Then, generate the initialized population for design variables.(2)
Calculate the fitness value for each candidate.(3)
Update the parametersa, r, k, l, p″, and q in the EWOA process and the parameters r1, r2, and wk in the PSO process.(4)
Generate the trial populationPOP_EWOA according to the EWOA process shown in Section 4.3, and calculate the fitness values APOP_EWOA of population POP_EWOA.(5)
Generate the trial populationPOP_PSO according to the PSO process shown in Section 4.1, and calculate the fitness values APOP_PSO of population POP_PSO.(6)
Select the trial populationPOP_trial according to equation (59) and the global best individual Gbest according to equation (60).(59)POP_triali=POP_EWOAiifAPOP_EWOAi≤APOP_PSOi,POP_PSOi,otherwise,(60)Gbest=Gbest_EWOAifAGbest_EWOA≤AGbest_PSO,Gbest_PSO,otherwise,wherePOP_triali is the ith candidate of the trial population, POP_EWOAi is the ith candidate of POP_EWOA, POP_PSOi is the ith candidate of POP_PSO, Gbest_EWOA is the global best individual of POP_EWOA, and Gbest_PSO is the global best individual of POP_PSO.(7)
Generate a new population for the next generation according to the DE process shown in Section4.4.(8)
Terminate the PSO-EWOA-DE algorithm by satisfying the termination condition. If so, go to Step 9; otherwise, go to Step 2.(9)
Output the optimal solution obtained in Step 6.The flowchart of the PSO-EWOA-DE algorithm is presented in Figure9(a).Figure 9
Flowcharts of the (a) PSO-EWOA-DE algorithm and (b) structural optimization procedure.
(a)(b)


## 5.2. Optimization Procedure

This section presents a structural optimization procedure of the ship lock head structure based on the PSO-EWOA-DE algorithm. The optimization procedure was implemented by using a self-developed Python script, which was worked in Abaqus 6.14. The heat of hydration was calculated by using a user subroutine (HETVAL). Furthermore, a user subroutine (UMAT) was developed to implement the MCC model and calculate the thermal creep stress.The structural optimization process of a ship lock head on soft soil, considering the time-varying effect of the structure and foundation, can be summarized as follows:(1)
The population sizeN and the value range of the design variables were determined. Then, the initialized population for 34 design variables was generated.(2)
The FEM model for a ship lock head was generated. Furthermore, the initial temperature, convection conditions, construction steps, and the heat of hydration were applied to calculate the temperature field.(3)
The element type was automatically modified to calculate the initial geostress field.(4)
The pore water pressure boundaries, displacement boundaries, and design loads were applied to calculate the stress field by calling Abaqus.(5)
Three types of constraints shown inSection 3.3 were evaluated, and the fitness value was calculated.(6)
A new population was generated using the PSO-EWOA-DE algorithm shown inSection 5.1.(7)
The stopping criterion was checked, and Steps 2 to 6 were repeated whenk<Itmax.(8)
The optimum solution with the highest fitness obtained in Step 6 was output.The flowchart of the structural optimization procedure is shown in Figure9(b).


## 6. Engineering Example of a Ship Lock Head

### 6.1. Finite Element Model Description

The structural optimization of the proposed ship lock head during the construction period considered only the static load. The design loads of the lock head-soft foundation-backfilled soil system were gravity loads, active earth pressure, live loads, and temperature. A lock head with a width of 53.8 m and a height of 12.2 m was presented here. Table1 lists the initial values for design variables of the lock head. The lower and upper limits for design variables are also shown in Table 1. To reduce the calculation time and improve optimization efficiency, we set value ranges for 34 design variables. The value ranges for 34 design variables were set according to the engineering example’s initial sizes introduced in this paper. It ensured that each searched ship lock head’s volume was smaller than the initial design during the optimization process.Table 1
Range of design variables.
Design variableLower limits (m)Upper limits (m)Initial values (m)Optimum values (m)x13.294.594.594.3x20.691.291.290.69x30.060.240.240.102x43.713.913.713.91x50.060.710.710.506x60.060.690.690.6x7−1.31−0.51−0.51−0.78x817.9218.1917.9218.139x90.1650.190.1650.19x103.753.893.753.89x113.793.93.793.88x1222.923.3922.923.39x130.060.690.690.69x1423.2923.4923.2923.30x150.060.390.390.37y10.060.690.690.456y20.060.590.590.201y30.490.890.890.506y40.061.811.810.348y50.111.161.160.152y611.5911.911.911.6y70.061.511.510.605y80.060.690.690.62y90.060.610.610.45y100.390.490.490.42y110.060.710.710.6y120.061.891.891.15z10.061.601.600.3z20.060.290.290.06z30.060.230.230.06z40.431.021.020.49z50.060.520.520.3z60.060.310.310.07z70.060.570.570.06The lock head density is 2400 kg/m3, and Poisson’s ratio is 0.167. According to the previous research [4], the elastic modulus Et′, the creep function Ct,t′, and the heat of hydration Qt′ of concrete are as follows:(61)Et′=34.251−exp−0.28τ0.52×109,Ct,τ=6.39×1+9.2τ−0.451−exp−0.3t−τ+14.5×1+1.7τ−0.451−exp−0.005t−τ×10−12,Qt′=98167.481−exp−0.0033t′1.525.Table2 shows the main material properties of the soft foundation and backfilled soil. Equation (62) gives the air temperature curve during the construction of the lock head:(62)Ta=18.87+10.83cosπτ−132180.Table 2
Main material properties of the soil foundation and backfilled soil.
DefinitionSoft foundationBackfilled soilMass density (kg/m3)19001900Poisson’s ratio0.30.3Initial void ratio0.8370.837Virgin compression index0.08530.082Unloading-reloading index0.007140.00697Permeability coefficient0.7760.003Critical-state parameter0.9331.502Preconsolidation pressure (KPa)450.587191.361The entire construction process of the ship lock head was 386 days. As shown in Table3, the construction schedule was divided into 11 steps [4]. According to Kennedy and Eberhart [10], the user-supplied coefficients C1 and C2 of PSO were set to 2 in this study. Meanwhile, the maximum value wmax and the minimum value wmin of the inertia weight were 0.9 and 0.4, respectively [57]. The scaling factor F and the crossover constant CR were 0.5 and 0.55, respectively [28]. The population size N was 40, and the maximum iteration Itmax was 30.Table 3
Ship lock construction process.
Construction sequenceConstruction processConstruction duration (day)Construction intermissions (day)C1Pouring bottom plate28C2Removing bottom plate formwork135C3Backfilling the first layer of soil2096C4Pouring the corridor layer110C5Removing corridor layer formwork125C6Backfilling the second layer of soil1926C7Pouring the empty-box layer19C8Removing empty-box layer formwork153C9Pouring the second-stage concrete18C10Removing the second-stage concrete formwork120C11Backfilling the third layer of soil1830


### 6.2. Optimization Results

In this section, the previously described procedure was used to optimize the ship lock head. The optimization work met the requirements of China. For 60 runs of PSO [10], WOA [15], DE [25], EWOA [23], PSO-DE [28], and PSO-EWOA-DE algorithms, each algorithm’s convergence history is depicted in Figure 10. The PSO-EWOA-DE converged after 13 iterations, which was faster than others. Furthermore, PSO-EWOA-DE produced the best solution when other algorithms were trapped into a local optimal solution, meaning PSO-EWOA-DE was more suitable for the structural optimization of the lock head.Figure 10
Evolutionary processes of PSO, WOA, and EWOA algorithms.Table1 lists the optimum design variables of the lock head. Compared with the initial design, it can be found that the optimal design reduced the height of the corridor layer’s coverslip and correspondingly increased the height of the empty box. Therefore, the volume of the corridor layer was significantly reduced, and the volume of the empty-box layer was slightly increased. The results of the optimization are listed in Table 4. As seen, the PSO-EWOA-DE leads to a 10.45% reduction in the structure volume when compared with the initial design.Table 4
Results of optimization.
Volume (m3)PSOWOADEEWOAPSO-DEPSO-EWOA-DEInitial design3338.043338.043338.043338.043338.043338.04Optimal design3031.813019.493011.392994.072992.552989.30Decreasing ratio9.17%9.54%9.79%10.30%10.35%10.45%The antisliding safety factors and the anti-overturning safety factors for the optimal results of the PSO, WOA, DE, EWOA, PSO-DE, and PSO-EWOA-DE are summarized in Figure11. The optimal result of the PSO-EWOA-DE had the smallest safety factors under safe conditions. It verified that the optimal volume of the ship lock head searched by the PSO-EWOA-DE was the smallest.Figure 11
Comparisons of antisliding safety factors and anti-overturning safety factors.Figure12 demonstrates the convergence histories of the PSO-EWOA-DE algorithm and the EWOA algorithm when the population sizes N were 10, 20, 30, 40, and 50, respectively. As shown in Figure 12(a), it could be seen that, with the increase in population sizes, the results searched by the PSO-EWOA-DE algorithm gradually approached the optimal solution. The optimum volume of the lock head was 2989.30 m3. When the population sizes were 10, 20, and 30, the population sizes were too small to find the optimum solution. When the population sizes were 40 and 50, the PSO-EWOA-DE converged at the 13th and 11th iteration, respectively. The number of numerical simulations of the latter was 30 runs more than that of the former. As shown in Figure 12(b), when the population sizes were 40 and 50, the EWOA converged at the 19th and 18th iteration, respectively. The number of numerical simulations of the latter was 140 runs more than that of the former. It could be concluded that population size 40 was more suitable for optimizing the lock head.Figure 12
Effects of population sizeN on the (a) PSO-EWOA-DE’s convergence history and (b) EWOA’s convergence history.
(a)(b)


### 6.3. Analysis of Stresses

Figure13 identifies four feature points A, B, C, and D, located at the side bottom plate, corridor layer, empty-box layer, and middle bottom plate. These four points were used to analyse the lock head’s stress fields and displacement fields at different construction steps. Meanwhile, stress and vertical displacement changes of the lock head before and after optimization were also analysed.Figure 13
Locations of four feature points.It can be found from Figure14 that the maximum principle stresses of the ship lock head were significantly affected by the temperature load at early ages. When the concrete was poured, due to the heat released by cement hydration, the tensile stresses increased rapidly with the increase of elastic modulus. The rise in the corridor layer’s tensile stress reached 0.33 MPa at early age, which was the largest. It can also be found from Figure 14 that when backfilling the soil, under the influence of the backfilling soil’s gravity and the soil foundation’s pore water pressure, the tensile stresses of the lock head suddenly increased. The increase of the side bottom plate’s tensile stress reached 2.27 MPa, which was the largest. It was concluded that the ship lock head should be protected against cracking when pouring the concrete and backfilling the soil.Figure 14
Comparisons of maximum (max.) principle stress before and after optimization at (a) pointA, (b) point B, (c) point C, and (d) point D.
(a)(b)(c)(d)Figure14(a) shows the maximum principle stress change of feature point A on the side bottom plate before and after optimization. When backfilling the soil, the optimal design’s side bottom plate’s tensile stress became more extensive due to reducing the corridor layer’s volume. The maximum increase value was 0.42 MPa when backfilling the third layer of soil. Figure 14(b) shows that the optimal design corridor layer’s tensile stress was reduced by 0.21 MPa at early age. The thinning of the corridor layer’s coverslip increased the heat dissipation rate and reduced the thermal creep stress. As shown in Figure 14(c), due to the increase in the empty-box layer’s volume, the tensile stress of the optimal design’s empty-box layer increased by 0.14 MPa at early age. However, the tensile stress was reduced by 0.17 MPa when backfilling the soil. As shown in Figure 14(d), structural optimization had little effect on the middle bottom plate’s stress field. It was concluded that the structural optimization method that comprehensively considered the concrete creep and soft soil consolidation could reflect the time-varying characteristics of the lock head’s stress field.


### 6.4. Displacement Fields’ Analysis

The selected feature points’ time-history curves of the vertical displacements before and after optimization are shown in Figure15. It can be found from Figure 15 that the lock head settled instantaneously when pouring the concrete. Subsequently, as the soft foundation’s pore water pressure dissipated, the lock head’s vertical displacements changed slowly. Since the side bottom plate has the largest volume among the five components that consisted of the lock head, the maximum settlement occurred at the side bottom plate, reaching 320 mm. After pouring the bottom plate, the side bottom plate and the middle bottom plate’s vertical displacements rebounded due to the exothermic expansion caused by a large amount of concrete hydration. The vertical rebound displacement at the side bottom plate reached 21 mm. It was concluded that the analysis method that comprehensively considered the concrete creep and soft soil consolidation could reflect the time-varying characteristics of the lock head’s displacement field.Figure 15
Comparisons of vertical displacements before and after optimization at (a) pointA, (b) point B, (c) point C, and (d) point D.
(a)(b)(c)(d)Figure15(a) shows the vertical displacement change of feature point A on the side bottom plate before and after optimization. It could be found that the settlement of the optimal design’s side bottom plate became smaller due to reducing the volume of the corridor layer. The settlement of the side bottom plate was reduced by 11 mm. Figure 15(b) shows that the optimal design corridor layer’s settlement was reduced by 11 mm due to the lock head’s size changes. As shown in Figure 15(c), due to the increase in the empty-box layer’s volume, the optimal design’s empty-box layer’s settlement was increased by 4 mm. As shown in Figure 15(d), the optimal design’s middle bottom plate’s settlement was decreased by 9 mm. In general, structural optimization could noticeably reduce the lock head’s settlement, which was beneficial to engineering safety.


## 6.1. Finite Element Model Description

The structural optimization of the proposed ship lock head during the construction period considered only the static load. The design loads of the lock head-soft foundation-backfilled soil system were gravity loads, active earth pressure, live loads, and temperature. A lock head with a width of 53.8 m and a height of 12.2 m was presented here. Table1 lists the initial values for design variables of the lock head. The lower and upper limits for design variables are also shown in Table 1. To reduce the calculation time and improve optimization efficiency, we set value ranges for 34 design variables. The value ranges for 34 design variables were set according to the engineering example’s initial sizes introduced in this paper. It ensured that each searched ship lock head’s volume was smaller than the initial design during the optimization process.Table 1
Range of design variables.
Design variableLower limits (m)Upper limits (m)Initial values (m)Optimum values (m)x13.294.594.594.3x20.691.291.290.69x30.060.240.240.102x43.713.913.713.91x50.060.710.710.506x60.060.690.690.6x7−1.31−0.51−0.51−0.78x817.9218.1917.9218.139x90.1650.190.1650.19x103.753.893.753.89x113.793.93.793.88x1222.923.3922.923.39x130.060.690.690.69x1423.2923.4923.2923.30x150.060.390.390.37y10.060.690.690.456y20.060.590.590.201y30.490.890.890.506y40.061.811.810.348y50.111.161.160.152y611.5911.911.911.6y70.061.511.510.605y80.060.690.690.62y90.060.610.610.45y100.390.490.490.42y110.060.710.710.6y120.061.891.891.15z10.061.601.600.3z20.060.290.290.06z30.060.230.230.06z40.431.021.020.49z50.060.520.520.3z60.060.310.310.07z70.060.570.570.06The lock head density is 2400 kg/m3, and Poisson’s ratio is 0.167. According to the previous research [4], the elastic modulus Et′, the creep function Ct,t′, and the heat of hydration Qt′ of concrete are as follows:(61)Et′=34.251−exp−0.28τ0.52×109,Ct,τ=6.39×1+9.2τ−0.451−exp−0.3t−τ+14.5×1+1.7τ−0.451−exp−0.005t−τ×10−12,Qt′=98167.481−exp−0.0033t′1.525.Table2 shows the main material properties of the soft foundation and backfilled soil. Equation (62) gives the air temperature curve during the construction of the lock head:(62)Ta=18.87+10.83cosπτ−132180.Table 2
Main material properties of the soil foundation and backfilled soil.
DefinitionSoft foundationBackfilled soilMass density (kg/m3)19001900Poisson’s ratio0.30.3Initial void ratio0.8370.837Virgin compression index0.08530.082Unloading-reloading index0.007140.00697Permeability coefficient0.7760.003Critical-state parameter0.9331.502Preconsolidation pressure (KPa)450.587191.361The entire construction process of the ship lock head was 386 days. As shown in Table3, the construction schedule was divided into 11 steps [4]. According to Kennedy and Eberhart [10], the user-supplied coefficients C1 and C2 of PSO were set to 2 in this study. Meanwhile, the maximum value wmax and the minimum value wmin of the inertia weight were 0.9 and 0.4, respectively [57]. The scaling factor F and the crossover constant CR were 0.5 and 0.55, respectively [28]. The population size N was 40, and the maximum iteration Itmax was 30.Table 3
Ship lock construction process.
Construction sequenceConstruction processConstruction duration (day)Construction intermissions (day)C1Pouring bottom plate28C2Removing bottom plate formwork135C3Backfilling the first layer of soil2096C4Pouring the corridor layer110C5Removing corridor layer formwork125C6Backfilling the second layer of soil1926C7Pouring the empty-box layer19C8Removing empty-box layer formwork153C9Pouring the second-stage concrete18C10Removing the second-stage concrete formwork120C11Backfilling the third layer of soil1830


## 6.2. Optimization Results

In this section, the previously described procedure was used to optimize the ship lock head. The optimization work met the requirements of China. For 60 runs of PSO [10], WOA [15], DE [25], EWOA [23], PSO-DE [28], and PSO-EWOA-DE algorithms, each algorithm’s convergence history is depicted in Figure 10. The PSO-EWOA-DE converged after 13 iterations, which was faster than others. Furthermore, PSO-EWOA-DE produced the best solution when other algorithms were trapped into a local optimal solution, meaning PSO-EWOA-DE was more suitable for the structural optimization of the lock head.Figure 10
Evolutionary processes of PSO, WOA, and EWOA algorithms.Table1 lists the optimum design variables of the lock head. Compared with the initial design, it can be found that the optimal design reduced the height of the corridor layer’s coverslip and correspondingly increased the height of the empty box. Therefore, the volume of the corridor layer was significantly reduced, and the volume of the empty-box layer was slightly increased. The results of the optimization are listed in Table 4. As seen, the PSO-EWOA-DE leads to a 10.45% reduction in the structure volume when compared with the initial design.Table 4
Results of optimization.
Volume (m3)PSOWOADEEWOAPSO-DEPSO-EWOA-DEInitial design3338.043338.043338.043338.043338.043338.04Optimal design3031.813019.493011.392994.072992.552989.30Decreasing ratio9.17%9.54%9.79%10.30%10.35%10.45%The antisliding safety factors and the anti-overturning safety factors for the optimal results of the PSO, WOA, DE, EWOA, PSO-DE, and PSO-EWOA-DE are summarized in Figure11. The optimal result of the PSO-EWOA-DE had the smallest safety factors under safe conditions. It verified that the optimal volume of the ship lock head searched by the PSO-EWOA-DE was the smallest.Figure 11
Comparisons of antisliding safety factors and anti-overturning safety factors.Figure12 demonstrates the convergence histories of the PSO-EWOA-DE algorithm and the EWOA algorithm when the population sizes N were 10, 20, 30, 40, and 50, respectively. As shown in Figure 12(a), it could be seen that, with the increase in population sizes, the results searched by the PSO-EWOA-DE algorithm gradually approached the optimal solution. The optimum volume of the lock head was 2989.30 m3. When the population sizes were 10, 20, and 30, the population sizes were too small to find the optimum solution. When the population sizes were 40 and 50, the PSO-EWOA-DE converged at the 13th and 11th iteration, respectively. The number of numerical simulations of the latter was 30 runs more than that of the former. As shown in Figure 12(b), when the population sizes were 40 and 50, the EWOA converged at the 19th and 18th iteration, respectively. The number of numerical simulations of the latter was 140 runs more than that of the former. It could be concluded that population size 40 was more suitable for optimizing the lock head.Figure 12
Effects of population sizeN on the (a) PSO-EWOA-DE’s convergence history and (b) EWOA’s convergence history.
(a)(b)


## 6.3. Analysis of Stresses

Figure13 identifies four feature points A, B, C, and D, located at the side bottom plate, corridor layer, empty-box layer, and middle bottom plate. These four points were used to analyse the lock head’s stress fields and displacement fields at different construction steps. Meanwhile, stress and vertical displacement changes of the lock head before and after optimization were also analysed.Figure 13
Locations of four feature points.It can be found from Figure14 that the maximum principle stresses of the ship lock head were significantly affected by the temperature load at early ages. When the concrete was poured, due to the heat released by cement hydration, the tensile stresses increased rapidly with the increase of elastic modulus. The rise in the corridor layer’s tensile stress reached 0.33 MPa at early age, which was the largest. It can also be found from Figure 14 that when backfilling the soil, under the influence of the backfilling soil’s gravity and the soil foundation’s pore water pressure, the tensile stresses of the lock head suddenly increased. The increase of the side bottom plate’s tensile stress reached 2.27 MPa, which was the largest. It was concluded that the ship lock head should be protected against cracking when pouring the concrete and backfilling the soil.Figure 14
Comparisons of maximum (max.) principle stress before and after optimization at (a) pointA, (b) point B, (c) point C, and (d) point D.
(a)(b)(c)(d)Figure14(a) shows the maximum principle stress change of feature point A on the side bottom plate before and after optimization. When backfilling the soil, the optimal design’s side bottom plate’s tensile stress became more extensive due to reducing the corridor layer’s volume. The maximum increase value was 0.42 MPa when backfilling the third layer of soil. Figure 14(b) shows that the optimal design corridor layer’s tensile stress was reduced by 0.21 MPa at early age. The thinning of the corridor layer’s coverslip increased the heat dissipation rate and reduced the thermal creep stress. As shown in Figure 14(c), due to the increase in the empty-box layer’s volume, the tensile stress of the optimal design’s empty-box layer increased by 0.14 MPa at early age. However, the tensile stress was reduced by 0.17 MPa when backfilling the soil. As shown in Figure 14(d), structural optimization had little effect on the middle bottom plate’s stress field. It was concluded that the structural optimization method that comprehensively considered the concrete creep and soft soil consolidation could reflect the time-varying characteristics of the lock head’s stress field.


## 6.4. Displacement Fields’ Analysis

The selected feature points’ time-history curves of the vertical displacements before and after optimization are shown in Figure15. It can be found from Figure 15 that the lock head settled instantaneously when pouring the concrete. Subsequently, as the soft foundation’s pore water pressure dissipated, the lock head’s vertical displacements changed slowly. Since the side bottom plate has the largest volume among the five components that consisted of the lock head, the maximum settlement occurred at the side bottom plate, reaching 320 mm. After pouring the bottom plate, the side bottom plate and the middle bottom plate’s vertical displacements rebounded due to the exothermic expansion caused by a large amount of concrete hydration. The vertical rebound displacement at the side bottom plate reached 21 mm. It was concluded that the analysis method that comprehensively considered the concrete creep and soft soil consolidation could reflect the time-varying characteristics of the lock head’s displacement field.Figure 15
Comparisons of vertical displacements before and after optimization at (a) pointA, (b) point B, (c) point C, and (d) point D.
(a)(b)(c)(d)Figure15(a) shows the vertical displacement change of feature point A on the side bottom plate before and after optimization. It could be found that the settlement of the optimal design’s side bottom plate became smaller due to reducing the volume of the corridor layer. The settlement of the side bottom plate was reduced by 11 mm. Figure 15(b) shows that the optimal design corridor layer’s settlement was reduced by 11 mm due to the lock head’s size changes. As shown in Figure 15(c), due to the increase in the empty-box layer’s volume, the optimal design’s empty-box layer’s settlement was increased by 4 mm. As shown in Figure 15(d), the optimal design’s middle bottom plate’s settlement was decreased by 9 mm. In general, structural optimization could noticeably reduce the lock head’s settlement, which was beneficial to engineering safety.


## 7. Conclusions

The structural optimization procedure for ship lock heads on soft soil, considering the time-varying effects of the structure and foundation, was established in this paper. By comprehensively considering the linear viscoelastic creep of concrete and the elastoplastic consolidation characteristic of soft soil, a perfect time-dependent analysis method for the lock head on soft soil was proposed. The MCC model was applied in the consolidation calculation of the soft foundation. To obtain the minimum volume, the PSO-EWOA-DE algorithm was proposed to optimize thirty-four design variables of a lock head. Among all the compared algorithms, the PSO-EWOA-DE algorithm showed promising results and had a faster convergence rate. The following conclusions are drawn by analysing the maximum principle stresses and vertical displacements of the lock head:(1)
The temperature load and the backfilling soil’s gravity are the main reasons for increasing the ship lock head’s tensile stress on the soft foundation. The lock head should be protected against cracking when pouring the concrete and backfilling the soil.(2)
The lock head settled instantaneously when pouring the concrete. Subsequently, as the soft foundation’s pore water pressure dissipated, the lock head’s vertical displacements changed slowly. The maximum settlement occurred at the side bottom plate.(3)
In general, structural optimization could noticeably reduce the lock head’s settlement during construction, which was beneficial to engineering safety. The analysis method that comprehensively considered the concrete creep and soft soil consolidation could reflect the lock head’s stress field and displacement field’s time-varying characteristics.This structural optimization method of the ship lock head considered the time-varying effects of the structure and foundation. Future work should consider the optimization of the ship lock head’s construction sequence.



---

*Source: 5517060-2021-04-13.xml*