Next Article in Journal
A Rule-Based System for Human Performance Evaluation: A Case Study
Next Article in Special Issue
The Modified Void Nucleation and Growth Model (MNAG) for Damage Evolution in BCC Ta
Previous Article in Journal
Evolving a Multi-Classifier System for Multi-Pitch Estimation of Piano Music and Beyond: An Application of Cartesian Genetic Programming
Previous Article in Special Issue
Benchmarking Numerical Methods for Impact and Cratering Applications
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Assimilation of Dynamic Combined Finite Discrete Element Methods Using the Ensemble Kalman Filter

by
Humberto C. Godinez
*,†,‡ and
Esteban Rougier
Los Alamos National Laboratory, Los Alamos, NM 87545, USA
*
Author to whom correspondence should be addressed.
Current address: Applied Mathematics and Plasma Physics, MSB284, Los Alamos, NM 87544, USA.
These authors contributed equally to this work.
Appl. Sci. 2021, 11(7), 2898; https://doi.org/10.3390/app11072898
Submission received: 22 February 2021 / Revised: 18 March 2021 / Accepted: 19 March 2021 / Published: 24 March 2021
(This article belongs to the Special Issue Fracture Mechanics – Theory, Modeling and Applications)

Abstract

:
Simulation of fracture initiation, propagation, and arrest is a problem of interest for many applications in the scientific community. There are a number of numerical methods used for this purpose, and among the most widely accepted is the combined finite-discrete element method (FDEM). To model fracture with FDEM, material behavior is described by specifying a combination of elastic properties, strengths (in the normal and tangential directions), and energy dissipated in failure modes I and II, which are modeled by incorporating a parameterized softening curve defining a post-peak stress-displacement relationship unique to each material. In this work, we implement a data assimilation method to estimate key model parameter values with the objective of improving the calibration processes for FDEM fracture simulations. Specifically, we implement the ensemble Kalman filter assimilation method to the Hybrid Optimization Software Suite (HOSS), a FDEM-based code which was developed for the simulation of fracture and fragmentation behavior. We present a set of assimilation experiments to match the numerical results obtained for a Split Hopkinson Pressure Bar (SHPB) model with experimental observations for granite. We achieved this by calibrating a subset of model parameters. The results show a steady convergence of the assimilated parameter values towards observed time/stress curves from the SHPB observations. In particular, both tensile and shear strengths seem to be converging faster than the other parameters considered.

1. Introduction

The accurate simulation of fracture initiation, propagation, and arrest in materials is an important problem for a number of applications, ranging from the study of micro-scale behavior of materials to the analysis of large-scale fracture initiation and propagation. Many numerical solutions have been proposed to simulate fractures; one of the most widely accepted is the Combined Finite-Discrete Element Method (FDEM) [1,2]. FDEM, which merges solutions developed under both finite elements and discrete element methods, is mainly aimed at modelling failing, fracturing, and/or fragmenting solids. This is achieved by the usage of Finite Element formulations to resolve stresses in the material combined with Discrete Element formulations to resolve fracture and contact problems [2]. A review of the current fracture models based on FDEM is presented by [3]. Furthermore, FDEM has been widely used in real-world problems, such as dam and rock slope stability [4,5], hydraulic fracture and fracture propagation on shale [6,7,8], fracture coalescence [9,10], rock mass strength [11], dynamic fault rupture [12,13,14], and granular gouge response [15,16].
As with many models, several material-specific parameters (i.e., material properties) are needed for a FDEM simulation, such as density, Young’s modulus, Poisson’s ratio, and strength parameters for failure modes I (normal) and II (tangential). Additionally, the specific fracture energies for these two modes are also needed. Closely related to the specific fracture energies is the use (within FDEM) of a displacement-softening curve to represent the strength of degradation as micro-fractures form in the material. In previous work, the sensitivity of the FDEM model to perturbations in these parameters was analyzed [17,18], where the authors found a relatively small number of parameters influencing the behaviour of the stress versus time curve obtained from FDEM.
To accurately simulate fracture behavior, the characterization of displacement-softening parameters need to be correctly specified. However, experimental information that could be used to determine these parameters is very often not easy to find, and modelers are forced to establish them in an ad hoc manner or based on their own experience. Data assimilation is a method that fuses observational data into models to improve the nowcast and/or forecast of the model [19,20]. There are several types of data assimilation methods, from statistical based on Bayesian theory, to optimization-based methods. Data assimilation methods are used in several geophysical sciences, including climate, atmospheric weather prediction, space weather, and ocean sciences, to name a few. Additionally, assimilation methods are increasingly being implemented in geomaterials and fracture propagation to improve the simulation of geophysical models (weather, climate, subsurface flow, etc.) [19,20,21].
The calibration of model parameters is of relevance for a number of applications. For example, it is relevant to the simulation of complex masonry construction problems [22,23] where non-linear analysis of the response of the structures are needed. In these cases, the problem itself is complex enough, and the uncertainty on the parameter determination can be quite high. This is also the case for the simulation of underground processes, such as hydraulic fracturing and earthquake ruptures, where knowledge about material properties is limited to areas where core sampling can be extracted and analyzed. This provides a starting point for analysts to populate their models, but in many cases this is not sufficient, and further calibration is needed. Sensitivity analysis [17] and assimilation methods [18] provide a robust and systematic way of analyzing the system under consideration as a whole, while identifying the key parameters and establishing a confidence interval for them.
The Split Hopkinson Pressure Bar (SHPB) is a standard laboratory experiment where the tensile strength of geologic materials, measured at high rates of strain, is obtained. It is worth noting that the simulation of the SHPB experiment is very demanding for the numerical method being used, mainly in the post-peak portion of the stress-time curve. The correct representation of damage after the material has reached its peak poses many numerical challenges, for example with the proper accounting of the energy being dissipated. In this work, it will be demonstrated that via this assimilation process, the Combined Finite Discrete Element Method is able to simulate the post-peak behavior of granite at high strain rates. Furthermore, the calibrated properties obtained from this assimilation exercise can be used to populate more complicated models, such as when trying to simulate underground mining explosions, or dynamic earthquake rupture problems. This paper presents the procedure to obtain such calibrated strength parameter sets, but the application to other types of problems is outside its scope.
This work presents the implementation of the ensemble Kalman filter to estimate key model parameters of a dynamic Combined Finite Discrete Element model, which simulates the evolution of fractures and cracks in different geomaterials. Specifically, the ensemble Kalman filter (EnKF) is applied to the Hybrid Optimization Software Suite (HOSS) to better simulate a set of laboratory experiments. The main objective is to estimate a set of key strength parameters for the model which have the most influence on the material response to the imparted load. For this, the FDEM numerical results are contrasted against observations obtained for a SHPB experiment on granite samples.
The paper is organized as follows: Section 2 provides a brief introduction of the HOSS model and the parameters to be estimated; Section 3 presents the SHPB experiment description, as well as the HOSS model setup to simulate an SHPB experiment; Section 4 provides a brief overview of the ensemble Kalman filter and the assimilation algorithm, and the results of the assimilation are presented in Section 5; and Section 6 provides the conclusions and discussion of the results.

