Numerical analysis of tailing dam with calibration based on genetic algorithm and geotechnical monitoring data

Prior to performing the optimization, it is required to specify the data to which numerical model should be fitted. In this work, it was assumed that measurement data are functions of time, for example changes of displacement in time (data from benchmarks) or changes of piezometric head in time (data from piezometers).

The numerical model should be created using all possible construction information, such as geological layers arrangement, history of loading or geotechnical tests. The size of the model largely determines the time needed for optimization. Selection of nodes or elements that correspond to actual sensors is necessary to perform the algorithm. The results needed from the model are only appropriate values (displacement, pressure) of the selected nodes or elements in text form. It is obligatory to use FEM software that allow to create such an output automatically, that is without running the postprocessor.

The genetic algorithm is based on population that consist of individuals. The individual represents the parameters of all optimized materials. As the Mohr-Coulomb model was adopted to numerical calculation, parameters for each material include: angle of internal friction, cohesion, Young‘s modulus, Poisson‘s ratio and hydraulic conductivity. Different sets of parameters can be selected for each material. The individual is a sequence of numbers, which determines all the parameters according to the following scheme: φ1,φ2,…,φn,c1,c2,…,cn,E1,E2,…,En,ν1,ν2,…,ν3,k1,k2,…,kn,\matrix{{{\varphi _1},{\varphi _2}, \ldots ,{\varphi _n},{c_1},{c_2}, \ldots ,{c_n},{E_1},{E_2}, \ldots ,{E_n},{\nu _1},{\nu _2}, \ldots ,} \cr {{\nu _3},{k_1},{k_2}, \ldots ,{k_n},} \cr } where: ϕ – angle of internal friction, c – cohesion, E – Young modulus, ν – Poisson ratio, k – hydraulic conductivity, n – number of optimized materials, and the subscript refers to material number.

The initial population of individuals is generated randomly. However, parameters are drawn from ranges specified by the designer. First reason for such an approach is the fact that some parameters cannot exceed certain values, for example, the Poisson’s ratio. Secondly, these ranges can be determined on the basis of geotechnical documentation, which can lead to a result consistent with engineering experience. Additionally, the calibration time can be reduced.

To perform the evaluation, numerical calculation must be carried out. The same task is calculated, but the optimized materials parameters are assigned according to the individuals. Evaluation of an individual involves calculation error defined as the difference between measurement data and numerical output. Since measurement data can have different physical sense and therefore different units and values, data should be normalized or relative error can be used. The latter was implemented, the error is based on mean relative error estimated on the whole dataset according to the following formula: ε=1∑i=1rni∑i=1r(∑j=1ni|xn,i,j−xm,i,j||xm,i,j|),\varepsilon = {1 \over {\sum\nolimits_{i = 1}^r {{n_i}} }}\sum\limits_{i = 1}^r {\left( {{{\sum\nolimits_{j = 1}^{{n_i}} {\left| {{x_{n,i,j}} - {x_{m,i,j}}} \right|} } \over {\left| {{x_{m,i,j}}} \right|}}} \right)} , where: r – number of gauges in calibration, n_i – number of measurements in i-th gauge, x_n –measurement value from numerical analysis, x_m –measurement value from monitoring.

Three individuals with the lowest error values are selected for reproduction. In the method (Whitley and Starkweather, 1990) only two best individuals had been selected, however, it was decided to increase this number due to the computing capabilities of computers, which allow to calculate several FEM models at the same time, while maintaining the crossover of the best individuals.

Crossover is the exchange of information between two individuals. The crossover point is selected, in relation to which the sequence of numbers is divided. Then the parts of individuals are swapped as shown in the following scheme: 101000∨110111xxyxxy∧yxxyyx→101000yxxyyxxxyxxy110111,\matrix{{101000 \vee 110111} \cr {xxyxxy \wedge yxxyyx} \cr } \to \matrix{{101000yxxyyx} \cr {xxyxxy110111} \cr } , where: ∨, ∧ – crossover point, 1, 0, x, y – individuals features.

As a result, two new individuals are created, but only one is randomly selected for the next population. Pairing of individuals for crossover occurs according to the following pattern: 1↔2,2↔3,3↔1,1 \leftrightarrow 2,2 \leftrightarrow 3,3 \leftrightarrow 1, where: 1, 2, 3 – three best ranked individuals.

In each generation, i.e. that is every next population, the crossover point is randomly selected for each crossover. It should be emphasized, that the way parameters are written as individuals is important, although the parameters are independent, several of them describe one material. Already presented encoding scheme enforces the exchange of parameters in all materials, however, the following pattern: φ1,c1,E1,ν1,k1,φ2,c2,…,kn,{\varphi _1},{c_1},{E_1},{\nu _1},{k_1},{\varphi _2},{c_2}, \ldots ,{k_n}, imposes replacement among the materials. Since it is highly time-consuming to establish which method deliver the best results selecting the encoding system is treated as an algorithm parameter.

