Pennes’ bioheat equation vs. porous media approach in computer modeling of radiofrequency tumor ablation

The objective of this study was to compare three different heat transfer models for radiofrequency ablation of in vivo liver tissue using a cooled electrode and three different voltage levels. The comparison was between the simplest but less realistic Pennes’ equation and two porous media-based models, i.e. the Local Thermal Non-Equilibrium (LTNE) equations and Local Thermal Equilibrium (LTE) equation, both modified to take into account two-phase water vaporization (tissue and blood). Different blood volume fractions in liver were considered and the blood velocity was modeled to simulate a vascular network. Governing equations with the appropriate boundary conditions were solved with Comsol Multiphysics finite-element code. The results in terms of coagulation transverse diameters and temperature distributions at the end of the application showed significant differences, especially between Pennes and the modified LTNE and LTE models. The new modified porous media-based models covered the ranges found in the few in vivo experimental studies in the literature and they were closer to the published results with similar in vivo protocol. The outcomes highlight the importance of considering the three models in the future in order to improve thermal ablation protocols and devices and adapt the model to different organs and patient profiles.

Spatial coordinate (m) Z Impedance (Ω Interest in thermal ablation as an anticancer technique has grown through the years since the clinical data supports its use against different types of tumor, especially liver tumors, and it is a minimally invasive technique with certain advantages, such as shorter hospital stays and lower costs 1,2 . Different types of energy can be used to raise tissue temperature over 50 °C, such as radiofrequency (RF), microwaves or ultrasound to completely destroy tumor tissue while saving the surrounding healthy tissue. In this context, the mathematical modeling of heat transfer in biological tissue during thermal ablation has a key role to play in predicting coagulation zone volume (i.e. thermal lesion size)according to different factors, such as tumor size and shape. The results thus obtained may also give rise to the development of new medical devices and protocols for thermal ablation. The importance of computational modeling in this field is due to the few in vivo studies carried out, which are in general expensive and time-consuming.
Various bioheat models have been developed since 1948, when Pennes proposed the simplest bioheat equation 3 which is still widely used for its simplicity but nevertheless presents certain shortcomings. As it assumes a uniform perfusion rate, it does not consider blood flow direction and also neglects the artery-vein countercurrent arrangement. It also assumes that venous blood is the only type in thermal equilibrium with the tissue and that arterial blood remains at a constant temperature of 37 °C. Some researchers have tried to overcome these limitations by proposing different models to provide more accurate results. An in-depth review of the different bioheat models applied to hyperthermia treatments can be found in Andreozzi et al. 4 .
In this study, we focused on modeling liver radiofrequency ablation (RFA) and compared the performance of Pennes' bioheat equation with other bioheat equations based on porous media theory in terms of coagulation zones and temperature distributions. According to the porous media theory, the entire biological medium can be divided into two distinct phases: the tissue phase, which is the solid part made up of cells and interstitial spaces, and the blood phase, i.e. the fluid part, which is represented by the blood flowing through the solid phase 5 . As will be explained below, the governing equations for the bioheat transfer and blood flow are thus averaged on a control volume which represents a fluid-saturated porous medium infiltrated by flowing blood. The blood volume fraction in the whole biological medium is represented by porosity, so the higher the vascularization of the tissue the more important will be its role in the model, such as in the liver, kidney, and tumors themselves.
Our motive for this study was that more realistic models based on porous media could predict tissue temperature more accurately, especially in highly porous organs such as the liver 6 . In addition, since the porous media theory includes fewer assumptions than other bioheat equations, it is reasonable to think that it is more www.nature.com/scientificreports/ suitable for treating heat transfer in biological tissues 5,[7][8][9][10][11] . In fact in the models based on porous media, tissue and blood temperatures are calculated separately, ignoring Pennes' assumption of uniform 37 ºC blood temperature throughout the tissue, which is unrealistic in thermal ablation, where tissue temperatures can go beyond 100 °C near the probe. The porous media theory also considers blood flow in different directions to represent the real vascular structure, so that the other strong assumption of uniform perfusion is not included. The differences between Pennes' equation and the porous media-based models described above are graphically illustrated in Fig. 1. We considered two different models based on porous media theory: the Local Thermal Equilibrium equation (LTE) and the Local Thermal Non-Equilibrium equations (LTNE) 5 . In both models we modified the governing equations to consider vaporization in both tissue and blood and assumed uniform blood velocity in all directions to simulate a more realistic blood flow in a vascular network. All the models built were designed for the particular case of in vivo liver RFA with a cooled electrode and an impedance-controlled pulsing protocol as described in 12 , which represents the usual clinical scenario in hepatic tumor RFA 13 .

Results and discussion
Analysis of outcomes. The results are given in terms of the minor diameter of coagulation zones d c (transverse diameter to the applicator shaft in r direction in Figs. 2 and 3), total energy delivered during the application E RF , maximum tissue temperature reached T t,max , total time in which the generator is "on" (t on ), time of first roll-off (t roll-off ) and number of roll-offs N roll-off , distinguishing between the Pennes' , LTNE and LTE results with four different porosity values.
All the models employed were compared at voltage values of 45, 65 and 90 V for720 s. As expected, only in the 45 V case did the absence of roll-offs allow RF power to be applied continuously for the entire 720 s. When roll-offs occurred at 65 and 90 V the pulsing protocol was activated 12 .
Low voltage RFA. Table 1 shows the results of the 45 V simulations. The first finding was that the coagulation diameters computed from the LTE and LTNE models were smaller than those from the Pennes' model for all porosity values. For instance, the difference between Pennes' and LTE ranged from only 1.30 × 10 -3 m for ε = 0.2 to 12.4 × 10 -3 m for ε = 0.6. The coagulation diameter reduced dramatically as porosity increased (from 0.2 to 0.6) and the value of the maximum temperature in tissue also dropped. In fact, the higher the porosity, the larger the convective contribution of the mass blood flux.
At lower temperatures electrical conductivity increased less, so impedance (Z) decreased less and RF applied power decreased too (P = V 2 /Z). This was confirmed by the reduced energy (~ 19.5 kJ to 16.5 kJ when porosity rose from 0.2 to 0.6). This highlights the importance of considering the blood volume fraction, (i.e. porosity) in bioheat thermal ablation models. The porosity of different organs ranges from negligible values such as the brain (0.03 to 0.05) to very high as in liver (about 0.3) and kidney (about 0.35) 6,14-19 .   www.nature.com/scientificreports/ The same organ can have different porosities in different physiological conditions, as happens for example in chronic liver hepatitis and cirrhosis, when porosity can be less than 0.2 20 . Figure 2 shows temperature distributions and coagulation contours at 45 V and 720 s computed from the Pennes' , LTE and LTNE models, respectively. The mean temperature dropped in the domain as porosity increased as given in Fig. 2.

Tissue cells and interstitial spaces
The results of the LTE and LTNE models at the same porosity value were almost identical. As all the cases referred to a tissue with infiltrating terminal arteries, as described in Section "Modified LTNE model", the values of blood vessel diameter and blood velocity in the LTNE model were small enough to validate the local thermal equilibrium assumption, in agreement with the results reported in 7,21,22 , which all confirmed that the LTE temperature distributions agree with those of LTNE only when small blood vessels are considered (up to 3.00 × 10 -5 m) and blood velocity is less than 4.00 × 10 -3 m s -1 , showing that blood does not act as a heat sink in these cases. Note that the LTNE and LTE equations are not limited to modeling vessels smaller than 1 mm, but can also model larger vessels, which were not taken into account in this work. Table 1. 45 V RFA results computed for different bioheat models (Pennes' , Local Thermal Equilibrium, Local Thermal Non-Equilibrium) and porosity values. d c : coagulation diameter; E RF : applied RF energy; t on : total time generator is "on"; t roll-off : time of first roll-off; N roll-off : number of roll-offs.  Table 2. Results for different bioheat models at 65 V. d c : coagulation diameter; E RF : applied RF energy; t on : total time generator is "on"; t roll-off : time at first roll-off; N roll-off : number of roll-offs.  Table 3. Results for different bioheat models at 90 V. d c : coagulation diameter; E RF : applied RF energy; t on : total time generator is "on"; t roll-off : time at first roll-off; N roll-off : number of roll-offs.  Tables 2 and 3 show the results for 65 and 90 V, respectively. Unlike 45 V, the coagulation diameters computed from the LTE and LTNE models with 65 and 90 V at all porosity values were greater than those from the Pennes' model, except for LTNE ε = 0.6. This can be explained by their different ways of applying power: at 45 V it was continuous for 720 s, so that the only physical phenomenon that affected coagulation zone size was the larger heat loss through blood as porosity increased, which also meant less energy delivered in the LTNE and LTE models. Instead, the higher voltage values involving alternating periods of rising (power on) and falling (power off) temperatures weighed on the different thermal inertia of the models. In fact, unlike Pennes' equation, they considered solid and fluid phases separately at different temperatures and water contents, which resulted in a better heat storage capability in the "off " periods. Consequently, the higher heat storage capability became determinant in achieving necrosis, especially away from the electrode, so larger coagulation diameters were finally obtained. This can also be seen in Fig. 5 concerning the comparison with experimental results. As for 45 V, at these voltages the larger the porosity (from 0.2 to 0.6), the smaller the coagulation diameter (from 4.09 × 10 -2 m to 3.17 × 10 -2 m for LTE, and 3.72 × 10 -2 m to 2.23 × 10 -2 m for LTNE at 65 V, from 4.13 × 10 -2 m to 3.16 × 10 -2 m for LTE, and 3.94 × 10 -2 m to 2.62 × 10 -2 m for LTNE at 90 V).

Pennes
Interestingly, the smaller diameter was associated with a slight increase in delivered energy: from 28.6 to 32.1 kJ for LTE and from 31.2 to 34.2 kJ for LTNE at 65 V, and 33.8 to 36.2 kJ for LTE and 38.1 to 39.8 kJ for LTNE at 90 V. This was probably due to roll-offs again, involving different thermal inertia in the LTE and LTNE models.
When porosity is increased, the larger blood volume removes heat from the tissue more effectively, which simultaneously delays roll-off, i.e. power can be applied longer (higher values of t on ), requiring slightly higher power in LTE and LTNE.
The maximum tissue temperature was quite similar in all cases at the same voltage level (~ 114 ºC for 65 V and ~ 120 ºC for 90 V) due to the fact that in all cases the on-off periods avoided the temperature rising above the limit value.
When the LTE and LTNE models were compared at the same porosity value, the LTNE coagulation diameters were similar to or even smaller, with longer t on and higher energy, possibly because of the effect of vaporization at different temperatures for tissue and blood, which is more significant for these applications than at 45 V.
For the same bioheat model, the results obtained at the two different voltage levels were almost identical, (differences in coagulation diameter from 1.00 × 10 -4 m to 7.00 × 10 -4 m in LTE and from 3.00 × 10 -4 m to 3.90 × 10 -3 m for LTNE), which suggests that 65 V is enough to obtain maximum coagulation diameter after 12 min. Figure 3 shows temperature distributions and coagulation contours at 65 V and 720 s computed from the Pennes' , LTE and LTNE models, respectively. Comparing Figs. 2 and 3 shows the different coagulation shapes at voltages higher than 45 V. In fact, at 65 V the zones are more spherical for both the LTE and LTNE models than Pennes' , especially at low porosities. Figure 4 summarizes the coagulation zone diameters computed from the three bioheat models at different porosities (models based on porous media approach).
As highlighted in Figs. 2, 3, and 4, the differences in terms of coagulation diameters and temperature distributions could differ significantly or not at all between the porous media-based and Pennes' models, according to applied voltage and porosity. For instance, at 45 V, Pennes' provides the same result as LTE and LTNE for ε = 0.2. In fact, the range of differences in coagulation diameters obtained is from only 1 mm at 45 V and LTE ε = 0.2 to about 1.60 × 10 -2 m for the highest voltage applied and LTE ε = 0.2. These differences could play a relevant role in predicting coagulation zones, since the risk could be either incomplete ablation and tumor recurrence or overestimating the ablated area and healthy tissue necrosis.  www.nature.com/scientificreports/ 3.42 × 10 -2 m and 3.47 × 10 -2 m respectively, as shown in Table 3. Note that the exact protocol is unknown for the specific case, but the predominant current peak ofabout 1600 ± 71 mA is comparable with the 1500 mA predominant peak obtained in our case. Although these differences could slightly affect the results, they did not affect the overall comparison of all the models. Figure 5 compared the tissue temperature evolution at 10 mm and 20 mm from the electrode between the experimental results in 23  It can be seen that the temperature difference between the porous media-based models and the experimental data are ~ 4 °C and 2 °C at 10 mm and 20 mm, respectively. However the slightly higher predominant current peak could partially justify this difference. At 10 mm the temperature differences between the two porous mediabased models and Pennes' equation are only about 2 °C, mostly in the last 200 s, so the three bioheat models come reasonably close to the experimental data. However, at 20 mm, the differences between both the LTE and LTNE models and Pennes' equation are ~ 8 °C, which justifies the difference of about 0.01 m in the final coagulation diameters obtained and suggests that the porous media-based models could be more accurate than Pennes' in this case. The LTE and LTNE model results confirmed that the thermal local equilibrium is maintained in this case, as explained previously, so that they also give very similar outcomes in terms of tissue temperature.

Limitations of the study. The main limitation of this study is the lack of accurate experimental data on
in vivo studies with which to compare our results, which means it is impossible to assess which model offers the best prediction of coagulation zone size in all conditions. By taking the currently available data into account, we found that LTNE and LTE (which must be considered more realistic bioheat models than Pennes) provided results in general within the range of values reported for experimental studies and approached the experimental data with similar protocols. As new values are reported in the future, it will be possible to determine which model best predicts thermal behavior during RFA and under what conditions. Uniform porosity was assumed in the LTNE and LTE models, even though some studies showed porosity changes from the core to the rim of the tumor and in the adjacent normal tissue 26 . However, the same changing porosity could be implemented in both LTE and LTNE equations, so the comparison between the models would give the same conclusions.
The electrical conductivity used for tissue was for healthy liver tissue. Although it is known that inside a tumor its value can be slightly higher 30 (0.45 S/m instead of 0.19 S/m at 37 °C), we consider that the effect of changing σ on the results would affect all three models equally, so that the conclusions would remain unchanged.
Finally, as the comparison was performed only for a single RF probe it does not seem likely that considering other RF applicator designs would change the conclusions in qualitative terms.
Overall, all these shortcomings should be solved in the future by means of developing experimental setups under strictly controlled conditions and with very accurate temperature measurements. From a computational point of view, the improvements could be focused on developing new models with variable porosity inside the tumor domain, i.e. from the center to the rim, and also considering differences in electrical conductivity between tumor and healthy tissues. All this would contribute to the development of realistic models tailored to specific organs and diseases.

Conclusions
The aim of this work was to compare three different heat transfer models for radiofrequency ablation (RFA) of in vivo liver tissue using a cooled electrode. It is important to point out that our goal was not to suggest the best model for predicting the coagulation zone in RF tumor ablation, but to point out the differences obtained by comparing them. The study was conducted implementing a low voltage and two high voltage levels in order to consider cases with and without roll-off. As far as we know, this has never been done for in vivo liver thermal ablation. When the computer results obtained were compared with the few in vivo experimental works available in the literature it was found that they generally covered all the value ranges. However, for a similar protocol (e.g. 90 V pulsed protocol with a predominant current peak of about 1500 mA and 15 s off periods), the porous media model achieved a better agreement with the experimental results, since Pennes' bioheat model tended to underestimate temperature distributions for the cases here investigated. Furthermore, from the theoretical point of view the porous media models represented a valid alternative approach that tried to overcome the constraints of Pennes' bioheat equation and could be useful when additional tissue characteristics (such as porosity) become available. Despite the paucity of in vivo studies and the need for further experimental studies, the three computer models should be taken into account in the future in order to go on improving both medical protocols and devices. Figure 6 shows the model geometry, which includes an RF applicator (comprised of metal electrode and plastic handle) completely inserted in the hepatic tissue. The needle-like electrode has a conical tip, 3.00 × 10 -2 m long and 1.50 × 10 -3 m outer diameter, and an internal cooling circuit with cold saline. The axial symmetry meant that a two-dimensional model could be implemented.

Model geometry and materials properties.
RFA is a monopolar procedure in which electrical current flows between the ablation electrode (see Fig. 6) and a large size dispersive electrode (in contact with the patient's skin). The RF generator frequency is 500 kHz, and the electric and thermal properties of the materials used in the models are shown in Table 4. Since the temperature dependence of thermal conductivity is very small this was considered constant.
Electrical conductivity in tissue and blood was defined as a temperature-dependent piecewise function to consider the drastically reduced water content loss as temperature increases and vaporization occurs 34 . Electrical problem equations. The heat source from RF power Q ext (Joule losses) was given by: where E is the electric field. E = − ∇V was obtained from the governing equation of the electrical problem ∇·(σ(T) ∇V = 0), σ being the electrical conductivity and V the voltage. www.nature.com/scientificreports/ A constant voltage was set on the ablation electrode during the entire protocol duration of 720 s except during roll-offs, when impedance rises by 30 Ω, in which case we followed the standard roll-off protocol of switching off the applied voltage for 15 s 12 .
The three voltage values employed were 45 V, 65 V and 90 V. This choice covers the range from the highest roll-off avoidance value (45 V) to the highest standard protocol value used in clinical practice (90 V) 12 .
For the electrical boundary conditions, an insulation electrical condition was applied to the contours of the tissue, except for top and bottom, which represent the dispersive electrode, where a condition V = 0 V was set. Thermal problem equations. Three different bioheat models were implemented: the Pennes' , LTNE and LTE. We first described the issues common to all three models. Initial and boundary tissue temperature were 37 °C. The thermal cooling effect inside the electrode is observable from the colder zone around the electrode in Fig. 7, where a detail of Fig. 3 is shown. This effect was modeled by applying a convective heat flux (q c ) to the internal boundary of the electrode as follows: where T t is the tissue temperature and T r the coolant temperature (5 °C), h r is the thermal convective coefficient (3127 W m -2 K -1 ), which was calculated by considering the electrode length (3.00 × 10 -2 m) and a flow rate of 45 ml/min through an area equivalent to half the cross section of the inner diameter of the electrode, as described in 12 . The thermal damage was computed using the Arrhenius model, in which tissue damage is increased with time of exposure, and is obtained as follows from an exponential relationship between tissue exposure temperature, time and the kinetic parameters that account for morphologic changes in tissue relating to the thermal degradation of proteins 35 :  www.nature.com/scientificreports/ where A is the frequency factor (7.39 × 10 39 s −1 ), ΔE is the activation energy for the irreversible damage reaction (2.577 × 10 5 J mol −1 ) 36 , and R is the universal gas constant. The thermal damage contour was estimated using the isoline Ω = 4.6, which encompasses the zone with 99% cell death probability.
Pennes' model. Pennes' bioheat model was our starting point for the comparison of the models, based on the following equation: where subscript t refers to the tissue, ρ is the density, c the specific heat, k t the thermal conductivity, T t the temperature, t the time, Q perf the blood perfusion term, Q met the metabolic heat source (which is neglected in RFA applications 12 ) and Q ext is the external RF power density imposed to the tissue during RF power application (see Eq. (2)). To introduce vaporization into the Pennes' equation, one alternative is to follow the enthalpy method as described in 31 , so that the term (ρ t c t ) in Eq. (5) becomes: where ρ l and c l are density and specific heat of tissue at temperature below 100 °C (liquid phase), ρ g and c g are density and specific heat of tissue at temperature above 100 °C (gas phase), h fg is the product of water latent heat of vaporization and water density at 100 °C (2.17·10 9 J m -3 ), and C w,t is the water content inside the liver tissue (68%) 32 . Furthermore, Q perf in Eq. (5) is defined as follows: where ρ b and c b are the density and specific heat of blood respectively, ω b is blood perfusion coefficient (0.019 s -1 ) 37 , T b is blood temperature (which is assumed to be constant throughout the tissue and taken as body temperature of 37 °C in Pennes' model), and β is a coefficient that is 0 or 1 depending on the value of thermal damage function Ω (see Eq. (4)), so, β = 0 for Ω ≥ 4.6 and β = 1 for Ω < 4.6.
Modified LTNE model. The second model was based on a modified form of the Local Thermal Non-Equilibrium equations, first developed in 1997 by Roetzel and Xuan 38 to model heat transfer in a porous medium. Here, the entire biological medium is divided into two distinct phases: the tissue phase, which is the solid part, made up of cells and interstitial spaces, and the blood phase, which is the fluid part, represented by the blood flowing through the solid phase. We thus have two energy equations for this model, one for the tissue temperature (T t ): and one for the blood temperature (T b ): www.nature.com/scientificreports/ where ɛ is porosity (representing the blood volume/total volume ratio), u the blood velocity vector, β is the coefficient that is 0 or 1 depending on the value of thermal damage function Ω employed previously, h c the interfacial heat transfer coefficient between tissue and blood phases, and a the volumetric heat transfer area between tissue and blood, which is computed from the definition of hydraulic diameter as a = 4·ɛ/d, where d is the blood vessel diameter. The second term on the right side of Eqs. (8) and (9) describes the interfacial convective heat transfer between blood and tissue phases across the vascular wall as defined by Newton's cooling law. Note that the Pennes' model does not consider any advective term such as the second term on the left side of Eq. (9), but an overall blood perfusion term as a heat sink for tissue. The LTE and LTNE models split the Pennes' equation perfusion term into a modified perfusion term and a convective term 39 . The modified perfusion term in the porous media-based models (third term on the right side of Eqs. (8) and (9)) differs from Eq. (7) for the value of the perfusion coefficient ω, which is ω = 0.0005 s -1 , instead of ω b = 0.019 s -1 , to consider the modification described above. In fact, Pennes obtained this coefficient for the volume flow of blood through tissue by fitting his model with experiments. He specified that 0.0005 s -1 is the most suitable value when the balance between blood and tissue is incomplete. Nakayama et al. also used this value in their work on LTNE equations 3,39 .
The volume averaging technique is employed to consider the volume average quantities of the variables 40 , so that the symbol < > refers to the average volume of a generic variable and is neglected from this point onwards. In the LTNE model, the most important modifications were to the phase change and blood velocity. In fact, the phase change was considered separately in both phases, while vaporization only refers to water in the tissue in the Pennes' model (see Eq. 7). Vaporization was thus included as in Eq. (6) for the tissue phase, while for the blood phase we have: In this way the evaporation of water in blood was included, choosing the value of water content in blood C w,b = 79% as in 41 . Four different blood velocity directions were considered to simulate a more realistic vascular network and the initial values were the same in all four directions u z = u -z = u r = u -r , where directions z and r are specified in Fig. 6. As for blood perfusion, blood velocity was considered zero when cell death probability was 99% according to the β coefficient. Terminal artery blood velocity was chosen following the experimental values reported in Crezee and Lagendijk 42 , a reasonable choice according to Chen and Holmes' LTNE model 43 , in which blood heat exchange is assumed to occur only downstream of the terminal arteries before the arterioles. Blood velocity value was thus assumed to be 4.00 × 10 -3 m·s -1 , with a blood vessel diameter d = 3.00 × 10 -5 m. We also considered four different porosity values, ɛ: 0.2, 0.3, 0.4 and 0.6 to consider the liver values given in clinical and numerical studies in the literature 6,14,20,24,25 .These, in fact show wide dispersity because of the different assessment methods considered 20,25 .The interfacial heat transfer coefficient h c was 170 W·m -2 ·K -1 , as in Yuan 21 , based on experimental measurements. Table 5 summarizes the specific surface areas a and volumetric convective coefficients h v for all the cases considered.
Modified LTE equation. The LTNE model described above considered the blood phase and the tissue phase at two distinct temperatures (i.e. T t ≠ T b ). However, when the local thermal equilibrium hypothesis is maintained, tissue and blood are really at the same temperature, so that T t = T b = T and Eqs. (8) and (9) can be combined in a single equation as follows: In this model the heat sink effect for tissue is only represented by the advective term related to the blood velocity. In fact, when Eqs. (8) and (9) were combined the perfusion term disappeared. Note that in this case local thermal equilibrium should be a good approximation for temperature distributions because of the small size of the vessels considered (as described in 2,7,22 ), while this assumption would not be valid in the presence of larger vessels. Even if the two phases are at the same temperature, they have different water contents, as in the LTNE model, so that vaporization was included as in Eqs. (6) and (10) by considering T t = T b = T. The assumptions on blood velocity for LTNE were also valid for this model. www.nature.com/scientificreports/ Numerical solution. The three models were numerically solved with the software Comsol Multiphysics (Burlington, MA, USA). A triangular mesh was employed with a finer size on the boundary between electrode and tissue domains as in 12 , where the highest temperature gradients take place. The grid convergence was verified on the maximum tissue temperature (T t,max ) obtained at first roll-off time. When the difference between simulations was less than 0.5% in T t,max we considered the former mesh size as appropriate. A similar convergence test was employed to estimate the optimal outer dimensions. The PARDISO direct solver was used with the implicit intermediate Backward Differentiation Formulas (BDF) time stepping method, where the intermediate configuration was chosen to fix the initial and maximum steps of the solver, in this case 0.01 s and 1 s, respectively.