2. Description of the Model

In this work, the 2D version of the Hybrid Optimization Software Suite (HOSS) [24,25,26,27], which is an in-house implementation of the Combined Finite-Discrete Element method (FDEM) [2,28,29], is used. For an in-depth description of FDEM, the following references are recommended: [2,3,5,17,30]. In a previous work [17], the authors performed a detailed sensitivity analysis using the Fourier Amplitude Sensitivity Test (FAST [31,32]), where seven HOSS parameters were explored. These parameters include tensile and shear strengths, σ n m a x and σ t m a x , the maximum relative displacements in the normal and tangential directions, δ n m a x and δ t m a x , and the shape of the displacement-softening curve (see [17]). The sensitivity analysis study showed that only two parameters had an impact on the stress versus time curve (for a point located at the center of the sample) obtained from HOSS for this type of problem. These parameters are the shear strength σ t m a x and the maximum relative displacements in the tangential direction δ t m a x . The sensitivity results are consistent with a previous work performed on the HOSS model by Osthus et al. [18], where only the σ t m a x is found to be of significant influence on the model. A combination of these parameters defines the specific fracture energy dissipated in mode II, E t , which is given by
E t = σ t max δ t max A 0 ,
where A 0 is given by A 0 = 0 1 z ( x ) d x (see [17]). The question remains: How do we select all these coefficients or material parameters, when there is not enough experimentally-based information, in a way that the model response is both physically sound and in agreement with available experimental data? As a way of circumventing the information limitation mentioned previously, in this work, an assimilation estimation of these parameters is proposed. In order to perform a comprehensive assimilation estimation of the HOSS model to the input parameters, both the shear strength σ t m a x and the maximum relative displacement in tangential direction δ t m a x are estimated through data assimilation. It is worth noticing that both parameters σ t m a x and δ t m a x are independent from each other. The way in which the δ t m a x is perturbed in this work is by changing the energy in tangential direction and therefore obtaining the maximum displacement using Equation (1). A more detailed description of the HOSS model and its parameters is presented in Appendix A and in [17].

3. SHPB Experiment and HOSS Model Setup

The SHPB experiment is used to indirectly determine the tensile strength of materials at high rates of strain. In this type of experimental setup, a cylindrical sample is positioned between two slender metal bars (see Figure 1). The load is generated by the impact of a projectile on the free end of the incident bar. This impact creates a stress wave that travels along the incident bar, and after a certain amount of time it reaches the sample, imposing a compressive load on it. A portion of the wave is transmitted to the sample, while the rest is reflected back into the incident bar [33,34]. Similar wave reflection behavior is observed between the sample and the transmission bar. In order to measure the amplitudes of the wave trains, strain gauges are placed in both bars. This allows for the determination of the initial load in the form of a velocity versus time series, and also of the contact force between the bars and the sample, f ( t ) , which in turn allows for the calculation of the stress at the middle of the sample as follows [30]:
σ n = 2 f ( t ) π d h ,
where h and d are the sample’s thickness and the diameter, respectively.
Observational data from the SHPB experiment described above is obtained for granite. The elastic properties and the density of the material are listed in Table 1. More details about the setups and the experiment in general are provided in [30,35]. For the simulation of the SHPB experiment with HOSS, the bars are discretized using around 5000 constant-strain triangular finite elements. For the granite sample in the middle of the experiment (circle in Figure 1), a mesh is generated using 17,500 finite elements with an average size of 0.5 mm. The action of the impactor was simulated by enforcing a velocity boundary condition on the x direction, v ( t ) , on the nodes that are located on the incident bar’s free end (left bar on Figure 1), as shown in Figure 2.

4. Data Assimilation

The simulation of complex physical phenomena is commonplace in many areas of science. As our knowledge of the physical phenomena increases, more complete models can be developed. Nevertheless, model errors and bias, resulting from uncertain parameters, incorrect initial and/or boundary conditions, and unaccounted physical processes have a significant influence on the accuracy of model forecasts. In this regard, the use of data assimilation methods, which combine data collected from the model, observations, and corresponding error statistics, is a viable method for estimating improved model forecasts [19,20]. This is accomplished by modifying the model parameters, as well as the model variables, in order to better agree with observational data using the Bayesian probability theory. In this way, data assimilation provides an improved model forecast, along with its probability distribution function to evaluate how trustworthy the resulting forecast is. Depending on the model, its dimension, and its complexity, there are different data assimilation methods that can be applied. The methods range from deterministic to Monte-Carlo-based methods which take advantage of the model dynamics. There are several assimilation methods, and the most widely used are the variational methods (3D-Var and 4D-Var) [36], Kalman filtering [37], as well as the Monte Carlo techniques [38,39,40,41].
In this work, the ensemble Kalman filter has been selected to assimilate data into the HOSS model. The main motivation for selecting an EnKF method to calibrate the model parameters in HOSS was because the EnKF is one of the most computationally inexpensive and easy to implement [42]. There are other methods that can potentially be used for the model calibration, such as control theory. Unfortunately, due to the complexity of the model and its simulation requirements, these methods might be hard to implement. The EnKF has been shown to be computationally feasible and efficient for parameter estimation problems, since it is a compromise between Monte Carlo methods and Bayesian theory for Gaussian distributions [42]. This compromise enables the EnKF to use statistical information efficiently and provide an estimate of the probability distribution function of the model parameter. Practically, the computational efficiency of the EnKF translates to the number of model simulations needed to achieve convergence. A typical number of model simulation that the EnKF needs is around 30, whereas other methods might need significantly more simulations, without even providing a probability distribution of the model parameters. This advantage makes the EnKF an attractive alternative to estimate the parameters for a complex model like HOSS. A description of the EnKF method is presented in the following subsection.

