# Parameter Estimation for Coupled Hydromechanical Simulation of Dynamic Compaction Based on Pareto Multiobjective Optimization

Authors: Wei Wang; Lu-Lu Zhang; Jin-Jian Chen; Jian-Hua Wang
Journal: Shock and Vibration
 (2015)

---

## Abstract

This paper presented a parameter estimation method based on a coupled hydromechanical model of dynamic compaction and the Pareto multiobjective optimization technique. The hydromechanical model of dynamic compaction is established in the FEM program LS-DYNA. The multiobjective optimization algorithm, Nondominated Sorted Genetic Algorithm (NSGA-IIa), is integrated with the numerical model to identify soil parameters using multiple sources of field data. A field case study is used to demonstrate the capability of the proposed method. The observed pore water pressure and crater depth at early blow of dynamic compaction are simultaneously used to estimate the soil parameters. Robustness of the back estimated parameters is further illustrated by a forward prediction. Results show that the back-analyzed soil parameters can reasonably predict lateral displacements and give generally acceptable predictions of dynamic compaction for an adjacent location. In addition, for prediction of ground response of the dynamic compaction at continuous blows, the prediction based on the second blow is more accurate than the first blow due to the occurrence of the hardening and strengthening of soil during continuous compaction.


---

## Body

## 1. Introduction

Dynamic compaction (DC) is a ground improvement technique in which a heavy hammer with a weightWt is dropped on the ground from a specific height Ht repeatedly to densify soft ground, as illustrated in Figure 1. This method has been widely used because it is (1) easy to implement, (2) economically competitive (no material added or replaced), and (3) environmentally safe (no chemicals injected or added). Originally, the predominant soil types considered for dynamic compaction treatment are well-drained sandy and granular materials [1]. Recently, the dynamic compaction method has been applied in various civil engineering projects at coastal areas, such as highway construction (e.g., [2, 3]), land reclamation (e.g., [4, 5]), port construction (e.g., [6, 7]), and airport construction (e.g., [8, 9]). As the groundwater table is near the ground surface in the coastal areas, a significant increase of pore water pressure is noticed after each impact, which results in a reduction of bearing capacity and liquefaction may occur during dynamic compaction in saturated sandy deposits [1]. Consequently, the densification effect of heavy weight decreases dramatically with further drops.Figure 1
Schematic of dynamic compaction.The generation of excess pore water pressure and soil deformation due to dynamic compaction in saturated soils are interacted [10, 11]. Numerical modeling provides a rational tool to describe the complex coupled behavior during dynamic compaction [10, 11]. A coupled hydromechanical simulation can examine the ground response to dynamic compaction for both solid and fluid phases. However, accurate prediction using numerical simulation depends greatly on the selection of constitutive models and the determination of soil parameters [12–14]. Firstly, the soil parameters in coupled modeling of dynamic compaction are usually estimated from the standard laboratory or field tests for static soil behavior due to the high cost of laboratory dynamic tests [11, 13–15]. Secondly, the laboratory or field measured soil parameters are subject to uncertainties of sample disturbance and measurement error. Even with accurately measured soil parameters, the predicted performance using numerical modeling may still deviate from the field observation due to inherent spatial variability of in situ soils [16–18] and model error [19].A rational approach to calibrate a numerical model is back estimation of soil parameters based on the field observations. For dynamic compaction problem in saturated soil, both soil displacement and pore water pressure are usually measured in situ. To consider different types of measurements simultaneously, the back analysis may be implemented within a multiobjective optimization framework. In recent years, the multiobjective optimization with the Pareto solutions [20–22], which represent trade-offs among the incommensurable and conflicting objectives, is gradually attracting research attention in geotechnical problems. Papon et al. [23] presented a multiobjective parameter identification method based on MOGA II and back-analyzed soil parameters using the data from two pressuremeter tests. Huang et al. [24] developed a back analysis method based on AMALGAM algorithm and the Pareto multiobjective optimization for deep excavation. To the best knowledge of the authors, parameter estimation of dynamic compaction based on multiple types of field data, especially using the Pareto multiobjective optimization method, is very limited in the literature.In this study, a parameter estimation method using multiple sources of field data based on a coupled hydromechanical model of dynamic compaction and the Pareto multiobjective optimization technique is presented. The hydromechanical model of dynamic compaction is established in the FEM program LS-DYNA. The multiobjective optimization algorithm, Nondominated Sorted Genetic Algorithm (NSGA-IIa), is integrated with the numerical model to identify soil parameters. The application of the method is illustrated by an example of highway construction project. The trade-offs of the two objectives are analyzed based on the Pareto optimal solutions. The robustness of the back-estimated parameters is illustrated by a forward prediction of the performance for an adjacent test point. The contribution of this study is the integration of a coupled hydromechanical model with the Pareto multiobjective optimization technique in back analysis of soil parameters for dynamic compaction problems. In addition, the proposed parameter estimation method for dynamic compaction problems may provide a rational tool to determine in situ dynamic properties of soils.


## 2. Coupled Hydromechanical Model for Dynamic Compaction

### 2.1. Governing Equations of Coupled Hydromechanical Analysis

The mechanical behavior of saturated soils subjected to dynamic loads is governed by the interaction of the solid skeleton with the pore fluid [25]. In this paper, the extended Biot dynamic equations [26] are adopted for modeling of dynamic compaction. The governing equations, which represent equilibrium of a porous medium, equilibrium of fluid phrase, and the mass conservation of fluid phase, are as follows:(1)σij,j-ρu¨i-ρfw¨i+ρgi=0,(2)p,i-ρfgikw˙i-ρfu¨i-ρfnw¨i+ρfgi=0,(3)u˙i,i+w˙i,i+nKfp˙=0,where σij is total stress which is further separated into the effective stress σij′ acting on the solid skeleton and pore pressure p for the fluid in the pores. gi is body force acceleration, ui is displacement of the solid skeleton, wi is displacement of pore fluid, k is permeability, and n is porosity of the mixture. The density ρf is that of the pore fluid and ρ is the density of the mixture written as ρ=1-nρs+nρf, where ρs is the density of the solid skeleton. Kf is bulk modulus of the fluid.These governing equations should be solved simultaneously to consider the interdependency of skeleton displacementui, fluid displacement wi, and pore pressure p. Assuming that the fluid acceleration relative to the solid skeleton w¨i is negligible and substituting w˙i in (2) into (3), a simplified U-P dynamic form [26] is taken in this study. Therefore, only skeleton displacement ui and pore pressure p are to be solved through the governing equations. The spatial discretized form of the governing equations in the U-P form can be expressed as follows:(4)MU¨+CU˙+KU-QP=FU,(5)GU¨+QTU˙+SP˙+HP=FP,where M is global mass matrix, C is viscous damping matrix, K is stiffness matrix, Q is coupling matrix, G is dynamic seepage force matrix, S is compressibility matrix, and H is permeability matrix. U and P are the solid displacement and pore pressure vectors, respectively. The vectors FU and FP include the effect of external forces and fluid flux, respectively.In (4), the term KU in spatial domain Ω also can be expressed as ∫ΩBTσ′dΩ, where B is strain-displacement matrix and σ′ is the effective stress tensor (determined by soil constitutive model which will be discussed later). In this study, the viscous damping of the soil has been modelled via the Rayleigh approach, in which the damping matrix C is computed as a linear combination of the mass (M) and stiffness (K) matrices:(6)C=αmM+βkK,where αm and βk are Rayleigh damping parameters.


### 2.2. Soil Constitutive Law

