Simulation of Small-Break Loss-of-Coolant Accident Using the RELAP5 Code with an Improved Wall Drag Partition Model for Bubbly Flow

Abstract

The RELAP5 code is a computational tool designed for transient simulations of light water reactor coolant systems under hypothesized accident conditions. The original wall drag partition model in the RELAP5 code has a problem in that the bubble velocity is predicted to be faster than the water velocity in the fully developed flow in a constant-area channel. The wall drag partition model, based on the wetted perimeter concept, proves insufficient for accurately modeling bubbly flows. In this study, the wall drag partition model was modified to account for the physical motion of fluid particles. After that, the modified RELAP5 code was applied to predict the SBLOCA of a full-scale nuclear power plant. Considering the SBLOCA scenario, the behavior change in the loop seal clearing phenomenon was clearly shown in the analysis by the model change. Upon the termination of natural circulation, the loop seals were cleared, allowing the steam trapped within the system to discharge through the break. The modified model was confirmed to have an impact at this time. It mainly affected the timing and shape of the loop seal clearing and delayed the overall progress of the accident. It was observed that the flow rate of the bubbly phase decreased as the modified model accounted for wall friction during dispersed flow in the horizontal section, impacting the two-phase flow behavior at the conclusion of the natural circulation phase.

1. Introduction

Most nuclear safety analysis codes are established based on the two-fluid model [1,2,3,4,5,6,7]. The two-phase flow equation is assumed to be a continuous phase even if it is a dispersed phase, and the same type of equation is applied by averaging each phase.

Among nuclear safety analysis codes, CATHARE [1], COBRA-TF [2], and TRACE [3,4] do not consider the wall drag on the bubble phase, based on the fact that bubbles do not contact the wall surface. RELAP5 [5,6] accounts for wall drag on the bubble phase by employing a wall drag partition model based on the wetted perimeter concept. The RELAP5 code is most commonly used to analyze thermal–hydraulic behavior in light-water nuclear reactors. However, the original wall drag partition model in the RELAP5 code has a problem: In fully developed flow within a constant-area channel, the bubble velocity is predicted to exceed the water velocity. The wall drag partition model, based on the wetted perimeter concept, is inadequate for accurately representing bubbly flows.

In this study, the wall drag partition model was modified for dispersed bubbly flow based on the equation of fluid particle motion. Kim et al. [8,9] theoretically showed that the overall friction pressure drop is partitioned to each phase proportional to its phase fraction. The present study adopted and implemented their theory into RELAP5. The modified wall drag partition model was applied to predict the SBLOCA [10] of a full-scale nuclear power plant.

2. Wall Drag Partition Model for Bubbly Flows

2.1. RELAP5 Model

When the flow is adiabatic, the momentum equations of the gas and liquid phases are given as follows:

$$\begin{array}{ll}{\alpha }_g{\rho }_g{\displaystyle \frac{\partial v_g}{\partial t}}+{\displaystyle \frac{1}{2}}{\alpha }_g{\rho }_g{\displaystyle \frac{\partial v_g^2}{\partial x}}& =-{\alpha }_g{\displaystyle \frac{\partial p}{\partial x}}+{\alpha }_g{\rho }_g{B}_x-{\alpha }_g{\rho }_g\mathrm{FWG}{v}_g\\ & -{\alpha }_g{\rho }_g\mathrm{FIG}(v_g-v_f)\\ & -C_{vm}{\alpha }_g{\rho }_f({\alpha }_g{\rho }_g+{\alpha }_f{\rho }_f)\left[{\displaystyle \frac{\partial (v_g-v_f)}{\partial t}}+v_f{\displaystyle \frac{\partial v_g}{\partial x}}-v_g{\displaystyle \frac{\partial v_f}{\partial x}}\right]\end{array}$$

(1)