Ensemble Kalman Filter

The ensemble Kalman Filter (EnKF) [41,43] is an approximation to the Kalman filter [37] designed for non-linear models. The EnKF is widely implemented in several fields of science for its effectiveness in producing good results and the overall ease of implementation. The model mean and error covariance matrix is approximated through an ensemble of simulation results, which are used in the EnKF formulation. The observations are used to update the ensemble, through the Kalman filter formulas, and then evolved through the model between assimilation cycles. Details on the EnKF method and algorithm can be found in [43,44].
The EnKF algorithm analysis equations are performed in the following way: For a vector of m observations y o and an ensemble of N model simulations valid at time t k , the EnKF analysis equations are given by:
x i a = x i f + K y i o H x i f , i = 1 , , N ,
K = P f H T H P f H T + R 1 ,
where K R n × m is the Kalman gain matrix; x i a R n is the analysis; P f R n × n is the model covariance matrix; R R m × m is the covariance matrix of the observations; H R m × n is the observation operator used to map the model state to the observations; and y i o R m is the observation vector with Gaussian white noise.
The EnKF approximates the model covariance matrix using the ensemble simulation. There are a number of problems with this approximation, namely the introduction of artificial cross-correlations between model state variables, which can introduce errors when performing assimilation with the Kalman filter equations [45]. To solve this particular problem, the correlation matrix of the model ( P f ) is modified by removing the cross-correlation between parameters and variables that we know are not correlated. That is, the off-diagonal terms in P f that correspond to correlation between parameters and variables that we know are not correlated is set to zero. This is called localization of the covariance matrix.
It is also well-known that using the EnKF to estimate model parameters can lead to ensemble collapse, that is, when the parameter estimate in all ensemble members is essentially the same. In order to avoid ensemble collapse in the model parameters, an inflation technique is implemented. The inflation, as the name suggests, increases the difference between the ensemble parameter average and each individual ensemble parameter value. That is, the parameters are inflated as
p ^ i = p ¯ + λ p i p ¯ ,
where p i is the HOSS parameter value from ensemble i, p ¯ is the average, p ^ i is the inflated parameter for ensemble i, and the inflation factor is denoted by λ . For our problem, λ is set to approximately 45 % of the standard deviation for the initial parameter spread.
For the HOSS model, the parameters considered cannot be changed mid-simulation. That is, the HOSS parameters are set at the start of the simulation, and once the simulation is underway, they cannot be modified. Due to this, the EnKF is implemented as a smoother [46], which considers the results of HOSS over time as a single state vector to be assimilated using the observations over time. Additionally, the EnKF is implemented in an iterative way to better approximate the HOSS parameter values. The full mathematical framework for the iterative EnKF can be found in [47,48,49]. EnKF algorithm implementation is done in the following steps:
  • Initialization: An ensemble of parameter values p i , i = 1 , , N is provided by the latin hypercube sampling strategy, where p i = σ t , i max , E t , i ;
  • Simulation: For each ensemble member p i simulate HOSS for the SHPB experiment, where the output is a stress versus time curve ( x i for i = 1 , , N );
  • Assimilation: Using the HOSS output from the ensemble x i for i = 1 , , N and observations y o from the SHPB experiment, assimilate using the Kalman analysis Equations (3) and (4), where this will provide updated parameter values p i a = σ t , i max , a , E t , i a for all ensemble members;
  • Cycle: Feedback the updated parameter values p i a to HOSS and repeat Steps 2 and 3 until we see convergence of the error between the HOSS output and SHPB experiment ( x ¯ y o < t o l , where x ¯ is the ensemble average);
  • Forecast: Once convergence is achieved, perform a “forecast” of the HOSS model (best simulation with estimated parameter values).
The workflow is illustrated in Figure 3, where there are the five basic steps as in the description above.

5. Assimilation Results

5.1. Data Assimilation Simulation Setup

For the EnKF experiments, the ensemble is generated by perturbing/assimilating each of the HOSS parameters of interest at the beginning of the simulation. The parameters of interest are the shear strength σ t m a x and the maximum relative displacement in tangential direction δ t m a x . Due to the set-up of HOSS, we cannot directly modify the maximum relative displacement, and instead, for the assimilation experiments, we use the relationship given by Equation (1) and modify the specific energy in the tangential direction E t . That is, the two parameters that are perturbed are σ t m a x and E t . The parameters are described in Table 2, where the interval indicates the range of values for each parameter, and are sampled using the Latin Hypercube Sampling method (LHS) [50] to generate the initial HOSS ensemble simulations. The initial distribution for each of the parameters is a uniform distribution with lower and upper limits provided in Table 2. The HOSS simulation uses the parameter values for σ t m a x and E t , given in Table 2, to define the properties of the HOSS model. An ensemble of 25 HOSS simulations is used for all EnKF assimilation experiments presented. Each HOSS simulation requires 31 CPU cores for approximately 40–50 min of CPU time; hence, each HOSS simulation requires 1240–1550 min in CPU time (number of cores × CPU time requested for all cores). In total, the ensemble then requires between 31,000 and 38,750 min in CPU time (number of cores × CPU time per core × number of HOSS simulations). The HOSS simulation was performed on a Intel Xeon Broadwell CPU on the Grizzly high-performance computing platform available at the Los Alamos National Laboratory.

5.2. Data Assimilation Validation Experiments