Most numerical models for dynamic compaction assumed elastic soil behavior [10, 27–29]. Considering permanent plastic deformations observed in the compaction process and the possibility of yield due to compression, the cap model [30] was adopted to model soil behavior under impact loads, for example, for dynamic compaction in granular soils by Gu and Lee [13], Lee and Gu [14], and Ghassemi et al. [11, 15]. The role of the cap, which is the yield surface due to compression, is to enhance a classic elastoplastic model by modeling the phenomenon of soil compaction where the soil mass yields under high compressive stress. In this paper, an extended two-invariant cap model based on the formulations of Sandler and Rubin [30] and Simo et al. [31] was used.The incremental effective stress changedσij′ of constitutive relations of the skeleton can be defined in terms of(7)dσij′=Dijkldεije=Dijkldεij-dεijp,where dεije is the incremental elastic strain, dεij is the incremental total strain, dεijp is the incremental plastic strain, and Dijkl is the elastic stress matrix coefficients with Dijkl=K-(2/3)Gδijδkl+Gδikδjl+δilδjk, where K is the bulk modulus, G is the shear modulus, and δij is the Kronecker delta.The cap model is an elastoplastic model. An associated law is assumed for the plastic deformation potential. The plastic yield surfaces of the cap model are composed of a shear failure surfacefE, a movable cap surface fC, and a tension cutoff surface fT, as shown in Figure 2.Figure 2
The cap model.The failure surface, denoted byfEI1,J2 in Figure 2, is a nonhardening modified Drucker-Prager form with a yield function defined as(8)fE=J2-α+γexp⁡-βI1-θI1=0,where I1 is the first stress invariant and J2 is the second invariant of deviatoric stress tensor. α, γ, β, and θ are the material parameters. If the material parameter γ=β=0, the failure surface, fE, becomes a simple Drucker-Prager surface.The movable cap, denoted byfCI1,J2,κ in Figure 2, is a hardening elliptical surface defined as(9)fC=J2-1RXκ-Lκ2-I1-Lκ2=0,where R is the curvature of the hardening cap and κ is hardening parameter related to the actual plastic volumetric change εvp(10)εvpXκ=W1-exp⁡-DXκ-X0,where W and D are the material parameters. X0 is the initial stress state of underground. X(κ) defines the intersection of the cap surface with the I1 axis(11)Xκ=κ+Rα-γexp⁡-βκ+θκand L(κ) is the value of I1 at the location of the start of cap(12)Lκ=κifκ>00ifκ≤0.The tension cutoff surface, denoted byfTI1, is represented by(13)fT=T-I1=0,where T is maximum hydrostatic tension. Since most soils cannot support significant tensile stress, a value of tension cutoff is assumed to be zero in this study.When the cap model is adopted for a dynamic compaction problem where the strain hardening behavior due to compaction load is dominant, model parameters related to expansion of the cap surface, such as the size of cap in elastoplastic model (W, D, and R) and the elastic moduli (K and G), are most important. The parameters, W, D, and R, depend on the plastic deformation and can be determined from laboratory triaxial tests [32]. However, the bulk modulus,K, and the shear modulus, G, may vary during dynamic compaction as the soil density increases significantly in the vicinity of the contact surface of impact load [28]. Previous research studies [11–15, 28, 33] had made efforts to estimate the varied elastic properties at the end of each impact. They found that the bulk modulus K and the shear modulus G during dynamic compaction are stress or strain dependent and deviate much from static moduli [11, 13, 14]. Therefore, K and G are selected as the target soil parameters to be identified in the back analysis of dynamic compaction in this study.


### 2.3. Modeling of Impact Load

There are two methods commonly adopted to simulate dynamic load in numerical modeling [12]. The first one is to impose the hammer with a vertical velocity to strike the ground surface, which can be determined from the free fall equation (as shown in Figure 1). Another method is to simplify the dynamic load as a normal distribution curve or a damped half-sine wave [12, 15]. In this study, the first method is used and a symmetric penalty function algorithm [34] is adopted to obtain the contact force in the numerical model.


## 2.1. Governing Equations of Coupled Hydromechanical Analysis

The mechanical behavior of saturated soils subjected to dynamic loads is governed by the interaction of the solid skeleton with the pore fluid [25]. In this paper, the extended Biot dynamic equations [26] are adopted for modeling of dynamic compaction. The governing equations, which represent equilibrium of a porous medium, equilibrium of fluid phrase, and the mass conservation of fluid phase, are as follows:(1)σij,j-ρu¨i-ρfw¨i+ρgi=0,(2)p,i-ρfgikw˙i-ρfu¨i-ρfnw¨i+ρfgi=0,(3)u˙i,i+w˙i,i+nKfp˙=0,where σij is total stress which is further separated into the effective stress σij′ acting on the solid skeleton and pore pressure p for the fluid in the pores. gi is body force acceleration, ui is displacement of the solid skeleton, wi is displacement of pore fluid, k is permeability, and n is porosity of the mixture. The density ρf is that of the pore fluid and ρ is the density of the mixture written as ρ=1-nρs+nρf, where ρs is the density of the solid skeleton. Kf is bulk modulus of the fluid.These governing equations should be solved simultaneously to consider the interdependency of skeleton displacementui, fluid displacement wi, and pore pressure p. Assuming that the fluid acceleration relative to the solid skeleton w¨i is negligible and substituting w˙i in (2) into (3), a simplified U-P dynamic form [26] is taken in this study. Therefore, only skeleton displacement ui and pore pressure p are to be solved through the governing equations. The spatial discretized form of the governing equations in the U-P form can be expressed as follows:(4)MU¨+CU˙+KU-QP=FU,(5)GU¨+QTU˙+SP˙+HP=FP,where M is global mass matrix, C is viscous damping matrix, K is stiffness matrix, Q is coupling matrix, G is dynamic seepage force matrix, S is compressibility matrix, and H is permeability matrix. U and P are the solid displacement and pore pressure vectors, respectively. The vectors FU and FP include the effect of external forces and fluid flux, respectively.In (4), the term KU in spatial domain Ω also can be expressed as ∫ΩBTσ′dΩ, where B is strain-displacement matrix and σ′ is the effective stress tensor (determined by soil constitutive model which will be discussed later). In this study, the viscous damping of the soil has been modelled via the Rayleigh approach, in which the damping matrix C is computed as a linear combination of the mass (M) and stiffness (K) matrices:(6)C=αmM+βkK,where αm and βk are Rayleigh damping parameters.


## 2.2. Soil Constitutive Law