$$\begin{array}{ll}{\alpha }_f{\rho }_f{\displaystyle \frac{\partial v_f}{\partial t}}+{\displaystyle \frac{1}{2}}{\alpha }_f{\rho }_f{\displaystyle \frac{\partial v_f^2}{\partial x}}& =-{\alpha }_f{\displaystyle \frac{\partial p}{\partial x}}+{\alpha }_f{\rho }_f{B}_x-{\alpha }_f{\rho }_f\mathrm{FWF}{v}_f\\ & -{\alpha }_f{\rho }_f\mathrm{FIF}(v_f-v_g)\\ & -C_{vm}{\alpha }_f{\rho }_g({\alpha }_g{\rho }_g+{\alpha }_f{\rho }_f)\left[{\displaystyle \frac{\partial (v_f-v_g)}{\partial t}}+v_g{\displaystyle \frac{\partial v_f}{\partial x}}-v_f{\displaystyle \frac{\partial v_g}{\partial x}}\right]\end{array}$$

(2)

This is based on [5,6]. The subscripts $g$ and $f$ denote the gas and liquid phases, respectively. $\alpha$, $\rho$, $v$, $p$, $B_x$, and $C_{vm}$ represent the phasic void fraction, density, velocity, pressure, body force, and virtual mass coefficient, respectively. FWG and FWF are the wall drag coefficients acting on the gas and liquid phases, respectively. FIG and FIF are the interphase drag coefficients for the gas and liquid phases, respectively.

The overall friction pressure drop is computed based on the two-phase friction multiplier as follows:

$${\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }={\displaystyle \frac{{{\lambda }^{\prime }}_f{\rho }_f{({\alpha }_f{v}_f)}^2+C{[{{\lambda }^{\prime }}_f{\rho }_f{({\alpha }_f{v}_f)}^2{{\lambda }^{\prime }}_g{\rho }_g{({\alpha }_g{v}_g)}^2]}^{1/2}+{{\lambda }^{\prime }}_g{\rho }_g{({\alpha }_g{v}_g)}^2}{2D_h}}$$

(3)

where ${{\lambda }^{\prime }}_g$ and ${{\lambda }^{\prime }}_f$ indicate the gas and liquid-alone Darcy friction factors, respectively, calculated at the respective Reynolds numbers $\mathrm{R}{{\mathrm{e}}^{\prime }}_g={\alpha }_g{\rho }_g{v}_g{D}_h/{\mu }_g$ and $\mathrm{R}{{\mathrm{e}}^{\prime }}_f={\alpha }_f{\rho }_f{v}_f{D}_h/{\mu }_f$. $D_h$ is the hydraulic diameter, and $\mu$ is the phasic viscosity. The correlation factor C is computed using the HTFS correlation [11].

Let us consider a fully developed flow in a constant-area pipe. Taking the sum of Equations (1) and (2) gives the following:

$${\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }={\alpha }_g{\rho }_g\mathrm{FWG}{v}_g+{\alpha }_f{\rho }_f\mathrm{FWG}{v}_f$$

(4)

The overall friction pressure drop, which consists of the gas and liquid wall frictions, is partitioned into gas and liquid phases as follows:

$${\alpha }_g{\rho }_g\mathrm{FWG}{v}_g={\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }{\displaystyle \frac{{\alpha }_{gw}{\lambda }_g{\rho }_g{v}_g^2}{{\alpha }_{gw}{\lambda }_g{\rho }_g{v}_g^2+{\alpha }_{fw}{\lambda }_f{\rho }_f{v}_f^2}}$$

(5)

$${\alpha }_f{\rho }_f\mathrm{FWF}{v}_f={\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }{\displaystyle \frac{{\alpha }_{fw}{\lambda }_f{\rho }_f{v}_f^2}{{\alpha }_{gw}{\lambda }_g{\rho }_g{v}_g^2+{\alpha }_{fw}{\lambda }_f{\rho }_f{v}_f^2}}$$

(6)

where ${\lambda }_g$ and ${\lambda }_f$ are the gas and liquid Darcy friction factors, respectively, calculated at the respective Reynolds numbers ${\mathrm{Re}}_g={\rho }_g{v}_g{D}_g/{\mu }_g$ and ${\mathrm{Re}}_f={\rho }_f{v}_f{D}_f/{\mu }_f$.