The data assimilation algorithm, described in Figure 3, is validated using model-generated synthetic observations. This is performed with a HOSS setup to simulate the SHPB experiment described in Section 3. The synthetic observations experiment, referred to as a twin experiment, consists of selecting a set of “reference” parameter values for σ t m a x and E t , from which a HOSS model simulation for the SHPB is generated and used to get synthetic observations. In this case, the observations are the stress versus time curve generated from this reference HOSS simulation. Afterwards, an EnKF assimilation experiment is performed using these SHPB synthetic observations. To validate the EnKF assimilation method if the estimated parameters provided by the EnKF are close to the “reference” parameters chosen to generate the synthetic data, then the assimilation method is validated, since it provided the correct parameter values for this synthetic observation experiment.
The parameter values given by the EnKF validation experiment show steady convergence from the very first iteration, and provides a very good approximation at three iterations (see Figure 4). This steady convergence proves that the EnKF can accurately estimate the correct parameter values for σ t m a x and E t for the current setup. In order to ensure that the EnKF is actually converging to the correct solution, we compared the synthetic observed SHPB with the EnKF HOSS ensemble. As expected, the HOSS ensemble time-stress curves converge to the synthetic SHPB time-stress curve, as shown in Figure 5.
To verify that the data assimilation performed consistently over different parameter values, an additional two experiments with different reference parameter values were performed. The assimilation results on those two additional experiments are consistent with the results presented in this section—that is, both parameters σ t m a x and E t converged to the reference parameter values in no more than three EnKF iterations.

5.3. Results with SHPB Observations

The EnKF assimilation estimates the most likely mean and covariance of a normal distribution for each parameter that produces the HOSS SHPB simulation that best matches the data. The assimilation experiment setup is the same as in Section 5.1, where we are estimating σ t m a x and E t .
From our twin experiments shown in Section 5.2, three EnKF assimilation iterations are needed to estimate a fairly accurate parameter value for σ t m a x and E t . The observations are composed of a series of stress versus time points obtained from a SHPB laboratory experiment. The results show a good convergence for σ t m a x , where the assimilation provides a parameter value of around 0.15 (see Figure 6). For σ t m a x , the assimilation is able to hone in to the most likely value. This might not be the case for E t , where the assimilation provides a monotonically decreasing parameter value. Given that no strong suggestion of convergence is seen in this parameter value, it is difficult to conclude that the final value provided in the third iteration of the assimilation is the most likely parameter value for this SHPB observations.
To provide a clearer assessment of the convergence of the assimilation, the time–stress curve for all ensemble members, from the initial estimation to the second iteration, is plotted against the observed SHPB stress versus time curve. Figure 7 shows the time–stress curve for the initial ensemble with uniform distribution (blue lines), the curves for the first EnKF assimilation iteration (green curves), and the second assimilation iteration (red lines). As can be seen, the curves are converging towards the observed SHPB curves, but with the risk of ensemble collapse, it is more prudent to avoid over-iterating the EnKF.

6. Conclusions

Accurate simulation of fracture propagation for a number of materials is an important problem in the engineering community in general. There are a number of state-of-the-art numerical models in the community for fracture simulations, and among them is the HOSS model. In the HOSS model, as well as in other models, a number of important processes are determined by the value of key model parameters, such as material softening or degradation behavior, and energy dissipation and strength in both normal and tangential directions. Oftentimes, the appropriate value of these parameters becomes elusive for a number of materials, which poses a difficulty for the accurate simulation of fracture behavior in these materials. In this work, we present a systematic and reliable framework to estimate two key model parameters for HOSS using data assimilation for a Split Hopkinson Pressure Bar experiment for granite material. This was achieved via the implementation of an ensemble Kalman filter assimilation method. The main advantage of the EnKF is its ease of implementation and reliability to provide an estimate of the parameter values for the SHPB observational data.
Given that the parameters are set only at the beginning of the HOSS simulation, and have an impact throughout the simulation, the assimilation approach is based on the comparison of the stress versus time curve from HOSS against the one obtained from SHPB experiments. In essence, the EnKF assimilation utilizes this stress versus time curve as a state vector, computes the correlation between the stress versus time curve and the parameters, and exploits that correlation to adjust the values of the parameters so as to agree with the SHPB stress versus time curve. The EnKF assimilation, as explained, was applied three times (called, in this work, EnKF iterations) to obtain the parameter values. A previous sensitivity analysis test [17] indicated that the two main parameters that play a significant role in determining the stress versus time curve in the HOSS experiment are the strength and energy parameters in the tangential direction ( σ t and E t ). Therefore, only these two parameters are estimated using the EnKF assimilation method. To determine the accuracy and reliability of the assimilation method, we performed a twin experiment test. In this twin experiment, we generated synthetic data by sampling from a stress versus time curve from a HOSS simulation, where we fix the parameter values to a chosen reference value. Afterwards, we assimilate the synthetic observations to determine whether the assimilation method can recover the reference values used to generate the observations. Our twin experiment indicates that the assimilation can, in fact, recover the reference values (see Section 5.2).
The stress versus time curve from a SHPB experiment is used to estimate the model parameters σ t and E t using EnKF data assimilation. We observed that after three EnKF iterations, there is convergence in σ t , but there is no such convergence on E t (see Figure 6). The lack of convergence in E t might be attributed to the need of more EnKF iterations, or, alternatively, that the EnKF might not be able to correctly estimate this parameter. Nevertheless, the resulting model parameter estimates provide a HOSS simulation which show good convergence to the SHPB observational data (see Figure 7). Since the parameters are estimated through the assimilation, the HOSS convergence in the stress versus time curve indicates that the EnKF is adjusting the parameter’s values towards the best HOSS simulation for the observations considered.
Future work is still to be done for the parameter estimation of HOSS, specifically the exploration of other model parameters, as well as different observation types for calibration. One of the advantages of the EnKF data assimilation method is that it can handle a moderate increase in the number of model parameters to estimate without changing the algorithm outlined in the paper and with no considerable increase in computational cost [42]. In this paper, we only considered stress versus time curves, which might provide enough information to estimate certain HOSS parameters, but other parameters, such as the softening curve parameters [17] were not analyzed. Additionally, multiple sets of SHPB time–stress curves can be considered to further constrain the possible parameter values, and provide a template for the parameter values that are possible for SHPB experiments for certain materials.