Most numerical models for dynamic compaction assumed elastic soil behavior [10, 27–29]. Considering permanent plastic deformations observed in the compaction process and the possibility of yield due to compression, the cap model [30] was adopted to model soil behavior under impact loads, for example, for dynamic compaction in granular soils by Gu and Lee [13], Lee and Gu [14], and Ghassemi et al. [11, 15]. The role of the cap, which is the yield surface due to compression, is to enhance a classic elastoplastic model by modeling the phenomenon of soil compaction where the soil mass yields under high compressive stress. In this paper, an extended two-invariant cap model based on the formulations of Sandler and Rubin [30] and Simo et al. [31] was used.The incremental effective stress changedσij′ of constitutive relations of the skeleton can be defined in terms of(7)dσij′=Dijkldεije=Dijkldεij-dεijp,where dεije is the incremental elastic strain, dεij is the incremental total strain, dεijp is the incremental plastic strain, and Dijkl is the elastic stress matrix coefficients with Dijkl=K-(2/3)Gδijδkl+Gδikδjl+δilδjk, where K is the bulk modulus, G is the shear modulus, and δij is the Kronecker delta.The cap model is an elastoplastic model. An associated law is assumed for the plastic deformation potential. The plastic yield surfaces of the cap model are composed of a shear failure surfacefE, a movable cap surface fC, and a tension cutoff surface fT, as shown in Figure 2.Figure 2
The cap model.The failure surface, denoted byfEI1,J2 in Figure 2, is a nonhardening modified Drucker-Prager form with a yield function defined as(8)fE=J2-α+γexp⁡-βI1-θI1=0,where I1 is the first stress invariant and J2 is the second invariant of deviatoric stress tensor. α, γ, β, and θ are the material parameters. If the material parameter γ=β=0, the failure surface, fE, becomes a simple Drucker-Prager surface.The movable cap, denoted byfCI1,J2,κ in Figure 2, is a hardening elliptical surface defined as(9)fC=J2-1RXκ-Lκ2-I1-Lκ2=0,where R is the curvature of the hardening cap and κ is hardening parameter related to the actual plastic volumetric change εvp(10)εvpXκ=W1-exp⁡-DXκ-X0,where W and D are the material parameters. X0 is the initial stress state of underground. X(κ) defines the intersection of the cap surface with the I1 axis(11)Xκ=κ+Rα-γexp⁡-βκ+θκand L(κ) is the value of I1 at the location of the start of cap(12)Lκ=κifκ>00ifκ≤0.The tension cutoff surface, denoted byfTI1, is represented by(13)fT=T-I1=0,where T is maximum hydrostatic tension. Since most soils cannot support significant tensile stress, a value of tension cutoff is assumed to be zero in this study.When the cap model is adopted for a dynamic compaction problem where the strain hardening behavior due to compaction load is dominant, model parameters related to expansion of the cap surface, such as the size of cap in elastoplastic model (W, D, and R) and the elastic moduli (K and G), are most important. The parameters, W, D, and R, depend on the plastic deformation and can be determined from laboratory triaxial tests [32]. However, the bulk modulus,K, and the shear modulus, G, may vary during dynamic compaction as the soil density increases significantly in the vicinity of the contact surface of impact load [28]. Previous research studies [11–15, 28, 33] had made efforts to estimate the varied elastic properties at the end of each impact. They found that the bulk modulus K and the shear modulus G during dynamic compaction are stress or strain dependent and deviate much from static moduli [11, 13, 14]. Therefore, K and G are selected as the target soil parameters to be identified in the back analysis of dynamic compaction in this study.


## 2.3. Modeling of Impact Load

There are two methods commonly adopted to simulate dynamic load in numerical modeling [12]. The first one is to impose the hammer with a vertical velocity to strike the ground surface, which can be determined from the free fall equation (as shown in Figure 1). Another method is to simplify the dynamic load as a normal distribution curve or a damped half-sine wave [12, 15]. In this study, the first method is used and a symmetric penalty function algorithm [34] is adopted to obtain the contact force in the numerical model.


## 3. Back Analysis of Dynamic Compaction Using Pareto Multiobjective Optimization

### 3.1. Multiobjective Optimization of Back Analysis

The match of simulated results to the observations is one of the most important indicators of how well a model represents an actual system. The difference between the simulated result and observed result is measured by an objective function. To determine the optimized soil parameters by back analysis, the soil parameters are updated repeatedly through an iteration process until the simulated result is close enough to the observed result, while the uncertainties of the initial and boundary conditions, the structural inadequacies of the prediction model, uncertainties of input parameters and measurement errors, and the residual errors are not expected to be equal to zero. The objective function of a back analysis is usually defined as a function of the residual errors. In this paper, the objective function is formulated as the root mean square error (RMSE):(14)fX=∑i=1dyi′-yi2,where X is n sized vector of unknown model input parameters, X={x1,x2,…,xn} and is constrained to variable bounds (xi(L)≤xi≤xi(U)), the parameter d is the number of measurement points that are used to back figure the target parameters in the optimization procedure, and yi and yi′ are the simulated and measured value at measurement point i, respectively. When multiple types of measurements are to be considered in the back analysis for dynamic compaction simultaneously, the objective function can be written as a multiobjective optimization form [35, 36]:(15)MinimizeFX=f1X,f2X,…,fmX,…,fMX,where FX is the vector of objective functions. fmX, m=1,2,…,M, is the mth objective function. For example, f1(X), f2(X), and f3(X) can be the values of the objective function in terms of the relative error of pore water pressure, crater depth, and lateral displacement, respectively.The objective functions of multiobjective optimization constitute a multidimensional objective space. Different objective functions are often conflicting against each other and it is not possible to satisfy all the goals at a time. Hence, the multiobjective optimization yields a set of Pareto optimal solutions [20–22]. The Pareto optimal solutions are defined as follows. Let X(1) and X(2) be two parameter vectors. X(1) is said to dominate X(2), if both the following conditions are true: (1)fm(X(1)) ≤ fm(X(2)) for all m=1,2,…,M and (2)fm(X(1)) < fm(X(2)) for at least one m. In other words, the solution X(1) is no worse than X(2) in all objectives and X(1) is strictly better than X(2) in at least one objective. A solution is said to be Pareto optimal if it is not dominated by any other solution in the solution space.The set of Pareto optimal solutions forms a front in the objective space and the front is called the Pareto front. In the simplest case of two objective functions, the Pareto front is a curve in a two-dimensional space. Figure3 illustrates the Pareto solution set for a simple problem where the aim is to simultaneously minimize two objectives (f1(X),f2(X)). The black square dots indicate an initial set of parameter estimates. Moving along the Pareto front results in smaller values of one objective but at the expense of increasing of the second objective [22, 37]. Therefore, all the solutions in the Pareto front are considered as optimal solutions. The shape of this Pareto front between two objectives is very informative. A sharp Pareto front (Figure 4(a)) indicates that the model can simultaneously meet both objectives adequately [37]. On the other hand, a rather linear Pareto front (Figure 4(b)) indicates that the model is not capable of simulating both objectives simultaneously [37].Figure 3
Illustration of the Pareto optimal points and the Pareto front for a problem with two objective functions.Figure 4
Illustration of (a) a sharp Pareto front and (b) a rather linear Pareto front.
(a)
(b)


### 3.2. Multiobjective Optimization Algorithm and Interaction with LSDYNA