$D_g=4A_g/P_g$, $D_f=4A_f/P_f$	: phasic hydraulic diameters;
$A_g={\alpha }_{g}A$$, A_f={\alpha }_{f}A$	: phasic flow areas;
$P_g={\alpha }_{gw}P$, $P_f={\alpha }_{fw}P$	: phasic wetted perimeters;
${\alpha }_{gw}$$, {\alpha }_{fw}$	: phasic wetted fraction;
$P$	: wetted perimeter of the channel;
$A$	: cross-sectional area of the channel.

The phasic wetted fractions are determined depending on the flow regimes.

For the bubbly flow regime, the wetted fraction of each phase is assumed to be each phasic fraction: ${\alpha }_{gw}={\alpha }_g$ and ${\alpha }_{fw}={\alpha }_f$. Then, the wall drag coefficients are calculated as follows:

$$\mathrm{FWG}={\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }{\displaystyle \frac{{\lambda }_g\left|v_g\right|}{{\alpha }_g{\lambda }_g{\rho }_g{v}_g^2+{\alpha }_f{\lambda }_f{\rho }_f{v}_f^2}}$$

(7)

$$\mathrm{FWF}={\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }{\displaystyle \frac{{\lambda }_f\left|v_f\right|}{{\alpha }_g{\lambda }_g{\rho }_g{v}_g^2+{\alpha }_f{\lambda }_f{\rho }_f{v}_f^2}}$$

(8)

Applying the wetted perimeter concept to the wall drag partitioning is reasonable for separated flows. However, the wetted perimeter concept is not reasonable for dispersed flows because the dispersed phase may not contact the wall. For an upward bubbly flow, most of the wall surface of the pipe is covered by water, and bubbles do not contact the wall surface. The wall drag acting on the bubble phase should be determined differently.

2.2. Modification to Wall Drag Partition Model

The equation of motion for small particles is as follows:

$$m_p{\displaystyle \frac{d{u}_p}{dt}}=-V_p\nabla p_c+V_p\nabla \cdot {\tau }_c+F_i+m_{p}g$$

(9)

where $m_p$, $u_p$, and $V_p$ are the particle mass, velocity, and volume, respectively. The first and second terms on the right represent the fluid dynamic forces arising from pressure and shear stress associated with the undisturbed flow. The subscripts $p$ and $c$ indicate the particle and continuous phase, respectively. $F_i$ is the fluid dynamic force acting on the particle and includes the steady drag, virtual mass, and history force [12]. $g$ is the gravitational acceleration.

Ref. [13] suggested using the following two-fluid momentum equations when the flow is bubbly:

$${\alpha }_g{\rho }_g{\displaystyle \frac{\partial u_g}{\partial t}}+{\alpha }_g{\rho }_g{u}_g\cdot \nabla u_g=-{\alpha }_g\nabla p+{\alpha }_g\nabla \cdot ({\tau }_l+{\tau }_l^{\mathrm{Re}})+f_i+{\alpha }_g{\rho }_{g}g$$

(10)

$${\alpha }_l{\rho }_l{\displaystyle \frac{\partial u_l}{\partial t}}+{\alpha }_l{\rho }_l{u}_l\cdot \nabla u_l=-{\alpha }_l\nabla p+{\alpha }_l\nabla \cdot ({\tau }_l+{\tau }_l^{\mathrm{Re}})-f_i+{\alpha }_l{\rho }_{l}g$$

(11)

where ${\tau }_l^{\mathrm{Re}}$ is the liquid Reynolds stress tensor, and $f_i$ is the interphase momentum transfer to the dispersed phase. Equations (10) and (11) are based on the equation of motion for small fluid particles. When the flow is one-dimensional, the wall shear stress plays an important role in determining the magnitude of $\nabla \cdot ({\tau }_c+{\tau }_c^{\mathrm{Re}})$. In this case, the term $\nabla \cdot ({\tau }_c+{\tau }_c^{\mathrm{Re}})$ can be interpreted as an overall friction pressure drop ${(-\partial p/\partial x)}_{2\varphi }$. Accordingly, the one-dimensional version of Equations (10) and (11) can be written as follows:

$${\alpha }_g{\rho }_g{\displaystyle \frac{\partial u_g}{\partial t}}+{\alpha }_g{\rho }_g{u}_g\cdot \nabla u_g=-{\alpha }_g\nabla p-{\alpha }_g{\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }+f_i+{\alpha }_g{\rho }_{g}g$$

(12)

$${\alpha }_l{\rho }_l{\displaystyle \frac{\partial u_l}{\partial t}}+{\alpha }_l{\rho }_l{u}_l\cdot \nabla u_l=-{\alpha }_l\nabla p-{\alpha }_l{\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }-f_i+{\alpha }_l{\rho }_{l}g$$

(13)

The terms ${\alpha }_g{(-\partial p/\partial x)}_{2\varphi }$ and ${\alpha }_f{(-\partial p/\partial x)}_{2\varphi }$ correspond to the wall drags acting on the gas and liquid phases, respectively. Consequently, the wall drag terms in Equations (1) and (2) can be modeled as follows:

$${\alpha }_g{\rho }_g\mathrm{FWG}{v}_g={\alpha }_g{\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }$$

(14)

$${\alpha }_f{\rho }_f\mathrm{FWF}{v}_f={\alpha }_f{\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }$$

(15)

The total frictional pressure drop is apportioned among the phases in proportion to their respective volumetric fractions. Subsequently, the wall drag coefficients are determined as follows:

$$\mathrm{FWG}={\displaystyle \frac{1}{{\rho }_g\left|v_g\right|}}{\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }$$

(16)

$$\mathrm{FWF}={\displaystyle \frac{1}{{\rho }_f\left|v_f\right|}}{\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}_{2\varphi }$$

(17)

3. Verification of Wall Drag Partition Model

The test condition was air–water bubbly flow at atmospheric pressure and room temperature. The wall drag partition model was compared with the experimental data and validated in previous studies [14,15,16]. The present study adopted the validated model, implemented it into the RELAP5 code, and applied it to safety analysis.

A horizontal bubbly flow simulation was conducted to check the validity of the wall drag partition model. Figure 1 shows a nodding diagram. The flow channel included the contraction and expansion regions. The total length was 2 m. The large diameter was 0.04 m, and the smaller was 0.02 m. The diameters of nodes 8, 15, and 16 were 0.03 m, 0.025 m and 0.035 m, respectively. The flow was an air–water bubbly flow at 800 kPa and 20 °C. The inlet void fraction was set to 0.01, and the inlet water and air velocities were set to 5 m/s.

Figure 2 shows the simulation result with the original wall drag partition model. The specific gravity of water was significantly higher than that of air. Thus, bubbles accelerated faster than water in the contraction region. Conversely, bubbles decelerated faster than water in the expansion region. These behaviors are observed in Figure 2. Conversely, air was predicted to flow faster than water in the constant-area sections. Figure 3 shows the simulation result with the modified wall drag partition model. The air velocity was predicted to be the same as the water velocity in the constant-area regions, which is physically valid.

4. SBLOCA Analysis of Nuclear Power Plant

The Advanced Power Reactor 1400 MWe (APR1400) is an advanced light water reactor in the Republic of Korea [17,18]. The overall configuration of the APR1400 is illustrated in Figure 4. The APR1400 consists of two loops. Each loop consists of one line of hot leg and two lines of cold leg. The steam generator (SG) and two reactor coolant pumps (RCPs) are installed in each loop, and a pressurizer (PZR) is connected to the hot leg. The safety injection system (SIS) of APR1400 consists of four safety injection pumps (SIPs) and four safety injection tanks (SITs). The safety injection flow from SIS is delivered by direct vessel injections (DVIs). The horizontal pipes include two hot legs, four cold legs, and four pump suction legs.

In this study, RELAP5/MOD3.3/s for calculations [19] of thermal–hydraulic behavior of the reactor coolant system and fuel rod cladding temperature behavior was used to investigate a small-break loss of coolant accident (SBLOCA) in the APR1400. SBLOCA is one of the design-basis accidents (DBAs), which refers to a hypothetical accident caused by the loss of reactor coolant due to the rupture of pipes belonging to the reactor coolant pressure boundary exceeding the supplementary capacity of the reactor coolant supplement system.