Mutation introduces small changes to the newly created individuals. Maximum values of changes are defined for each type of parameter, for example for all cohesion values mutation parameter is the same. Mutation can be conducted as addition or multiplication. The second option is useful for hydraulic conductivity, as its values can differ by several orders of magnitude. Parameters which resulted from mutation cannot exceed the ranges defined for initial population. A certain number of parameters (defined as algorithm parameter) are randomly selected for mutation. Each parameter can be increased (addition or multiplication by a mutation parameter), decreased (subtraction or division) or remain unchanged due to the random choice.

Each subsequent generation consists of individuals from the previous generation without the three worst ones; in their place three individuals created by crossbreeding and mutation are added. In the first generation, FEM analysis must be carried out for all the individuals, while in each subsequent generation, for only three new ones.

In this chapter the results of calibration carried out for synthetic data are presented. A numerical model with some material parameters was created and calculated. Selected nodes were treated as sensors and displacement and piezometric head values over the time were recorded. The data prepared in this way were noised by adding values from normal distribution with mean equal to 0 and standard deviation adapted to values of the data. Figure 1 presents an example of the noised synthetic data.

Figure 1

The same model was used both to create the data and to perform calibration to ensure that a solution of optimization exists. Moreover, the impact of model inaccuracies was eliminated. The calculation problem concerned embankment of flotation reservoir. The height of embankment increased in time as the water and sedimentation particles filled the reservoir. Numerical model was created using FEM software ZSoil (Zimmermann et al. 2016) with assumption of plain strain. The first order EAS (enhanced assumed strain) (Simo and Rifai, 1990) quadrilateral elements were used. After 205 computational step the height of the dam was equal to 65 m, while the depth of subsoil was 265 m. In the figures 2 to 6 subsequent stages of construction history, geometry, mesh and model's dimensions are presented.

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

On both right and left edges the horizontal displacement was blocked, while on the bottom edge horizontal and vertical displacement was fixed. On the left edge, water pressure was raised to simulate reservoir filling, and the right edge water pressure had constant value throughout the whole calculation. As it was already mentioned elastic-perfectly plastic Mohr-Coulomb with unassociated plastic flow rule was assumed. Coupled analysis, i.e. that is consolidation, with adoption of van Genuchten's model of unsaturated zone (Van Genuchten, 1980) was selected to calculate the distribution of pore pressure and deformation of the structure. The model consisted of almost 12,000 elements and the calculation took less than 30 minutes (on 3.2 GHz CPU and 32GB RAM computer). Nodes selected as sensors are shown in Figure 7. Sensors with letter U are benchmarks, which gather information about vertical and horizontal displacement, while sensors starting with P are piezometers and these collect piezometric head values.

Figure 7

Initial population consisted of 15 individuals and 15 generations were assumed. Materials selected for optimization were: 1A, 2B, 3C, 4D, 5E. Their location in the model is shown in Figures 2 to 6. It was decided that 10 parameters could be mutated. Ranges used for the assessment of initial population as well as mutation parameters are presented in Table 1. Such a calibration took about 10 hours, of which only about 5 minutes were needed for the operations related to the performing genetic algorithm (not FEM analysis).

The calibration was performed 5 times for different random initial populations, yet for the same initial ranges and mutation values. Parameter values selected to create synthetic data and optimized material parameters are presented in Table 1. Optimizations are described with Roman numbers I–V. It has to be mentioned that exact parameters have not been reproduced in any of the 5 optimizations. However good approximation was achieved for materials 2B, 3C and 5E. Surely, it is connected to role the played by each of these material in the model. Layer 2B corresponds to about 80% of elements in subsoil, layer 3C is the weak layer that enforces the shape of potential failure zone and layer 5E is the stiffest one. The hydraulic conductivity values have been optimized satisfactorily. Probably narrowing the ranges of Poisson ratio could improve results, as these ranges were the same for all materials and relatively wide.

Figure 8 shows the change of error depending on the generation number, that is how the error in each of the five calibrations was reduced. It can be observed that in all five cases error decreased, which means that the process gradually improves the initial individuals, which were created randomly. However, there is a possibility that there will be no significant improvement during the optimization process if the first population contains very good individual, similar to calibration IV. From a practical point of view, such a situation cannot be treated as proof of algorithm failure because the parameters provide a reasonable fitting to monitoring data.

Figure 8

Figures 9 to 13 show the fitting diagrams for all the sensors. Data from measurement is marked with black color and in these plots they are shown without noise. On the X-axis the time as number of computational step is presented and on the Y-axis sensor quantity: displacement or piezometric head. It is difficult to decide which calibration (from I to V) is the best, but all of them should be acceptable in real case. It can be observed that for different sensors different sets of parameters correspond to the best fitting obtained.

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13