The multiobjective optimization is challenging because a reliable and exhaustive set of Pareto solutions should be determined with limited computation costs [37]. In recent years, a number of multiobjective optimization algorithms based on the artificial intelligence and the evolutionary computational techniques emerge and are applied in various engineering problems [22, 38, 39]. Some widely used algorithms are the Multiobjective Genetic Algorithm (MOGA) [40], the Niched Pareto Genetic Algorithm (NPGA) [41], the Multiobjective Complex Evolution (MOCOM-UA) method [42], the Nondominated Sorted Genetic Algorithm (NSGA-II) [43], Strength Pareto Evolutionary Algorithm (SPEA2) [44], and the multiobjective Evolutionary Algorithm (FFEA) by Fonseca and Fleming [45].Among all multiobjective optimization algorithms, Nondominated Sorted Genetic Algorithm (NSGA-II) has been reported to perform well for a large variety of problems of varying complexities because NSGA-II result in a diverse set of trade-off solutions in a single simulation run. However, this algorithm uses a finite population size; there may be the Pareto drift problem [46]. To eliminate this problem, in this study, the Nondominated Sorted Genetic Algorithm (NSGA-II) [43] augmented with an archiving strategy [46] is adopted. The archived NSGA-II method (NSGA-IIa) combines two concepts, that is, a population-based elitism search and maintaining an external archive to save potential Pareto optimal set. A random initial population of size s and an empty archive are first generated. Then, each parent from the initial population is assigned a rank using the Fast Nondominated Sorting algorithm [43]. The population of offspring (with size s) is subsequently created using crossover and mutation operators. The offsprings and parents are then ranked together and members for the next population are chosen from the Pareto front based on their rank and crowding distance [43]. Simultaneously, all rank = 1 solutions are added to an external archive, and the archive is updated by removing all dominated and duplicate solutions. The new population is then used to generate offspring and the aforementioned algorithmic steps are repeated until convergence is achieved [43]. The performance of NSGA-IIa is evaluated by two criteria [47]. The first is convergence to a global Paretooptimal front. The second is diversity (or spread) on the Paretooptimal front. More details about the descriptions of NSGA-IIa can be found in Goel et al. [46] and hence are not repeated here.In this study, the multiobjective back analysis for DC is proposed by integrating the NSGA-IIa algorithm with the coupled hydromechanical model of dynamic compaction. The coupled hydromechanical model is established in the commercial FEM package LSDYNA. The interaction between multiobjective optimization and the finite element modeling of LSDYNA in the back analysis for DC is shown in the flowchart in Figure5. In the NSGA-IIa code, the parent population is selected and the offspring population is developed through mutation and crossover. After that, the two populations are imported into LSDYNA as input parameters for the FEM analysis. The results of numerical analysis are used to calculate the objective function. The population of the next generation is selected and is then used to generate offsprings. When the change of objective function is less than the fixed tolerance (10−4), the inverse analysis procedure is terminated.Figure 5
Flow chart of multiobjective back analysis.


## 3.1. Multiobjective Optimization of Back Analysis

The match of simulated results to the observations is one of the most important indicators of how well a model represents an actual system. The difference between the simulated result and observed result is measured by an objective function. To determine the optimized soil parameters by back analysis, the soil parameters are updated repeatedly through an iteration process until the simulated result is close enough to the observed result, while the uncertainties of the initial and boundary conditions, the structural inadequacies of the prediction model, uncertainties of input parameters and measurement errors, and the residual errors are not expected to be equal to zero. The objective function of a back analysis is usually defined as a function of the residual errors. In this paper, the objective function is formulated as the root mean square error (RMSE):(14)fX=∑i=1dyi′-yi2,where X is n sized vector of unknown model input parameters, X={x1,x2,…,xn} and is constrained to variable bounds (xi(L)≤xi≤xi(U)), the parameter d is the number of measurement points that are used to back figure the target parameters in the optimization procedure, and yi and yi′ are the simulated and measured value at measurement point i, respectively. When multiple types of measurements are to be considered in the back analysis for dynamic compaction simultaneously, the objective function can be written as a multiobjective optimization form [35, 36]:(15)MinimizeFX=f1X,f2X,…,fmX,…,fMX,where FX is the vector of objective functions. fmX, m=1,2,…,M, is the mth objective function. For example, f1(X), f2(X), and f3(X) can be the values of the objective function in terms of the relative error of pore water pressure, crater depth, and lateral displacement, respectively.The objective functions of multiobjective optimization constitute a multidimensional objective space. Different objective functions are often conflicting against each other and it is not possible to satisfy all the goals at a time. Hence, the multiobjective optimization yields a set of Pareto optimal solutions [20–22]. The Pareto optimal solutions are defined as follows. Let X(1) and X(2) be two parameter vectors. X(1) is said to dominate X(2), if both the following conditions are true: (1)fm(X(1)) ≤ fm(X(2)) for all m=1,2,…,M and (2)fm(X(1)) < fm(X(2)) for at least one m. In other words, the solution X(1) is no worse than X(2) in all objectives and X(1) is strictly better than X(2) in at least one objective. A solution is said to be Pareto optimal if it is not dominated by any other solution in the solution space.The set of Pareto optimal solutions forms a front in the objective space and the front is called the Pareto front. In the simplest case of two objective functions, the Pareto front is a curve in a two-dimensional space. Figure3 illustrates the Pareto solution set for a simple problem where the aim is to simultaneously minimize two objectives (f1(X),f2(X)). The black square dots indicate an initial set of parameter estimates. Moving along the Pareto front results in smaller values of one objective but at the expense of increasing of the second objective [22, 37]. Therefore, all the solutions in the Pareto front are considered as optimal solutions. The shape of this Pareto front between two objectives is very informative. A sharp Pareto front (Figure 4(a)) indicates that the model can simultaneously meet both objectives adequately [37]. On the other hand, a rather linear Pareto front (Figure 4(b)) indicates that the model is not capable of simulating both objectives simultaneously [37].Figure 3
Illustration of the Pareto optimal points and the Pareto front for a problem with two objective functions.Figure 4
Illustration of (a) a sharp Pareto front and (b) a rather linear Pareto front.
(a)
(b)


## 3.2. Multiobjective Optimization Algorithm and Interaction with LSDYNA

The multiobjective optimization is challenging because a reliable and exhaustive set of Pareto solutions should be determined with limited computation costs [37]. In recent years, a number of multiobjective optimization algorithms based on the artificial intelligence and the evolutionary computational techniques emerge and are applied in various engineering problems [22, 38, 39]. Some widely used algorithms are the Multiobjective Genetic Algorithm (MOGA) [40], the Niched Pareto Genetic Algorithm (NPGA) [41], the Multiobjective Complex Evolution (MOCOM-UA) method [42], the Nondominated Sorted Genetic Algorithm (NSGA-II) [43], Strength Pareto Evolutionary Algorithm (SPEA2) [44], and the multiobjective Evolutionary Algorithm (FFEA) by Fonseca and Fleming [45].Among all multiobjective optimization algorithms, Nondominated Sorted Genetic Algorithm (NSGA-II) has been reported to perform well for a large variety of problems of varying complexities because NSGA-II result in a diverse set of trade-off solutions in a single simulation run. However, this algorithm uses a finite population size; there may be the Pareto drift problem [46]. To eliminate this problem, in this study, the Nondominated Sorted Genetic Algorithm (NSGA-II) [43] augmented with an archiving strategy [46] is adopted. The archived NSGA-II method (NSGA-IIa) combines two concepts, that is, a population-based elitism search and maintaining an external archive to save potential Pareto optimal set. A random initial population of size s and an empty archive are first generated. Then, each parent from the initial population is assigned a rank using the Fast Nondominated Sorting algorithm [43]. The population of offspring (with size s) is subsequently created using crossover and mutation operators. The offsprings and parents are then ranked together and members for the next population are chosen from the Pareto front based on their rank and crowding distance [43]. Simultaneously, all rank = 1 solutions are added to an external archive, and the archive is updated by removing all dominated and duplicate solutions. The new population is then used to generate offspring and the aforementioned algorithmic steps are repeated until convergence is achieved [43]. The performance of NSGA-IIa is evaluated by two criteria [47]. The first is convergence to a global Paretooptimal front. The second is diversity (or spread) on the Paretooptimal front. More details about the descriptions of NSGA-IIa can be found in Goel et al. [46] and hence are not repeated here.In this study, the multiobjective back analysis for DC is proposed by integrating the NSGA-IIa algorithm with the coupled hydromechanical model of dynamic compaction. The coupled hydromechanical model is established in the commercial FEM package LSDYNA. The interaction between multiobjective optimization and the finite element modeling of LSDYNA in the back analysis for DC is shown in the flowchart in Figure5. In the NSGA-IIa code, the parent population is selected and the offspring population is developed through mutation and crossover. After that, the two populations are imported into LSDYNA as input parameters for the FEM analysis. The results of numerical analysis are used to calculate the objective function. The population of the next generation is selected and is then used to generate offsprings. When the change of objective function is less than the fixed tolerance (10−4), the inverse analysis procedure is terminated.Figure 5
Flow chart of multiobjective back analysis.