Table 1. Steady-state prediction result of the main parameters of the APR1400.

Component	Parameter	Prediction	Difference
Vessel	Core power (102%), MWt	4063	0.00%
	Core flow rate, kg/s	20,332	−0.14%
	Core inlet temperature, K	564	0.00%
	Core outlet temperature, K	599	0.17%
	Total bypass flow, kg/s	629	−0.16%
Steam generator	Pressure, kPa	6914	0.28%
	Steam flow rate, kg/s	1149	−0.43%
	Water volume, m3	117	0.86%
Loop	Vessel flow rate, kg/s	20,960	−0.15%
	Pressurizer pressure, kPa	15,514	0.01%
	Pressurizer water volume, m3	33	0.00%

shows the steady-state prediction result of the main parameters of the APR1400. The plant input data for both models were identical. Under the steady-state condition, a single-phase flow occurred in the primary loop, and a low-speed bubbly flow occurred in the secondary side of the steam generator. Therefore, the wall drag partition model did not affect the steady-state result. The difference means the relative difference between the prediction and design values. The reactor power was assumed to be 102% from the conservative viewpoint. The RELAP5 results agree well with the design values.

Transient simulations were initiated by breaking the direct vessel injection line (Refer to the X mark in Figure 4). The transient calculation results of SBLOCA are shown in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11. The break location and size were selected as direct vessel injection line and 0.01 m², respectively, considering the SBLOCA limiting break. A time step of 0.05 s was used for the 2000 s simulation until the SBLOCA was terminated by the emergency core cooling system (ECCS).

Figure 5 shows the pressure in the primary side of the reactor in the event of an SBLOCA. Overall, the pressure trend of the original wall drag partition model is similar to that of the modified wall drag partition model. However, it started to change slightly from around 300 s after the pressure quasi-equilibrium state ended, and it can be seen that it changed distinctly after 750 s. The pressure of the original wall drag partition model tended to decrease faster than the modified wall drag partition model.

Figure 6 shows the break flow rate. In SBLOCA, the break flow rate was released as a supercooled or saturated phase at the time of the accident. The break flow rate also showed a slight difference from around 300 s. Figure 7 shows the flow rate of coolant injected into the core. The cooling water injected into the core rapidly decreases due to the shutdown of the reactor coolant pump. The core inlet flow rate also slightly changed from around 300 s.

Figure 8 shows the collapsed water level in the core. A slight difference started to occur around 300 s, and after 500 s, the core level in the original wall drag partition model tended to decrease rapidly, unlike the result of the modified wall drag partition model. There was a remarkable water level difference between 700 and 1200 s. When the original wall drag partition model was used, the core water level decreased significantly in the period between 700 and 1000 s, and the core level recovered after 1000 s. On the other hand, when the modified wall drag partition model was used, the core water level decreased slowly and recovered after 1200 s. The water level significantly affected the peak cladding temperature.

Figure 9 and Figure 10 show the steam flow rate indicating loop seal clearing [20]. Before loop seal clearing, steam was trapped in the system, and the core water level decreased rapidly, leading to core exposure, or the pressure continued to drop after the loop seal clearing, exposing the core before additional emergency core cooling water was supplied, resulting in peak cladding temperature. Therefore, the peak cladding temperature depends on the occurrence of loop seal clearing. If loop seal clearing occurs in one loop, the amount of coolant supplied to the core is limited, the core’s uncover time is quickened, and the peak cladding temperature occurs. Loop seal clearing is one of the uncertain phenomena, and it is unclear which of the four pump suction pipes is cleared first near the boundary break, so it is currently interpreted as a conservative approach.

One loop seal was cleared in the case of the original wall drag partition model, as shown in Figure 9. However, in the case of the modified wall drag partition model in Figure 10, it can be confirmed that two loop seals were cleared around 300 s. This is because the thermal–hydraulic behavior in the loop seal was changed as the speed of the bubble phase was limited in the modified wall drag partition model. Due to this difference, it was confirmed that the pressure, break flow rate, core inlet flow, and core level in the reactor exhibited differences from around 300 s. Due to this effect, as shown in Figure 11, the peak cladding temperature of the modified wall drag partition model that occurred in SBLOCA decreased by about 160 °C.

