.jpg)
1. Introduction
In the geological disposal system consisting of engineered barrier and natural barrier for radioactive wastes, various processes such should be considered to assess the performance of the system. In particular, for the radionuclide transport analysis, geochemical reactions that can accelerate or delay the radionuclide transport should be considered. In order to simulate the radionuclide transport considering coupled thermal-hydraulic-mechanical-chemical processes in a geological disposal system, a modularized process-based total system performance assessment framework, called APro [1], has been developed in Korea Atomic Energy Research Institute. In APro, the multi-physics simulation is basically conducted inside COMSOL [2], but the geochemical reaction simulation relies on an external code, PHREEQC [3]. Within this framework, the reactive transport simulation is done by coupling COMSOL and PHREEQC via Operator Splitting (OS) scheme. Instead of solving governing equations including the mass transport and chemical reaction at once, the two equations are solved by COMSOL and PHREEQC as follows. First, the equation for the mass transport without chemical reaction is solved by COMSOL. Then, the COMSOL solution is transmitted to PHREEQC to solve the equation for the chemical reaction. Finally, the PHREEQC solution is returned back to COMSOL as initial condition for next time step calculation.
The reactive transport simulation codes coupling PHREEQC and COMSOL via OS scheme have been widely developed in the previous studies, but all the codes including APro have connected two programs by exchanging data through other externally operated interface. iCP [4] (interface COMSOL-PHREEQC) of Amphos 21 in Spain has used Java, but most of the codes such as Ecole Polytechnique Federale de Lausanne [5] in Switzerland, SCKCEN [6] in Belgium, University of Ottawa [7] in Canada, and Osaka [8] University in Japan have used MATLAB. APro has also adopted MATLAB as a coupling interface for COMSOL and PHREEQC.
One matter for the coupling techniques used in the previous studies is that the COMSOL simulation should be stopped to deliver the solution to PHREEQC and restarted after the PHREEQC simulation. In order to restart COMSOL simulation, several processes such as stiffness matrix setup for the linear solver are required, so it is inefficient to repeat such processes for restarting COMSOL at every coupling time step.
In this study, a new technique to call PHREEQC and update the solution during the COMSOL simulation is proposed to improve the coupling efficiency. The proposed technique uses the MATLAB function to execute PHREEQC and pointwise constraint to update solutions in the COMSOL interface. Thus, it can be expected to reduce the computation time by omitting unnecessary processes for restarting COMSOL in the conventional coupling techniques.
2. METHODOLOGY
2.1 Operator Splitting Method
Operator Splitting (OS) method has been widely used as a method to more easily solve the reactive transport problem by separating the chemical reaction term, which is burdensome to numerical analysis due to its high nonlinearity, from the transport terms. In the OS method for the reactive transport modeling, the transport equation and chemical reaction equation are separated and solved sequentially. The OS method can be categorized into two approaches. One is Sequential Non-Iterative Approach (SNIA) that the separated equations are solved only once at each time step, and the other is Sequential Iterative Approach (SIA) that the calculations and data exchanges are iteratively conducted until the solutions are converged for each time step [9]. APro adopted SNIA considering computational efficiency as many other codes. In APro, COMSOL solves the non-reactive mass transport problem with multi-physics processes while PHREEQC computes the geochemical reaction.
SNIA in APro can be simply expressed as follows:
where f is the COMSOL operator simulating for one time step, and g is the PHREEQC operator simulating for one time step. is a intermediate state variable that is not updated by g operator at time step k, and uk is an updated state variable at time step k.
2.2 Conventional Coupling Technique Implemented in APro Using MATLAB
The coupling technique implemented in APro is similar to the conventional technique adopted by other codes using MATLAB. The mesh nodes in COMSOL and SOLUTION numbers of PHREEQC are matched based on coordinates to link the two simulators by MATLAB. The non-reactive transport calculation for one time step is performed by COMSOL and the results are extracted and delivered to PHREEQC to compute the geochemical reactions. After finishing the geochemical reaction calculation for one time step in PHREEQC, the results are delivered to COMSOL again for starting next time step calculation. All these data extractions and deliveries are conducted in the MATLAB interface. When using this technique, COMSOL should be stopped and restarted for the PHREEQC calculation, resulting in redundant delay time to setup COMSOL such as stiffness matrix setup for the linear solver, so it is inefficient to repeat such processes just to restart COMSOL at every coupling time step
2.3 New Coupling Technique Calling PHREEQC as MATLAB Function Inside COMSOL
The technique newly proposed in this study is based on an idea that an external MATLAB function can be called during COMSOL calculation using “Livelink for MATLAB” feature supported by COMSOL. As MATLAB functions can be called for any mesh points and every time step in COMSOL, this feature allows users to do various jobs by accessing the variables in COMSOL. In this technique, the MATLAB function in COMSOL is used as a receiver of the solution from COMSOL and a controller to run PHREEQC.
Another main idea for the implementation of the new technique is that the data transmission from PHREEQC to COMSOL is done by “Pointwise constraint” feature in COMSOL. By using this feature, the COMSOL solution is forced to meet the pointwise constraint reflecting the solution from PHREEQC. During the short time when the constraint is applied, COMSOL adjusts the solution to match with the solution from PHREEQC, so the COMSOL simulation can continue with the updated solution. This technique is specially devised because COMSOL cannot have two different results at the same time. Although this technique does not completely reproduce the SNIA concept, it is a kind of approximate technique for updating the results from PHREEQC to COMSOL without stopping COMSOL.
In addition, “State variable” feature in COMSOL, which is supposed to be calculated once for each time step, is used to avoid redundant execution of PHREEQC in this technique. Since the time step used in the transport calculation of COMSOL is shorter than that in coupling COMSOL and PHREEQC, it is necessary to have a technique that allows PHREEQC to be called with a certain waiting time instead of being called for every time step. By using conditional statement, turning on and off the state variables after the waiting time is implemented, and then PHREEQC defined in state variable can be called with specific time interval. That is, the MATLAB function to call PHREEQC is defined in the state variable to represent the PHREQC results, and then the pointwise constraint is defined with the values in the state variable.
A minor technical issue for implementing the new technique is that it needs temporary data storage for saving the results inside the MATLAB function during COMSOL calculation. MATLAB workspace is used for the data storage which can be accessed by COM server. The proposed technique is described in Fig. 1, and compared with the conventional technique.
The main difference of the proposed technique with the conventional one is that COMSOL is not stopped and restarted to run PHREEQC, as in Fig. 2. It is expected that the proposed technique can reduce the redundant time required to restart COMSOL in the conventional technique.
3. Numerical Results
3.1 Verification With the 1D ARD Benchmark Problem
In order to verify that the COMSOL and PHREEQC are properly coupled by the newly proposed approach for the OS scheme, the reactive transport calculation is conducted for a simple benchmark problem to compare with the conventional approach. In this study, a 1D Advection-Reaction-Dispersion (ARD) problem defined as “Example 11-Transport and Cation Exchange” [10] is chosen as the benchmark problem. This problem describes that the chemical composition of the effluent changes when 0.6 mM CaCl2 solution flows into a column containing 1 mM NaNO3 and 0.2 mM KNO3 solution in equilibrium with the 1.1 mM immobile cation exchanger. In this problem, the composition of the effluent changes since K and Na are exchanged with Ca, and the advection and dispersion phenomena also give a great influence on the changes in the chemical composition.
The governing partial differential equation for this problem is defined as
where v is 87.6 m·yr−1, and De is 3×10−10 m2·s−1.
Three ion exchange reactions such as Na+ + X− ↔ NaX, K+ + X− ↔ KX, and Ca2+ + 2X− ↔ CaX2 are considered for the reaction term R in the governing equation.
For the spatial and temporal discretization, the mesh size Δx is set as 0.1 mm for the column length 80 mm, and the time interval Δt is set as 360 seconds for the total simulation time of 72,000 seconds, resulting in 800 cells to compute PHREEQC and 200 time steps for the coupling.
The benchmark problem is described in Fig. 3.
The reactive transport calculation using COMSOL and PHREEQC is conducted by both proposed and conventional techniques. As shown in Fig. 4, the calculation results from the proposed technique agree well with those from conventional technique.
Fig. 4
Concentration of the effluent with time calculated by the proposed and conventional techniques.
3.2 Performance Analysis for the Proposed Technique
Additional calculations for the benchmark problem with varying the mesh size and time step are conducted to measure the computation time. The objective of this analysis is to observe the performance improvement by adopting the proposed technique. There is no significant difference between the calculation results from the proposed technique and conventional technique when varying the mesh size and time step in the verification calculation.
In the first test, the number of cells that need to calculate PHREEQC is doubled or halved by changing the mesh size. This test is to see the effect of the proposed technique according to the increase or decrease of the amount of PHREEQC computation in one time step. The calculation conditions and measured computation times are summarized in Table 1. As shown in the results, at least 30% reduction in computation time is expected in this benchmark problem when using the proposed technique. The performance improvement seems to be more noticeable when the number of cells decreases. This tendency can be explained that the proportion of setting time to restart COMSOL becomes larger as the computation time for PHREEQC decreases.
In the second test, the number of time steps is doubled or halved by changing the time interval. This test is to observe the performance improvement when the number of coupling steps for COMSOL and PHREEQC is changed. The measured results of computation time are summarized in Table 2. As shown in the table, it is confirmed that the degree of performance improvement increases as the number of coupling steps increases, because the proposed technique reduces the delay time in the coupling steps. From the results, it is expected that the proposed technique has a greater effect in the situation where the coupling between PHREEQC and COMSOL is more frequent.
As a third test, the performance analysis is performed on another benchmark problem [4], which is included in the previous study for verification of iCP. The new benchmark problem is also a 1D ARD problem, similar to the previous one, but more complex geochemical reactions are considered in the PHREEQC calculation. The detailed specification is described in the paper [4]. The ratio of computation time from proposed technique and conventional technique are summarized in Table 3. As shown in the table, similar trends are observed in the new benchmark problem when varying the number of meshes and time steps as in the previous performance analysis tests.
4. Conclusion
In this study, a new coupling technique of the OS scheme for COMSOL and PHREEQC to solve reactive transport problems is proposed by calling PHREEQC through the MATLAB function in COMSOL. Unlike the conventional coupling technique using MATLAB outside COMSOL, the proposed technique doesn’t require unnecessary setup time for restarting COMSOL. The accuracy of the proposed technique is verified by 1D ARD benchmark problem, a typical verification problem for the code coupling COMSOL and PHREEQC. In the performance analysis for the proposed technique, it shows that the computation time is reduced by more than 30% in the benchmark problem. It is also confirmed that the degree of performance improvement increases when the coupling process becomes more frequent or the proportion of the coupling process in total calculation time increases. The proposed technique can be applied to APro as well as other codes coupling COMSOL and PHREEQC. Also, it can be used for other coupling objectives, which rely on external codes like PHREEQC. As an example of the application, the proposed method can be useful for developing a machine learning model which is to be integrated into COMSOL as a surrogate model.