## 4. Example Application: A Field Test of Dynamic Compaction Treatment in a Highway

In this paper, a field case of dynamic compaction treatment in a highway construction reported by Miao et al. [2, 3] is used as an illustrative example of the multiobjective back analysis method. The dynamic compaction case is selected because a comprehensive field monitoring program was conducted in this site and various types of field observations were obtained.


### 4.1. Project Description and Monitoring Program

The ground treatment project is one section of the Lianyungang-Xuzhou Highway, which is located in Jiangsu Province, China. The site ground is composed of the liquefiable silty soil and soft clay. The soil layers and the corresponding soil properties are presented in Table1. For compaction test, a 19-ton hammer with a base radius of 1.3 m was used. Two applied tamping energy (Wt×Ht) levels for single blow, 1500 kN m (drop height Ht is 8.05 m) and 2500 kN m (Ht is 13.42 m), were used at two adjacent testing locations (point A and point B) as shown in Figure 6. The number of drops was set at five for each impact location.Table 1
Soil profile at the site (sourced from Miao et al. 2006 [3]).
Soil layer
Depth range (m)
Colour
Nature water content (%)
Plasticity index
Note
Clayed silt, silt
0.0~4.0
Yellow
25.2
8.9
Liquefaction
Silt
4.0~9.0
Grey
24.7
7.2
Liquefaction
Silty clay
9.0~10.0
Grey
25.6
15.4
Clay
10.0~11.0
Dark grey
28.6
17.1
Clay
>11.0
Red, drab
33.6
18.6Figure 6
The geology profile and plan layout of monitoring system.To investigate the effectiveness of dynamic compaction for improving liquefiable soil, a monitoring system was established. The crater depth (surface settlement), lateral displacements, and pore water pressure are measured (Figure6). The measurement points for pore pressure were at depths of 3, 5, 7, and 9 m and at a distance of 4 m from the compaction point.The oscillating pore water pressure right after high energy impact is difficult to measure. In this field test, the recorded data are available only at seconds after impacts, as shown in Figure7. Herein, crater depth is the best control parameter of ground improvement, and the generated pore water pressure can represent the coupled behavior of solid and fluid phase. Thus, in this example, two objective functions, f1(X) and f2(X) which are RMSEs of pore water pressure and crater depth, respectively, are used to back-analyze soil parameters.Figure 7
Illustration of the measured data of (a) pore water pressure with each blow of dynamic compaction and (b) lateral displacement at different depth with each blow of dynamic compaction.
(a)
(b)


### 4.2. Numerical Model and Material Parameters

Figure8 shows the 3D finite element mesh used in this analysis. To reduce the computational load, the soil body was modeled as a quarter of cylinder and consisted of three soil layers. The first layer is corresponding to the liquefiable clayey silt and silt soil layer. The second layer consisted of the silt layer. The third layer is composed of the silty clay and clay, since the silty clay is less thick and is close to the clay in soil characteristics. The overall model size is 15 meters in radius and 19 meters in height. The boundary on the cylindrical side is fixed along a horizontal direction, and the vertical direction is free. The horizontal and rotational movements are fixed on the two sides of symmetric boundary. The bottom boundary is fixed along both horizontal and vertical directions. Moreover, nonreflecting boundaries are used on the exterior boundaries around the soil model to prevent artificial stress wave reflections generated at the model boundaries from reentering the model and contaminating the results. The size of the model is selected to avoid interference from the constrained boundary conditions. The total number of brick elements to model soil layers is 180,768. The hammer is modeled as a circular disc using rigid elements and is constrained in all degrees of freedom (translational and rotational) except the vertical one. The cap model is adopted for the first two soil layers. The third soil layer is modeled as an elastic material as limited plastic deformation may occur in this soil layer. The Young’s modulus E of the third layer is estimated from the SPT value and assumed to be 7 × 103 kPa. The Poisson ratio ν is assumed to be 0.3, which is a typical value for clayey soils.Figure 8
Schematic diagrams of (a) FEM model and boundary conditions and (b) profile of the A-A section.
(a)
(b)The elastic bulk modulusK and shear modulus G are identified as the most important parameters of dynamic compaction modeling and have significant effects on the hardening behavior of soils during compaction [11, 13]. Thus, adequate information on dynamic soil properties, especially dynamic elastic modulus, is essential for accurate computations of ground response [48]. Since the third layer affects performance of dynamic compaction less significantly compared with the first two layers in this project, the soil parameters of the first two soil layers are estimated in the back analysis. In this case study, a total of four soil parameters are chosen to be estimated in the back analysis. The two parameters, that is, K and G for the ith soil layer, are denoted as Ki and Gi, respectively. Table 2 lists the upper and lower bounds of these two parameters to cover most possible values [11]. The coefficient of permeability and the bulk modulus of pore water are assumed to be 1 × 10−7m/s and 1 × 104 kPa based on Ghassemi et al. [11]. In addition, Rayleigh damping is set at 8% of critical damping, with lower and upper reference frequencies of 2 Hz and 30 Hz, this being consistent with the findings about vibrational frequency in the dynamic compaction test by Lukas [49, 50]. Therefore, input parameters of the numerical model are listed in Table 2.Table 2
Input data for dynamic compaction analysis.
(a)
Coupled hydromechanical model
k (m/s)
n
ρ
f (kg/m3)
ρ
s (kg/m3)
K
f (kPa)
1 × 10−7
0.5
1000
2650
1 × 104
(b)
Cap model
Soil layer
K (kPa)
G (kPa)
α (kPa)
γ
β
θ
W
D (kPa−1)
R
T
1
Blow 1
[8000, 30000]
[2000, 8000]
17.8
0
0
0.25
0.5
0.0003
4.33
0
Blow 2
[8000, 500000]
[1000, 10000]
2
Blow 1
[8000, 30000]
[2000, 8000]
21.2
0
0
0.21
0.5
0.0003
4.33
0
Blow 2
[8000, 500000]
[1000, 10000]
(c)
Elastic model
Soil layer
E (kPa)
ν
3
7 × 103
0.3
(d)
Rayleigh damping parameters
α
k
β
m
1.884
0.0008


### 4.3. Procedures of Parameter Estimation and Prediction

The main algorithmic parameter of NSGA-IIa is the population size,s. Previous multiobjective benchmark examples [46] have demonstrated that NSGA-IIa works well with a population size of about 50. In this example, we use a population size of s=80 for each generation. The maximum number of generation is set as 100. For each blow of dynamic compaction, a maximum of 8000 simulations (100 generations) are evaluated in the FE analysis program. Due to the high computation load, only the measured pore water pressures and the crater depths of the first and second blows (denoted as blow 1 and blow 2) at impact location A are used for back-estimate soil parameters. The back-analyzed parameters are then used to predict lateral displacement and the change of soil properties during subsequent blows at impact location A. Prediction for the performance at impact location B is also conducted.


## 4.1. Project Description and Monitoring Program