Table 2. Comparison of sequential event times between the original and modified wall drag partition models. SIP: safety injection pump, LSC: loop seal clearing, SIT: safety injection tank, PCT: peak cladding temperature; ⓐ, ⓑ, ⓒ, ⓓ: event location indicated in Figure 4.

Wall Drag Partition Model	SIP ⓐ	LSC ⓑ	Core Uncovery ⓒ	SIT ⓓ	PCT
Original	78	284	732	1014	1016
Modified	78	286	880	1124	1126

compares the time of major events in the original and modified wall drag partition models in SBLOCA. The important point is that in the modified wall drag partition model, all events were delayed, and in particular, a clear difference occurred after the loop seal clearing occurred. This is understood to be due to the slow thermal–hydraulic behavior and the slow pressure drop due to the consideration of friction on the bubble phase within the horizontal tube in the case of the modified wall drag partition model.

The bubble velocity is governed by the interfacial drag, wall drag, and gravity. In a horizontal pipe, such as intermediate legs, the wall drag on the bubble phase significantly affects the bubble velocity. In the original wall drag partition model, the total wall drag was partitioned into each phase based on the wetted perimeter concept. However, in the modified wall drag partition model, the total wall drag was partitioned into each phase in proportion to its fraction in each phase. As a result, compared to the original model, the wall drag on the bubble phase decreased, and the wall drag on the water phase increased in the modified model. The change in the wall partition significantly affected loop seal clearing in the intermediate legs and the amount of coolant in the core.

A typical sequence of events for SBLOCA was blowdown, natural circulation, loop seal clearing, core uncovery, and core recovery, leading to the accident. Of them, loop seal clearing strongly influenced core exposure, especially regarding the peak cladding temperature. In conventional SBLOCA scenarios, a single loop seal clearing typically resulted in the most limiting conditions for peak cladding temperature due to minimal coolant return to the core and rapid pressure drop, leading to a rapid decrease in core coolant level. In contrast, the modified wall drag model exhibited two loop seal clearings, increasing coolant recovery within the core and reducing gas release after loop seal clearing, thereby slowing the pressure decrease in the primary system.

5. Conclusions

In this study, the wall drag partition model in the RELAP5 code was modified in consideration of fluid particle motion. Next, the modified model was verified for a one-dimensional bubbly flow in a horizontal channel with contraction and diffuser sections. The modified wall drag partition model showed physically valid bubble and water velocities. After that, the effects of the modified model were reviewed through the analysis of the SBLOCA calculation of a full-scale nuclear power plant.

The break caused a rapid depressurization of the reactor coolant system (RCS) at the onset of the SBLOCA calculation. When the RCS pressure reached the reactor trip setpoint, the reactor was tripped by the insertion of control rods. During the blowdown phase, the RCS was predominantly filled with the liquid phase, and the subcooled or saturated coolant was discharged into the containment building through the break. As the RCS continued to depressurize and reached a quasi-equilibrium state with the secondary side of the SG, the blowdown phase was concluded. No significant differences were observed during the blowdown phase when the modified wall drag partition model was applied.

In the natural circulation phase after the blowdown phase, the thermal quasi-equilibrium state could last for hundreds of seconds. Since the RCP was shut down due to the assumption of loss-of-offsite power, there was no forced circulation flow during this period. As the loop seal in the form of a horizontal tube was filled with the coolant, it could not form the effective flow path for steam. When natural circulation was terminated, the loop seals were cleared, and the steam isolated in the system was discharged through the break near 300 s. The modified model was confirmed to have an impact at this time. The modified model affected the thermal–hydraulic behavior from the natural circulation phase rather than the forced circulation. It mainly affected the timing and shape of the loop seal clearing and delayed the overall progress of the accident. It was understood that the flow rate of the bubbly phase decreased as the modified model considered wall friction during dispersed flow in the horizontal section, which affected the two-phase flow from the end of the natural circulation phase.

In the future, we will also explore the applicability of this approach in other fields to further enhance its development [21].