Author Contributions

All authors contributed equally to the performance of the research and writting of the paper. All authors have read and agreed to the published version of the manuscript.

Funding

The funding for the project came from the Laboratory Directed Research and Development program from Los Alamos National Laboratory, LDRD # 20170103DR.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data used in this research can be shared from the authors upon request.

Acknowledgments

This research was performed as part of the project titled “Advancing Predictive Capability for Brittle Failure Using Dynamic Graphs project”. Authors would like to thank the LANL Institutional Computing program for their support in generating data used in this work. The authors would also like to thank three anonymous reviewers whose comments help improve the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Brief FDEM Overview and HOSS Model Parameters

The combined finite–discrete element method was conceived to represent the transition from continuum to discontinuum when a material breaks or fractures. The method formulation is based on the assumption that a material is composed of a group of deformable discrete entities, called discrete elements, that interact with each other via cohesive and contact forces. The governing equation that FDEM resolves, in its more general form, is given by
M 2 x t 2 + C x t = f ,
where M and C are the lumped mass and the damping matrices, respectively, while x is the vector of displacements and f is the vector of forces acting on the nodes of the finite element mesh. Vector f groups the forces due to material deformability, contact interaction, and cohesion with other elements, boundary conditions, and so forth. FDEM implementations use explicit time integration schemes [28,51] to resolve Equation (A1). Therefore, there is no need to assemble global matrices for the whole system. In turn, Equation (A1) is resolved on an element-by-element basis.
In order to accurately simulate fracture processes, HOSS uses a number of parameterizations in its hybrid finite-discrete element scheme. The parameterizations include a softening curve to simulate bond breakage (i.e., fracture) behavior between finite elements’ interfaces in tangential and normal directions. This is achieved via a combined single and smeared crack model [52] implemented at the finite element’s interfaces. There are a number of key model parameters that play a significant role in this crack model. The parameters we focus on in this work include tensile and shear strength ( σ n m a x and σ t m a x ), maximum relative displacements in the normal and tangential direction ( δ n m a x and δ t m a x ), and the parameters that determine the shape of the softening curve ( a , b , c ). The specific fracture energy dissipated at the interface during fracturing is calculated as
G = 0 1 z x σ m a x δ m a x d x ,
where
z x = 1 a + b 1 a + b exp x a + c b a + b 1 a b a 1 x + b 1 x c
is the functional form for the displacement-softening evolution [52], and
x = 0 if δ δ e 1 if δ > δ e + δ m a x δ δ e δ m a x if δ e < δ δ e + δ m a x .
The displacement-softening curve is defined by three parameters (a, b, and c), see Equation (A3). In the current work, the displacement-softening curve is utilized to describe both Mode I and Mode II fracture behavior. A description of the variables listed in Equations (A2)–(A4) can be found in [18].

Model Parameters

The specific energy in Mode I is given by
E n = σ n max δ n max A 0 ,
where A 0 is
A 0 = 0 1 z ( x ) d x .
Given that the shape of the displacement-softening curve is hard to obtain from laboratory experiments, it is often derived in an ad hoc manner. The expression for the specific energy for Mode II is given by
E t = σ t max δ t max A 0 .
In order to perform a comprehensive assimilation estimation of the HOSS model, the parameter set is completed by also considering the strength of the material in shear, σ t max , and the corresponding specific energy, E t , (Equation (A7)). Given these, the maximum displacements (i.e., in opening and in sliding modes) can be calculated as follows:
δ n max = E n σ n max A 0 ,
δ t max = E t σ t max A 0 .
It is worth noting that the normal and tangential modes of failure are linked with each other. This means that if an interface breaks in tension (i.e., normal mode), then it also loses its cohesion in the shear (i.e., tangential) direction and vice versa. The way this is implemented into HOSS is via a connection between δ t / δ t m a x and δ n / δ n m a x , defined as shown in Equation (Figure A1). In this way, the damage metric, which affects the interfaces’ strength degradation, is obtained as a combination of the post-peak displacements in both normal and tangential directions.
Figure A1. Connection between the normal and tangential failure modes, as implemented in HOSS.
Figure A1. Connection between the normal and tangential failure modes, as implemented in HOSS.
Applsci 11 02898 g0a1
For more information on these parameters, and the FDEM model used in this study, we refer the reader to the works from [2,24,28,29].