The ground treatment project is one section of the Lianyungang-Xuzhou Highway, which is located in Jiangsu Province, China. The site ground is composed of the liquefiable silty soil and soft clay. The soil layers and the corresponding soil properties are presented in Table1. For compaction test, a 19-ton hammer with a base radius of 1.3 m was used. Two applied tamping energy (Wt×Ht) levels for single blow, 1500 kN m (drop height Ht is 8.05 m) and 2500 kN m (Ht is 13.42 m), were used at two adjacent testing locations (point A and point B) as shown in Figure 6. The number of drops was set at five for each impact location.Table 1
Soil profile at the site (sourced from Miao et al. 2006 [3]).
Soil layer
Depth range (m)
Colour
Nature water content (%)
Plasticity index
Note
Clayed silt, silt
0.0~4.0
Yellow
25.2
8.9
Liquefaction
Silt
4.0~9.0
Grey
24.7
7.2
Liquefaction
Silty clay
9.0~10.0
Grey
25.6
15.4
Clay
10.0~11.0
Dark grey
28.6
17.1
Clay
>11.0
Red, drab
33.6
18.6Figure 6
The geology profile and plan layout of monitoring system.To investigate the effectiveness of dynamic compaction for improving liquefiable soil, a monitoring system was established. The crater depth (surface settlement), lateral displacements, and pore water pressure are measured (Figure6). The measurement points for pore pressure were at depths of 3, 5, 7, and 9 m and at a distance of 4 m from the compaction point.The oscillating pore water pressure right after high energy impact is difficult to measure. In this field test, the recorded data are available only at seconds after impacts, as shown in Figure7. Herein, crater depth is the best control parameter of ground improvement, and the generated pore water pressure can represent the coupled behavior of solid and fluid phase. Thus, in this example, two objective functions, f1(X) and f2(X) which are RMSEs of pore water pressure and crater depth, respectively, are used to back-analyze soil parameters.Figure 7
Illustration of the measured data of (a) pore water pressure with each blow of dynamic compaction and (b) lateral displacement at different depth with each blow of dynamic compaction.
(a)
(b)


## 4.2. Numerical Model and Material Parameters

Figure8 shows the 3D finite element mesh used in this analysis. To reduce the computational load, the soil body was modeled as a quarter of cylinder and consisted of three soil layers. The first layer is corresponding to the liquefiable clayey silt and silt soil layer. The second layer consisted of the silt layer. The third layer is composed of the silty clay and clay, since the silty clay is less thick and is close to the clay in soil characteristics. The overall model size is 15 meters in radius and 19 meters in height. The boundary on the cylindrical side is fixed along a horizontal direction, and the vertical direction is free. The horizontal and rotational movements are fixed on the two sides of symmetric boundary. The bottom boundary is fixed along both horizontal and vertical directions. Moreover, nonreflecting boundaries are used on the exterior boundaries around the soil model to prevent artificial stress wave reflections generated at the model boundaries from reentering the model and contaminating the results. The size of the model is selected to avoid interference from the constrained boundary conditions. The total number of brick elements to model soil layers is 180,768. The hammer is modeled as a circular disc using rigid elements and is constrained in all degrees of freedom (translational and rotational) except the vertical one. The cap model is adopted for the first two soil layers. The third soil layer is modeled as an elastic material as limited plastic deformation may occur in this soil layer. The Young’s modulus E of the third layer is estimated from the SPT value and assumed to be 7 × 103 kPa. The Poisson ratio ν is assumed to be 0.3, which is a typical value for clayey soils.Figure 8
Schematic diagrams of (a) FEM model and boundary conditions and (b) profile of the A-A section.
(a)
(b)The elastic bulk modulusK and shear modulus G are identified as the most important parameters of dynamic compaction modeling and have significant effects on the hardening behavior of soils during compaction [11, 13]. Thus, adequate information on dynamic soil properties, especially dynamic elastic modulus, is essential for accurate computations of ground response [48]. Since the third layer affects performance of dynamic compaction less significantly compared with the first two layers in this project, the soil parameters of the first two soil layers are estimated in the back analysis. In this case study, a total of four soil parameters are chosen to be estimated in the back analysis. The two parameters, that is, K and G for the ith soil layer, are denoted as Ki and Gi, respectively. Table 2 lists the upper and lower bounds of these two parameters to cover most possible values [11]. The coefficient of permeability and the bulk modulus of pore water are assumed to be 1 × 10−7m/s and 1 × 104 kPa based on Ghassemi et al. [11]. In addition, Rayleigh damping is set at 8% of critical damping, with lower and upper reference frequencies of 2 Hz and 30 Hz, this being consistent with the findings about vibrational frequency in the dynamic compaction test by Lukas [49, 50]. Therefore, input parameters of the numerical model are listed in Table 2.Table 2
Input data for dynamic compaction analysis.
(a)
Coupled hydromechanical model
k (m/s)
n
ρ
f (kg/m3)
ρ
s (kg/m3)
K
f (kPa)
1 × 10−7
0.5
1000
2650
1 × 104
(b)
Cap model
Soil layer
K (kPa)
G (kPa)
α (kPa)
γ
β
θ
W
D (kPa−1)
R
T
1
Blow 1
[8000, 30000]
[2000, 8000]
17.8
0
0
0.25
0.5
0.0003
4.33
0
Blow 2
[8000, 500000]
[1000, 10000]
2
Blow 1
[8000, 30000]
[2000, 8000]
21.2
0
0
0.21
0.5
0.0003
4.33
0
Blow 2
[8000, 500000]
[1000, 10000]
(c)
Elastic model
Soil layer
E (kPa)
ν
3
7 × 103
0.3
(d)
Rayleigh damping parameters
α
k
β
m
1.884
0.0008


## 4.3. Procedures of Parameter Estimation and Prediction

The main algorithmic parameter of NSGA-IIa is the population size,s. Previous multiobjective benchmark examples [46] have demonstrated that NSGA-IIa works well with a population size of about 50. In this example, we use a population size of s=80 for each generation. The maximum number of generation is set as 100. For each blow of dynamic compaction, a maximum of 8000 simulations (100 generations) are evaluated in the FE analysis program. Due to the high computation load, only the measured pore water pressures and the crater depths of the first and second blows (denoted as blow 1 and blow 2) at impact location A are used for back-estimate soil parameters. The back-analyzed parameters are then used to predict lateral displacement and the change of soil properties during subsequent blows at impact location A. Prediction for the performance at impact location B is also conducted.


## 5. Results and Discussion

### 5.1. Coupled Hydromechanical Response in Numerical Modeling

In this section, we present the FEM results including soil displacements and generated pore water pressure to illustrate the coupled hydromechanical responses due to dynamic compaction. The variation of lateral soil displacement with time is shown in Figure9. In this figure, the soil deformation near the edge of the hammer increases progressively, as the wave front propagates away from the hammer. In Figure 10, the generation of pore water pressure and its dissipation are illustrated. At the first stage of loading, during the beginning of the impact, the pore water pressure is increased instantaneously below the loading zone. As the wave progresses, the excess pore pressure moves downwards and decreases gradually below the loading zone.Figure 9
Contour lines of lateral displacement for the first blow: (a) 0.05 s, (b) 0.1 s, (c) 0.2 s, and (d) 0.3 s.
(a)
(b)
(c)
(d)Figure 10
Contour lines of pore water pressure for the first blow: (a) 0.05 s, (b) 0.1 s, (c) 0.2 s, and (d) 0.3 s.
(a)
(b)
(c)
(d)


### 5.2. Characteristics of the Pareto Front

