Phase-Field Discrete Element Modelisation with solute generation applied to multiple grains.
This work is presented in the following paper : Alexandre Sac-Morane, Manolis Veveakis, Hadrien Rattez. (2024). A Phase-Field Discrete Element Method to study chemo-mechanical coupling in granular materials. Computer Methods in Applied Mechanics and Engineering. https://doi.org/10.1016/j.cma.2024.116900
The goal is to simulate the pressure solution phenomena (described by the following scheme) between multiple grains.
The PFDEM is a discussion between phase-field (PF) simulation and discrete element modelization (DEM). Even if those formulations are different, they can collaborate with data exchange.
The phase-field theory is based on the mesh of the sample and grains are represented by a phase variable (= 1 if the node is a grain, = 0 else). This formulation transmit to the DEM the shape of the grain. Hence, a grain detection algorithm is used to extract the shape from the phase field. In this case, no information about the solute is needed as it is assumed the solute is instantaneously evacuated.
The discrete element model does not need mesh and grains are represented by a list of vertices (here polygonal grains are assumed). The role of the DEM is to compute the position of the grains. Then, an algorithm is used to convert the grain position in a new phase-field. The mechanical energy at the contact level is not considered as the acid is the motor term of the dissolution.
A back and forth is done between phase-field (PF) simulation and discrete element modelisation (DEM). Hence, once a steady-state is reached with the DEM, the phase-field theory is called. The shape of the grains will change and the steady state is not verified and must be computed by a new DEM step. The following schemes illustrate the algorithm.
Phase-field approach is computed with MOOSE (the solver must be correctly installed) and DEM is computed with an own script based on python 3.9.6. The launch of this project is done by executing the python script main.py (multiple run can be asked with multiple_main.py). The user can change variables inside the python script User.py. A section in the following helps to understand the different inputs.
Saves are done before the phase-field simulation and at the end of the DEM. If crash (it seems to happend often during the phase-field simulation with MOOSE) occurs during the simulation, the function main_after_crash.py can be called to come back before the crash.
A documentation for the python script is available.
The goal of the phase-field step is to compute the shape of the grains as it evolves with the precipitation/dissolution. An Allen-Cahn formulation is used coupled with a solute generation and diffusion. Hence, the conservation of the mass should be verified as precipitation/dissolution occurs.
The goal of the DEM step is to compute a steady-state configuration. Grains are assumed as polygonal. The shape is computed from an interpolation of the phase variable. At the end of the step, interactions (grain-grain and grain-wall) are computed to highlight the evolution of the chain force and k0.
Some next releases are presented here.
- add some quick tests to check the script
- update mesh with strain
- increase the user experience (homogenize names, correct english)
- correct the functions P_is_inside() because if the coordinate y = the coordinate of the vertex, the point is considered outside (even if it is not the case)
- create AI or a guide to select algorithm parameters (as neighborhood update frequency, time step, stop criteria, etc...)
- N_grain : the number of the grains
- L_R : is the list of the radius. It defines the particle size distribution of the grains with the variable L_percentage_R
- L_percentage_R : is the percentage of total number of grains with a linked radius. It defines the particle size distribution of the grains with the variable L_R
- x_box_min : the minimum on the x axis of the sample
- x_box_max : the maximum on the x axis of the sample
- n_x : the number of nodes in the x direction
- y_box_min : the minimum on the y axis of the sample
- n_y : the number of nodes in the y direction
- grain_discretization_square : is the discretization of the grains
- Mobility : The mobility of the phase variables
- kappa_eta : the gradient coefficient for the phase variables
- kappa_c : the gradient coefficient for the solute
- Y : the Young modulus
- nu : the Poisson ratio
- rho : the density of the grain. A surface density rho_surf is determined as 4/3*rho*R_mean
- mu_friction_gg : the friction coefficient for grain-grain contact
- mu_friction_gw : the friction coefficient for grain-wall contact
- coeff_restitution : the restitution coefficient for contact
- np_proc : number of processor used for the simulation
- n_t_PFDEM : the total number of PFDEM iteration. A PFDEM iteration is composed by one DEM and one phase-field simulations.
- n_t_PF : approximatively the time step of the phase-field simulation. It defines with dt_PF the total duration of the phase-field simulation
- dt_PF : the time step used for the phase-field simulation. This value is defined with different values. The selection of the value depend on the total absolute energy energy. Criterias to switch level are defined with Ed_level
- factor_etai : is a factor related to the minimal distance between grains with same eta
- dt_DEM : the time step for the DEM simulation
- i_DEM_stop : the maximum iterations done during a DEM simulation
- i_update_neighborhoods: the frequency of neighborhood update
- factor_neighborhood : defines the size of the neighborhood
- Ecin_ratio : defines a stop criteria for the DEM simulation. The steady-state is assumed when the mean displacement during one DEM step is lower than the ratio of the mean radius
- n_window_stop : defines a steps window to detect the steady state
- dy_box_max_stop : in the window defined by n_window_stop, the steady-state is detected if the difference between the maximum and the minimum upper wall position is lower than dy_box_max_stop
- sphericity_margin : is a margin (because of numerical estimation) for grain sphericity compute
- n_spatial_inscribing : is a precision (because of numerical estimation) for grain sphericity compute
- c_min and c_max : are minimum and maximum values for plotting solute concentration
- L_flag_plot : is the list of plots to do. The different options are available. - Config : Grain et solute configuration - DEM_tracker : Tracker during the DEM step - DEM_txt : A .txt file with DEM information at the end of the DEM step - Diff_Solute : Diffusion of the solute (not available for the moment!) - dt : Evolution of the time step used during phase-field iterations. This one depend on the energy introduced in the sample - Ed : Plot the initial energy map - Etai_distribution : Plot the initial distribution of the etai - Eta_c : Plot the evolution of the sum of the etai and the solute - Init_Current_Shape : Plot the shape comparaison (current vs initial) - Kc : Diffusion map - Movie : Movie of the different configurations - Mesh : Plot the mesh on the initial configuration - Porosity : Plot the evolution of the sample porosity (Sg/Sbox) - Sphericity : Evolution of the 2D sphericity of the grains - sum_Ed : Distribution of the energy in the sample - YBoxMax : Evolution of the upper wall position
- struct_element : element used in dilation step to compute node available map
- chi : chemical energy coefficient in the phase-field formulation
- gravity : gravity applied to grains
- Vertical_Confinement_Force : confinement force applied on the upper wall
The coefficient alpha applied to the mechanical energy term is computed from a calibration simulation with the function Add_variables_needed
- n_generation : number of grain generation
- N_test_max : the maximum number of try to generate a grain without overlap with grains already created
- factor_ymax_box : define the upper wall position for the grain generation
- i_update_neighborhoods_gen: frequency of neighborhood update during one generation step.
- i_update_neighborhoods_com: frequency of neighborhood update during multi generation step.
- factor_neighborhood_IC : defines the size of the neighborhood
- i_DEM_stop_IC : the maximum DEM iterations
- dt_DEM_IC : the DEM time step
- Ecin_ratio_IC : defines a stop criteria for the DEM simulation. The steady-state is assumed when the mean displacement during one DEM step is lower than the ratio of the mean radius