References

  1. Munjiza, A. Discrete Elements in Transient Dynamics of Fractured Media. Ph.D. Thesis, Swansea University, Swansea, UK, 1992. [Google Scholar]
  2. Munjiza, A. The Combined Finite-Discrete Element Method; John Wiley and Sons Ltd.: New York, NY, USA, 2004. [Google Scholar]
  3. Lisjak, A.; Grasselli, G. A review of discrete modeling techniques for fracturing processes in discontinuous rock masses. J. Rock Mech. Geotech. Eng. 2014, 6, 301–314. [Google Scholar] [CrossRef] [Green Version]
  4. Tatone, B.; Lisjak, A.; Mahabadi, O.; Grasselli, G.; Donnelly, C. A preliminary evaluation of the combined finite element-discrete element method as a tool to assess gravity dam stability. In Proceedings of the CDA 2010 Annual Conference, Niagara Falls, ON, Canada, 2–7 October 2010. [Google Scholar]
  5. Lisjak, A.; Grasselli, G. Combined finite-discrete element analysis of rock slope stability under dynamic loading. In Proceedings of the Pan-Am CGS Geotechnical Conference, Toronto, ON, Canada, 2–6 October 2011. [Google Scholar]
  6. Carey, J.W.; Lei, Z.; Rougier, E.; Mori, H.; Viswanathan, H. Fracture-permeability behavior of shale. J. Unconv. Oil Gas Resour. 2015, 11, 27–43. [Google Scholar] [CrossRef]
  7. Hyman, J.D.; Jiménez-Martínez, J.; Viswanathan, H.S.; Carey, J.W.; Porter, M.L.; Rougier, E.; Karra, S.; Kang, Q.; Frash, L.; Chen, L.; et al. Understanding hydraulic fracturing: A multi-scale problem. Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 2016, 374. [Google Scholar] [CrossRef] [PubMed]
  8. Rougier, E.; Knight, E.; Munjiza, A. Fluid Driven Rock Deformation via the Combined FEM/DEM Methodology. In Proceedings of the 46th US Rock Mechanics/Geomechanics Symposium, Chicago, IL, USA, 24–27 June 2012; American Rock Mechanics Association: Alexandria, VA, USA, 2012. [Google Scholar]
  9. Moore, B.A.; Rougier, E.; O’Malley, D.; Srinivasan, G.; Hunter, A.; Viswanathan, H. Predictive modeling of dynamic fracture growth in brittle materials with machine learning. Comput. Mater. Sci. 2018, 148, 46–53. [Google Scholar] [CrossRef]
  10. Euser, B.; Rougier, E.; Lei, Z.; Knight, E.; Frash, L.; Carey, J.; Viswanathan, H.; Munjiza, A. Simulation of fracture coalescence in granite via the combined finite-discrete element method. Rock Mech. Rock Eng. 2019, 52, 3213–3227. [Google Scholar] [CrossRef] [Green Version]
  11. Elmo, D.; Stead, D. An integrated numerical modelling–discrete fracture network approach applied to the characterisation of rock mass strength of naturally fractured pillars. Rock Mech. Rock Eng. 2010, 43, 3–19. [Google Scholar] [CrossRef]
  12. Okubo, K.; Bhat, H.S.; Rougier, E.; Lei, Z.; Knight, E.E.; Klinger, Y. Dynamic fracture network around faults: Implications for earthquake ruptures, ground motion and energy budget. In Proceedings of the AGU Fall Meeting Abstracts, New Orleans, LA, USA, 11–15 December 2017. [Google Scholar]
  13. Okubo, K.; Bhat, H.S.; Rougier, E.; Marty, S.; Schubnel, A.; Lei, Z.; Knight, E.E.; Klinger, Y. Dynamics, Radiation, and Overall Energy Budget of Earthquake Rupture With Coseismic Off-Fault Damage. J. Geophys. Res. Solid Earth 2019, 124, 11771–11801. [Google Scholar] [CrossRef] [Green Version]
  14. Klinger, Y.; Okubo, K.; Vallage, A.; Champenois, J.; Delorme, A.; Rougier, E.; Lei, Z.; Knight, E.E.; Munjiza, A.; Satriano, C.; et al. Earthquake Damage Patterns Resolve Complex Rupture Processes. Geophys. Res. Lett. 2018, 45, 10279–10287. [Google Scholar] [CrossRef] [Green Version]
  15. Gao, K.; Euser, B.J.; Rougier, E.; Guyer, R.A.; Lei, Z.; Knight, E.E.; Carmeliet, J.; Johnson, P.A. Modeling of Stick-Slip Behavior in Sheared Granular Fault Gouge Using the Combined Finite-Discrete Element Method. J. Geophys. Res. Solid Earth 2018, 123, 5774–5792. [Google Scholar] [CrossRef] [Green Version]
  16. Gao, K.; Guyer, R.A.; Rougier, E.; Johnson, P.A. Plate motion in sheared granular fault system. Earth Planet. Sci. Lett. 2020, 548, 116481. [Google Scholar] [CrossRef]
  17. Godinez, H.; Rougier, E.; Osthus, D.; Lei, Z.; Knight, E.; Srinivasan, G. Fourier amplitude sensitivity test applied to dynamic combined finite-discrete element methods–based simulations. Int. J. Numer. Anal. Methods Geomech. 2019, 43, 30–44. [Google Scholar] [CrossRef] [Green Version]
  18. Osthus, D.; Godinez, H.; Rougier, E.; Srinivasan, G. Calibrating the stress-time curve of a combined finite-discrete element method to a Split Hopkinson Pressure Bar experiment. Int. J. Rock. Mech. Min. 2018, 106, 278–288. [Google Scholar] [CrossRef]
  19. Kalnay, E. Atmospheric Modeling, Data Assimilation, and Predictability; Cambridge University Press: Cambridge, UK, 2003. [Google Scholar]
  20. Daley, R. Atmospheric Data Analysis; Cambridge University Press: Cambridge, UK, 1991. [Google Scholar]
  21. Asch, M.; Bocquet, M.; Nodet, M. Data Assimilation: Methods, Algorithms, and Applications; SIAM: Philadelphia, PA, USA, 2016. [Google Scholar]
  22. Smoljanović, H.; Živaljić, N.; Željana, N. A combined finite-discrete element analysis of dry stone masonry structures. Eng. Struct. 2013, 52, 89–100. [Google Scholar] [CrossRef]
  23. Pulatsu, B.; Erdogmus, E.; Lourenço, P.B.; Lemos, J.V.; Tuncay, K. Simulation of the in-plane structural behavior of unreinforced masonry walls and buildings using DEM. Structures 2020, 27, 2274–2287. [Google Scholar] [CrossRef]
  24. Knight, E.E.; Rougier, E.; Lei, Z. Hybrid Optimization Software Suite (HOSS)—Educational Version; Technical Report, LA-UR-15-27013; Los Alamos National Laboratory: Los Alamos, NM, USA, 2015.
  25. Rougier, E.; Knight, E.E.; Munjiza, A. LANL-CSM: HOSS—MUNROU Technology Overview; Presentation, LA-UR-13-23422; Los Alamos National Laboratory: Los Alamos, NM, USA, 2013.
  26. Knight, E.E.; Rougier, E.; Munjiza, A. LANL-CSM: Consortium Proposal for the Advancement of HOSS; Presentation, LA-UR-13-23409; Los Alamos National Laboratory: Los Alamos, NM, USA, 2013.
  27. Knight, E.E.; Rougier, E.; Lei, Z.; Euser, B.; Chau, V.; Boyce, S.H.; Gao, K.; Okubo, K.; Froment, M. HOSS: An implementation of the combined finite-discrete element method. Comput. Part. Mech. 2020, 7, 765–787. [Google Scholar] [CrossRef]
  28. Munjiza, A.; Knight, E.; Rougier, E. Computational Mechanics of Discontinua; John Wiley and Sons Ltd.: Chishester, UK, 2011. [Google Scholar]
  29. Munjiza, A.; Rougier, E.; Knight, E. Large Strain Finite Element Method: A Practical Course; John Wiley and Sons: New York, NY, USA, 2015. [Google Scholar]
  30. Rougier, E.; Knight, E.; Broome, S.; Sussman, A.; Munjiza, A. Validation of a three-dimensional Finite-Discrete Element Method using experimental results of the Split Hopkinson Pressure Bar test. Int. J. Numer. Meth. Eng. 2014, 70, 101–108. [Google Scholar] [CrossRef]
  31. Cukier, R.; Fortuin, C.; Shuler, K.; Petschek, A.; Schaibly, J. Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. I. Theory. J. Chem. Phys. 1973, 59, 3873–3878. [Google Scholar] [CrossRef]
  32. Cukier, R.; Schaibly, J.; Shuler, K.E. Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. III. Analysis of the approximations. J. Chem. Phys. 1975, 63, 1140–1149. [Google Scholar] [CrossRef]
  33. Brynson, J. Impact Response of Polyurethane. Ph.D. Thesis, School of Mechanical and Materials Engineering, Washington State University, Pullman, WA, USA, 2009. [Google Scholar]
  34. Chen, W.; Song, B. Split Hopkinson (Kolsky) Bar: Design, Testing and Applications; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2010. [Google Scholar]
  35. Broome, S.T.; Lee, M.Y. Dynamic Brazilian Tension Results on Core from Borehole U-15n NNSS in Support of SPE; Technical Report SAND2018-13771R; Sandia National Lab.: Albuquerque, NM, USA, 2018.
  36. Le Dimet, F.; Talagrand, O. Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects. Tellus 1986, 38A, 97–110. [Google Scholar] [CrossRef]
  37. Kalman, R. A new approach to linear filtering and prediction problems. Trans. ASME Ser. D J. Basic Eng. 1960, 82, 35–45. [Google Scholar] [CrossRef] [Green Version]
  38. Metropolis, N.; Ulam, S. The Monte Carlo Method. J. Am. Stat. Assoc. 1949, 44, 335–341. [Google Scholar] [CrossRef]
  39. Jazwinski, A. Stochastic Processes and Filtering Theory; Academic Press: Cambridge, MA, USA, 1970. [Google Scholar]
  40. Doucet, A.; de Freitas, N.; Gordon, N. Sequential Monte Carlo Methods in Practice; Springer: Berlin, Germany, 2001. [Google Scholar]
  41. Evensen, G. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res. 1994, 99, 10143–10162. [Google Scholar] [CrossRef]
  42. Evensen, G. The ensemble Kalman filter for combined state and parameter estimation. IEEE Control Syst. Mag. 2009, 29, 83–104. [Google Scholar] [CrossRef]
  43. Evensen, G. The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean Dyn. 2003, 53, 343–367. [Google Scholar] [CrossRef]
  44. Houtekamer, P.; Mitchell, H. Data assimilation using an ensemble Kalman filter technique. Mon. Weather Rev. 1998, 126, 796–811. [Google Scholar] [CrossRef]
  45. Hamill, T.; Whitaker, J.; Snyder, C. Distance-Dependent Filtering of Background Error Covariance Estimates in an Ensemble Kalman Filter. Mon. Weather Rev. 2001, 129, 2776–2790. [Google Scholar] [CrossRef] [Green Version]
  46. Evensen, G.; van Leeuwen, P. An Ensemble Kalman Smoother for Nonlinear Dynamics. Mon. Weather Rev. 2000, 128, 1852–1867. [Google Scholar] [CrossRef] [Green Version]
  47. Gu, Y.; Oliver, D.S. An iterative ensemble Kalman filter for multiphase fluid flow data assimilation. SPE J. 2007, 12, 438–446. [Google Scholar] [CrossRef]
  48. Emerick, A.A.; Reynolds, A.C. History matching time-lapse seismic data using the ensemble Kalman filter with multiple data assimilations. Comput. Geosci. 2012, 16, 639–659. [Google Scholar] [CrossRef]
  49. Emerick, A.A.; Reynolds, A.C. Ensemble smoother with multiple data assimilation. Comput. Geosci. 2013, 55, 3–15. [Google Scholar] [CrossRef]
  50. Iman, R.L.; Conover, W. Small sample sensitivity analysis techniques for computer models. with an application to risk assessment. Commun. Stat. Theory Methods 1980, 9, 1749–1842. [Google Scholar] [CrossRef]
  51. Rougier, E.; Munjiza, A. MRCK 3D Contact Detection Algorithm. In Proceedings of the Discrete Element Methods: 5th International Conference on Discrete Element Methods, London, UK, 25–26 August 2010; Munjiza, A., Ed.; [Google Scholar]
  52. Munjiza, A.; Andrews, K.; White, J. Combined single and smeared crack model in combined finite-discrete element analysis. Int. J. Numer. Meth. Eng. 1999, 44, 41–57. [Google Scholar] [CrossRef]