Due to the page limit, we present only the result based on the measured data of blow 1 at impact location A. Figure11 shows the convergence of optimization with an increasing number of generations. The objective function f2 converges to less than 2 × 10−3mm at the 33th generation. The objective function f1 converges at 3.29 kPa at the 37th generation. The back analysis is terminated at the 68th generation.Figure 11
Evolution of the minimum objective function values as a function of the number of generations.Figure12 presents the trade-off of the f1-f2 biobjective Pareto front for the model. All the numerical simulated solutions are marked as gray square points. The Pareto solutions of the last generation are marked as blue circles. The Pareto front exhibits a rectangular shape, indicating that it is possible to minimize both objectives simultaneously. Here, the most informative Pareto point is selected from the Pareto solutions. This point is the knee point in the Pareto front, where a small improvement in one objective would lead to a large deterioration in at least one other objective [51]. Such a point can be regarded as the optimal solution of the multiobjective optimization problem and is more relevant to the decision maker than the entire Pareto front. In the graph, we denote the point as Pc and the solution is subsequently referred to as the compromise solution.Figure 12
Pareto front and optimal solutions.


### 5.3. Identification of Soil Parameters

Figure13 shows normalized parameter ranges corresponding to the Pareto solutions that appeared in previous section. The back-estimated parameters, K and G, of the first two soil layers are listed along the abscissa. The subscript i of each parameter corresponds to the ith soil layer. The ordinate presents the parameter values normalized to the prior bounds (as given in Table 2). Each blue line in the graph represents one Pareto solution. It can be seen that the Pareto solutions are well-constrained and cover a small portion, suggesting that these parameters are well-identified by back analysis. Meanwhile, corresponding to the values of back-estimated parameters, K1, G1, K2, and G2, for the compromise solution Pc, are 1.0321 × 104, 2.455 × 103, 1.771 × 104, and 3.338 × 103 kPa, respectively.Figure 13
Normalized range of optimal parameter values for the soil parameters.


### 5.4. Calibrated and Predicted Responses at Impact Location A

The pore water pressure and crater depth of blow 1 using the back-analyzed soil parameters are simulated. Simulated and measured pore water pressures are plotted in Figure14. The coefficient of determination R2 is 0.87 and the RMSE is 3.29 kPa. It can be seen that the estimated pore water pressure is in good agreement with the measurement. The simulated crater depth for blow 1 is 0.8299 m which is almost the same as the measured 0.83 m in the field.Figure 14
Measured and simulated pore water pressure in depth at impact location A.The estimated soil parameters of the compromise solution are then used to predict the lateral displacement of ground at the same impact location as shown in Figure15. The overall trend of lateral displacements is adequately described. The predicted lateral displacements at depth greater than 5 m are very close to the measured field data. Therefore, the back-analyzed soil parameters obtained using the information of pore water pressure and crater depth can predict the lateral displacement of ground reasonably.Figure 15
Comparison of measured and predicted lateral displacement of ground at impact location A.A practice goal of back analysis is to predict the ground response to dynamic compaction based on the early stage observations and achieve a better understanding of soil behavior under dynamic compaction construction. In this section, the input parameters from the compromise solution of blow 1 are used to predict the compaction effects of dynamic compaction at subsequent blows (blows 2 to 5) at impact location A. Figure16 shows the observed and calculated crater depth with blow count. The crater depth can be generally well-predicted using the back-analyzed soil parameters for blows 2 and 3. The predicted deformations are slightly larger than the field observations. The differences between the predicted and measured value become large for blows 4 and 5.Figure 16
Comparison of field observed and calculated crater depth with blow count using compromise solution from 1st blow and 2nd blow at impact location A.To examine hardening effects of dynamic compaction in continuous blows at impact location A, the soil parameters are estimated based on measurements from blow 2. The back-estimated parameters,K1, G1, K2, and G2, are 1.01 × 105, 3.31 × 103, 1.503 × 104, and 1.43 × 103 kPa, respectively. Comparing with the back-analyzed parameters from blow 1, the elastic parameters of the first soil layer increase significantly while those of the second soil layer are only slightly increased. These results illustrate that hardening and strengthening of soil occurs due to continuous compaction. Figure 16 also presents the predicted crater depth based on the estimated parameters from blow 2. It can be seen that the prediction based on the second blow is more accurate than the first blow. Most of studies also show that the relative density of the soil increases significantly with the growth of blow numbers in dynamic compaction [11–15, 49, 50, 52].


### 5.5. Prediction at Impact Location B

In this section, the input parameters from the compromise solution of blow 1 at impact location A are used to predict the performance of the dynamic compaction for the first blow under the 2500 kN·m tamping energy at impact location B. The predicted crater depth is 1.02 m which is 8.51% larger than its observed value equal to 0.94 m. Predicted and measured pore water pressures are presented in Figure17. The RMSE is 7.618 kPa and the coefficient of determination R2 is 0.759. The predicted pore water pressure is almost the same with the observed results, except the depth that is close to 5 m has a little distinction, as shown in Figure 17. It is mainly because this depth is located at the interface of different soil layers, which significantly affected into the predicted value of pore water pressure, while effects of the relative friction between different soil layers are not taken into account by the calculation model. The estimated lateral displacements are compared with the field observation as shown in Figure 18. In addition, the RMSE is 5.399 mm and the coefficient of determination R2 is 0.991. The predicted lateral displacements agree well with the field observations. The soil properties before ground treatment at impact locations A and B are not quite different. Hence, the back-estimated soil parameters from the compromise solution of blow 1 at impact location A can be used to predict the response of the dynamic compaction for the first blow at impact location B. Therefore, the above results show that the compromise solution from biobjective back analysis can give generally acceptable predictions of dynamic compaction for an adjacent location.Figure 17
Comparison of measured and predicted pore water pressure at impact location B using compromise solution from back analysis of 1st blow at impact location A.Figure 18
Comparison of measured and predicted lateral displacement at impact location B using compromise solution from back analysis of 1st blow at impact location A.


## 5.1. Coupled Hydromechanical Response in Numerical Modeling

In this section, we present the FEM results including soil displacements and generated pore water pressure to illustrate the coupled hydromechanical responses due to dynamic compaction. The variation of lateral soil displacement with time is shown in Figure9. In this figure, the soil deformation near the edge of the hammer increases progressively, as the wave front propagates away from the hammer. In Figure 10, the generation of pore water pressure and its dissipation are illustrated. At the first stage of loading, during the beginning of the impact, the pore water pressure is increased instantaneously below the loading zone. As the wave progresses, the excess pore pressure moves downwards and decreases gradually below the loading zone.Figure 9
Contour lines of lateral displacement for the first blow: (a) 0.05 s, (b) 0.1 s, (c) 0.2 s, and (d) 0.3 s.
(a)
(b)
(c)
(d)Figure 10
Contour lines of pore water pressure for the first blow: (a) 0.05 s, (b) 0.1 s, (c) 0.2 s, and (d) 0.3 s.
(a)
(b)
(c)
(d)


## 5.2. Characteristics of the Pareto Front

Due to the page limit, we present only the result based on the measured data of blow 1 at impact location A. Figure11 shows the convergence of optimization with an increasing number of generations. The objective function f2 converges to less than 2 × 10−3mm at the 33th generation. The objective function f1 converges at 3.29 kPa at the 37th generation. The back analysis is terminated at the 68th generation.Figure 11
Evolution of the minimum objective function values as a function of the number of generations.Figure12 presents the trade-off of the f1-f2 biobjective Pareto front for the model. All the numerical simulated solutions are marked as gray square points. The Pareto solutions of the last generation are marked as blue circles. The Pareto front exhibits a rectangular shape, indicating that it is possible to minimize both objectives simultaneously. Here, the most informative Pareto point is selected from the Pareto solutions. This point is the knee point in the Pareto front, where a small improvement in one objective would lead to a large deterioration in at least one other objective [51]. Such a point can be regarded as the optimal solution of the multiobjective optimization problem and is more relevant to the decision maker than the entire Pareto front. In the graph, we denote the point as Pc and the solution is subsequently referred to as the compromise solution.Figure 12
Pareto front and optimal solutions.


## 5.3. Identification of Soil Parameters

Figure13 shows normalized parameter ranges corresponding to the Pareto solutions that appeared in previous section. The back-estimated parameters, K and G, of the first two soil layers are listed along the abscissa. The subscript i of each parameter corresponds to the ith soil layer. The ordinate presents the parameter values normalized to the prior bounds (as given in Table 2). Each blue line in the graph represents one Pareto solution. It can be seen that the Pareto solutions are well-constrained and cover a small portion, suggesting that these parameters are well-identified by back analysis. Meanwhile, corresponding to the values of back-estimated parameters, K1, G1, K2, and G2, for the compromise solution Pc, are 1.0321 × 104, 2.455 × 103, 1.771 × 104, and 3.338 × 103 kPa, respectively.Figure 13
Normalized range of optimal parameter values for the soil parameters.


## 5.4. Calibrated and Predicted Responses at Impact Location A

The pore water pressure and crater depth of blow 1 using the back-analyzed soil parameters are simulated. Simulated and measured pore water pressures are plotted in Figure14. The coefficient of determination R2 is 0.87 and the RMSE is 3.29 kPa. It can be seen that the estimated pore water pressure is in good agreement with the measurement. The simulated crater depth for blow 1 is 0.8299 m which is almost the same as the measured 0.83 m in the field.Figure 14
Measured and simulated pore water pressure in depth at impact location A.The estimated soil parameters of the compromise solution are then used to predict the lateral displacement of ground at the same impact location as shown in Figure15. The overall trend of lateral displacements is adequately described. The predicted lateral displacements at depth greater than 5 m are very close to the measured field data. Therefore, the back-analyzed soil parameters obtained using the information of pore water pressure and crater depth can predict the lateral displacement of ground reasonably.Figure 15
Comparison of measured and predicted lateral displacement of ground at impact location A.A practice goal of back analysis is to predict the ground response to dynamic compaction based on the early stage observations and achieve a better understanding of soil behavior under dynamic compaction construction. In this section, the input parameters from the compromise solution of blow 1 are used to predict the compaction effects of dynamic compaction at subsequent blows (blows 2 to 5) at impact location A. Figure16 shows the observed and calculated crater depth with blow count. The crater depth can be generally well-predicted using the back-analyzed soil parameters for blows 2 and 3. The predicted deformations are slightly larger than the field observations. The differences between the predicted and measured value become large for blows 4 and 5.Figure 16
Comparison of field observed and calculated crater depth with blow count using compromise solution from 1st blow and 2nd blow at impact location A.To examine hardening effects of dynamic compaction in continuous blows at impact location A, the soil parameters are estimated based on measurements from blow 2. The back-estimated parameters,K1, G1, K2, and G2, are 1.01 × 105, 3.31 × 103, 1.503 × 104, and 1.43 × 103 kPa, respectively. Comparing with the back-analyzed parameters from blow 1, the elastic parameters of the first soil layer increase significantly while those of the second soil layer are only slightly increased. These results illustrate that hardening and strengthening of soil occurs due to continuous compaction. Figure 16 also presents the predicted crater depth based on the estimated parameters from blow 2. It can be seen that the prediction based on the second blow is more accurate than the first blow. Most of studies also show that the relative density of the soil increases significantly with the growth of blow numbers in dynamic compaction [11–15, 49, 50, 52].


## 5.5. Prediction at Impact Location B

In this section, the input parameters from the compromise solution of blow 1 at impact location A are used to predict the performance of the dynamic compaction for the first blow under the 2500 kN·m tamping energy at impact location B. The predicted crater depth is 1.02 m which is 8.51% larger than its observed value equal to 0.94 m. Predicted and measured pore water pressures are presented in Figure17. The RMSE is 7.618 kPa and the coefficient of determination R2 is 0.759. The predicted pore water pressure is almost the same with the observed results, except the depth that is close to 5 m has a little distinction, as shown in Figure 17. It is mainly because this depth is located at the interface of different soil layers, which significantly affected into the predicted value of pore water pressure, while effects of the relative friction between different soil layers are not taken into account by the calculation model. The estimated lateral displacements are compared with the field observation as shown in Figure 18. In addition, the RMSE is 5.399 mm and the coefficient of determination R2 is 0.991. The predicted lateral displacements agree well with the field observations. The soil properties before ground treatment at impact locations A and B are not quite different. Hence, the back-estimated soil parameters from the compromise solution of blow 1 at impact location A can be used to predict the response of the dynamic compaction for the first blow at impact location B. Therefore, the above results show that the compromise solution from biobjective back analysis can give generally acceptable predictions of dynamic compaction for an adjacent location.Figure 17
Comparison of measured and predicted pore water pressure at impact location B using compromise solution from back analysis of 1st blow at impact location A.Figure 18
Comparison of measured and predicted lateral displacement at impact location B using compromise solution from back analysis of 1st blow at impact location A.


## 6. Limitations

Dynamic compaction in saturated soils involves complex physical process such as energy transition from the falling hammer to the soil body through the contact surface, wave propagation, dynamic consolidation, and so forth. In this study, the back analysis is carried out using a specified numerical model with certain assumptions about geometry, boundary conditions, and soil constitutive models. The back analysis can only minimize the discrepancy between the simulated and observed performance but is unable to overcome inherent limitations of the numerical model.In the illustrative example, we considered only two objective functions. However, the methodology of Pareto multiobjective back analysis is applicable to complex problems with multiple objective functions. In the back analysis for dynamic compaction, the measured lateral displacement can be also used as the additional objective function. The effect of different types of objective functions is out of the scope of this study. Further research study may be conducted on this aspect.


## 7. Conclusions

This paper presented a parameter estimation method based on a coupled hydromechanical model of dynamic compaction and the Pareto multiobjective optimization technique. The hydromechanical model of dynamic compaction is established in the FEM program LS-DYNA. The multiobjective optimization algorithm, Nondominated Sorted Genetic Algorithm (NSGA-IIa), is integrated with the numerical model to identify soil parameters using multiple sources of field data. A real field case of dynamic compaction treatment was presented to illustrate the methodology. The following conclusions can be made:(1)
For the demonstrated dynamic compaction case, thef1-f2 biobjective Pareto front exhibits a rectangular shape, illustrating that it is possible to minimize the RMSEs of both pore water pressure and crater depth simultaneously. In addition, the Pareto solutions are well-constrained and cover a small portion, suggesting that these parameters are well-identified by back analysis.(2)
The back-analyzed soil parameters obtained using the information of pore water pressure and crater depth can predict the lateral displacement of ground reasonably.(3)
The hardening and strengthening of soil occurs due to continuous compaction. For the predictions of the ground response of the dynamic compaction at continuous blows, the prediction based on the second blow is more accurate than the first blow.(4)
The back-analyzed soil parameters of the compromise solution from the biobjective back analysis can give generally acceptable predictions of dynamic compaction for an adjacent location.



---

*Source: 127878-2015-07-13.xml*