Figure 1. Generic setup of a SHPB experiment.
Figure 1. Generic setup of a SHPB experiment.
Applsci 11 02898 g001
Figure 2. Time velocity curve of the SHPB experiment used in the HOSS simulation.
Figure 2. Time velocity curve of the SHPB experiment used in the HOSS simulation.
Applsci 11 02898 g002
Figure 3. Data assimilation cycle for HOSS with EnKF for estimating key model parameters. The default parameters are passed to HOSS. If the error between the model output against observations is large, then the assimilation is performed, which provides new parameter values. The assimilation cycle is performed until the error between the model output and observations is small enough (a certain tolerance).
Figure 3. Data assimilation cycle for HOSS with EnKF for estimating key model parameters. The default parameters are passed to HOSS. If the error between the model output against observations is large, then the assimilation is performed, which provides new parameter values. The assimilation cycle is performed until the error between the model output and observations is small enough (a certain tolerance).
Applsci 11 02898 g003
Figure 4. Ensemble average parameter value (blue line) for the initial sample, generated from LHS using Table 2, and three EnKF assimilation iterations, as well as the reference parameter value (black line) also taken from Table 2. As can be seen, with each iteration, the EnKF provides an estimated parameter value (blue line) that is closer to the reference parameter value (black line). From these experiments, three iterations of the EnKF seem to be sufficient to estimate the correct value for both parameters σ t m a x and E t .
Figure 4. Ensemble average parameter value (blue line) for the initial sample, generated from LHS using Table 2, and three EnKF assimilation iterations, as well as the reference parameter value (black line) also taken from Table 2. As can be seen, with each iteration, the EnKF provides an estimated parameter value (blue line) that is closer to the reference parameter value (black line). From these experiments, three iterations of the EnKF seem to be sufficient to estimate the correct value for both parameters σ t m a x and E t .
Applsci 11 02898 g004
Figure 5. Time–stress curves of the reference simulation (black curve) and for each EnKF iteration, where each plot shows the evolution of the EnKF ensemble, starting with the initial ensemble of 25 members without assimilation (top left, blue curves), ensemble of the first EnKF iteration (top right, green curves), second iteration (bottom left, yellow curves), and third iteration (bottom right, magenta curves). At each iteration it is clearly seen how the EnKF iteration provides an ensemble that is getting closer or tighter around the reference solution.
Figure 5. Time–stress curves of the reference simulation (black curve) and for each EnKF iteration, where each plot shows the evolution of the EnKF ensemble, starting with the initial ensemble of 25 members without assimilation (top left, blue curves), ensemble of the first EnKF iteration (top right, green curves), second iteration (bottom left, yellow curves), and third iteration (bottom right, magenta curves). At each iteration it is clearly seen how the EnKF iteration provides an ensemble that is getting closer or tighter around the reference solution.
Applsci 11 02898 g005
Figure 6. Estimated parameter values for σ t m a x and E t provided by the EnKF data assimilation using SHPB time–stress curve observations, where the blue line is the average parameter value at each iteration, and the green dashed line is the final value at Iteration 3. The EnKF is applied three times, each time providing an estimated parameter value. The assimilation provides a consistent estimate for σ t m a x , where the parameter value is seen to converge around 0.15. In the case of E t , the assimilation seems to be providing monotonically decreasing values of the parameter value, which might not be reliable.
Figure 6. Estimated parameter values for σ t m a x and E t provided by the EnKF data assimilation using SHPB time–stress curve observations, where the blue line is the average parameter value at each iteration, and the green dashed line is the final value at Iteration 3. The EnKF is applied three times, each time providing an estimated parameter value. The assimilation provides a consistent estimate for σ t m a x , where the parameter value is seen to converge around 0.15. In the case of E t , the assimilation seems to be providing monotonically decreasing values of the parameter value, which might not be reliable.
Applsci 11 02898 g006
Figure 7. Time–stress curves for each assimilation iteration, starting with the initial ensemble of 25 members (top-left plot, blue lines), first assimilation iteration (top-right plots, green lines), second iteration (bottom-left plot, yellow lines), and third iteration (bottom-right plot, magenta lines), as well as the observations (black line in all plots).
Figure 7. Time–stress curves for each assimilation iteration, starting with the initial ensemble of 25 members (top-left plot, blue lines), first assimilation iteration (top-right plots, green lines), second iteration (bottom-left plot, yellow lines), and third iteration (bottom-right plot, magenta lines), as well as the observations (black line in all plots).
Applsci 11 02898 g007
Table 1. Elastic properties and density for granite.
Table 1. Elastic properties and density for granite.
ParameterDescriptionUnitsValue
EYoung’s ModulusGPa35.0
ν Poisson’s ratio-0.16
ρ Densitykg/m 3 2550.0
Table 2. Table with the parameter use for data assimilation, with  their corresponding interval values, as well as a “reference value” for the validation of the assimilation algorithm. The  parameters are sampled using the LHS strategy, where each parameter is assumed to have a uniform distribution. The intervals for each parameter are derived in an ad hoc manner from previous HOSS simulations for granite materials.
Table 2. Table with the parameter use for data assimilation, with  their corresponding interval values, as well as a “reference value” for the validation of the assimilation algorithm. The  parameters are sampled using the LHS strategy, where each parameter is assumed to have a uniform distribution. The intervals for each parameter are derived in an ad hoc manner from previous HOSS simulations for granite materials.
ParameterDescriptionUnitsIntervalReference val. (Validation)
σ t max shear strength 10 9 Pa 4.0 × 10 3 , 6.0 × 10 2 4.5 × 10 2
E t specific energy in tangential direction 10 6 J/m 2 2.0 × 10 4 , 24.0 × 10 4 5.0 × 10 4
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Godinez, H.C.; Rougier, E. Assimilation of Dynamic Combined Finite Discrete Element Methods Using the Ensemble Kalman Filter. Appl. Sci. 2021, 11, 2898. https://doi.org/10.3390/app11072898

AMA Style

Godinez HC, Rougier E. Assimilation of Dynamic Combined Finite Discrete Element Methods Using the Ensemble Kalman Filter. Applied Sciences. 2021; 11(7):2898. https://doi.org/10.3390/app11072898

Chicago/Turabian Style

Godinez, Humberto C., and Esteban Rougier. 2021. "Assimilation of Dynamic Combined Finite Discrete Element Methods Using the Ensemble Kalman Filter" Applied Sciences 11, no. 7: 2898. https://doi.org/10.3390/app11072898

APA Style

Godinez, H. C., & Rougier, E. (2021). Assimilation of Dynamic Combined Finite Discrete Element Methods Using the Ensemble Kalman Filter. Applied Sciences, 11(7), 2898. https://doi.org/10.3390/app11072898

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop