Void distribution and bubble motion in bubbly flows in a 4×4 rod bundle. Part II: numerical simulation

Abstract

Numerical simulations of bubbly flows in a four by four rod bundle are carried out using a multi-fluid model to examine effects of the numerical treatment of phase distribution and drag model. The transport equations of bubble number density and void fraction are used as the continuity equation of the gas phase. Two drag models are tested: one of them accounts for the bubble deformation (aspect ratio), whereas the other does not. The rod diameter, the rod pitch and the hydraulic diameter of the rod bundle are 10, 12.5 and 9.1 mm, respectively. The gas and liquid volume fluxes are J_G = 0.06 m/s and J_L = 0.9 and 1.5 m/s, respectively. The bubble diameter ranges from 1 to 5 mm. Comparisons between the numerical and measured data show that (1) the restriction on bubble lateral motion due to the presence of rods can be taken into account by using the transport equation of bubble number density, whereas that of the void fraction cannot deal with the restriction and causes large errors in the distribution of void fraction and (2) the reduction in the bubble-relative velocity near the wall is predictable by using the drag model accounting for the bubble deformation effect.

Keywords:

• bubbly flow
• rod bundle
• CFD
• two-phase flow
• multi-fluid model
• drag force

1. Introduction

Numerical methods for simulating two-phase flows in rod bundles of nuclear power reactors play an important role in designing fuel assemblies and evaluating the safety of nuclear power plants. Although multi-dimensional simulations of flows in a rod bundle have been carried out [1–8], their validation is still a difficult task due to the lack of detailed experimental data such as distributions of void fraction and velocities. In our previous study [9], we therefore carried out experiments on air–water bubbly flows in a four by four rod bundle to obtain experimental data applicable to the validation of the numerical methods. The experiments showed that the presence of rod walls restricts the motion of bubbles near the walls and affects the bubble shape. The sphere–volume equivalent bubble diameter, d, observed in the experiments ranged from 1 to 5 mm. These bubbles moved from the center of a subchannel toward the rod wall due to a shear-induced lift force [10]. The bubbles in the near wall region changed from oblate spheroidal to spherical shape due to the wall effect and mainly accumulated in the region d/2 away from the rod surface. Bubbles were apt to take distorted prolate spheroidal shape in the gap between the rods. The bubble relative velocities in the near wall and gap regions were smaller than those in the subchannel center region. This implies that the drag force increases as the bubble shape changes from oblate to prolate spheroids.

Since these wall effects on bubble shape, bubble lateral motion and bubble velocity have not been taken into account in multi-fluid models, they may fail in predicting the above-mentioned experimental facts.

Most of the incompressible multi-fluid models use the volume conservation of the gas phase for the continuity equation to predict the time evolution of the phase distribution. When the bubble size is larger than the computational cell size, the volume fractions in several cells surrounding the center of the bubble are non-zero though the bubble center is not located in these cells. Therefore, the interfacial momentum transfer, such as the drag and lift forces, in the momentum conservation equation must be modeled as forces acting on the volume fraction in each cell. However, in available constitutive equations, it has been modeled as forces acting on the bubble center. Hence the volume conservation and momentum conservation equations are inconsistent with the models of forces acting on the bubble center [11]. Tomiyama et al. [11] demonstrated that this inconsistency causes large errors in multi-fluid simulations. They also pointed out that since the forces act on the bubble center, they must be used to calculate not the gas phase-averaged velocity but the bubble velocity and the bubble velocity should be used to transport the bubble number density. They therefore proposed to use the transport equation of the bubble number density as the continuity equation of the gas phase, instead of the volume conservation equation, to recover the consistency. In their method, the shape of a bubble larger than the cell size is taken into account using a scalar function, which is called the shape factor. The use of this method would enable us not only to recover the consistency but also to account for the bubble shape deformation and the restriction of bubble lateral motion near the rod walls.

Tomiyama et al. [12,13] pointed out that the bubble aspect ratio plays an important role in the bubble velocity and proposed a drag model as a function of the bubble aspect ratio. This drag model was applied to single bubbles rising through stagnant water in a two by two rod bundle, and the estimation agreed well with their experimental data [14]. The velocity reduction in the near wall region observed in our recent experiments [9] might therefore be predicted by making use of this drag model.

In this study, numerical simulations of bubbly flows in the four by four rod bundle were carried out using a multi-fluid model utilizing the transport equation of the bubble number density and the drag model accounting for bubble deformation effect to numerically examine the above expectations. Distributions of liquid and gas velocities and void fraction of air–water bubbly flows obtained in the previous experimental study [9] were used for comparison with predictions.

2. Numerical method

2.1. Field equations

The gas and liquid are assumed isothermal incompressible fluids. Bubbles are classified into N classes of different bubble diameters. The continuity equations of the liquid and gas phases are given by (1) (2) where α is the volume fraction, V is the velocity, t is the time and the subscripts L and G denote the liquid and gas phases, respectively, and m represents the bubble class, m = 1, 2, …, N. The volume fractions satisfy (3) The momentum equations of the liquid and gas phases are given by (4) (5) where ρ is the density, P is the pressure, g is the acceleration of gravity, τ is the viscous stress tensor, τ^Re is the Reynolds stress tensor, and M is the interfacial momentum transfer. The viscous stress tensor is given by (6) where μ is the viscosity. The momentum transfer, M, consists of the drag, M_D, lift, M_L, virtual mass, M_V, turbulence dispersion, M_T, and wall forces, M_W, i.e. (7)

2.2. Time evolution of phase distribution

Most of the incompressible multi-fluid models use the volume conservation of the gas phase as the continuity equation and solve Equation (2) to obtain the time evolution of phase distribution. Hereafter this conventional method is referred to as the α-based method. Tomiyama et al. [11] suggested that the transport equation of the bubble number density, n_G, should be used instead of Equation (2) when dealing with bubbles larger than computational cell sizes. Substituting α_G = n_Gθ_B, where θ_B is the bubble volume, into Equation (2) yields (8) To deal with bubbles larger than the computational cell size, the shape factor, f, which is a scalar function representing a fraction of the bubble volume included in each computational cell as shown in Figure 1, is introduced. The volume fraction in a cell, IJK, required for computing M is calculated by summing up the contribution from the surrounding cells: (9)

Figure 1. Shape factor.

where θ is the computational cell volume, the subscripts I, J and K are the cell indexes in the x, y and z directions, and L, M and N are the indexes of cells surrounding the cell IJK. The subscript LMN ⇒ IJK denotes the contribution from the cell LMN to IJK. This method is referred to as the n-based method. Details of the latter method are given in [11].

2.3. Restriction on bubble lateral motion due to rod and channel walls

2.3.1. Wall force in α-based method

The following wall force [15] is utilized in the α-based method to account for the geometrical restriction to the bubble lateral motion due to the presence of rod and channel walls: (10) where C_W is the wall force coefficient, L is the distance from the wall, D_e is the hydraulic equivalent diameter of a subchannel, t_W is the unit tangent to the wall, and n_W is the unit normal to the wall. The following correlation (Hosokawa & Tomiyama [16]) is used for C_W: (11)

where Re, Eo and M are the bubble Reynolds number, Eötvös number and Morton number defined by (12) (13) (14) where g is the magnitude of g, σ is the surface tension and Δρ = ρ_L − ρ_G.

2.3.2. Geometrical restriction on bubble motion in n-based method

In the n-based method, the geometrical restriction to the bubble lateral motion due to the presence of rod and channel walls can be accounted for by preventing the migration of bubbles into the near wall region [17]. Figure 2 shows areas, in which the bubble center can lie, i.e. n_G > 0, in inner and side subchannels. Fluxes of n_G through the boundary between the shaded and white regions are forced to be zero, which ensures that n_G = 0 in the white (near wall) region throughout a simulation. Since bubbles near the rod wall were apt to be spherical [9], the minimum distance between the wall and the bubble center is set at d/2. Bubbles in the gap region took distorted prolate shapes and the maximum aspect ratio was 1.35, where the aspect ratio E is the ratio of the maximum vertical dimension to the maximum horizontal dimension of a bubble. Hence the minimum distance between the wall and the bubble center in the gap is set at the half of the minor axis of the bubble of the maximum aspect ratio, E = 1.35.

Figure 2. Restriction of bubble lateral motion in n-based method. (a) channel box and (b) inner and side subchannels.

2.4. Drag model

The drag force is given by (15) where C_D is the drag coefficient. The following drag correlation [18] has been often used in bubbly flow simulations [19,20]: (16) It should be noted that this correlation was originally developed for single bubbles rising through infinite stagnant liquids. The shapes of bubbles in the rod bundle, however, differ from those in the infinite liquid due to the wall effect. Therefore, Equation (16) would not be appropriate for bubbles in rod bundles.

In contrast to drag models without the bubble aspect ratio, drag models accounting for the effects ofbubble aspect ratio have been rarely proposed [12,13,21,22]. A drag model proposed by Tomiyama et al. [12,13] takes into account the effects of the bubble aspect ratio on the drag force: (17)

where F and G are functions of E given by (18) (19) Though this model was also proposed for free-rising bubbles, its applicability to single bubbles in a rod bundle has been confirmed [14]. Hence, Equation (17) can be utilized to account for the wall effect on drag force by substituting an appropriate value of E.

The drag multiplier, 1 + 0.55Sr², which accounts for the effect of liquid velocity gradient about a bubble [23], is applied to C_D, where Sr is the dimensionless shear rate defined by (20) where ω is the magnitude of the liquid velocity gradient. The effect of bubble swarm is also taken into account by making use of the multiplier α_G^3−2φ to C_D, where φ = 1.75 [24].

2.5. Other closure relations

The lift force is given by (21) where C_L is the lift coefficient. The following C_L correlation proposed by Tomiyama et al. [10] is used: (22) where Eo_H is the Eötvös number using the maximum horizontal dimension of a bubble as the characteristic length, i.e. Eo_Hm = Δρ_Lga_m²/σ and a_m = d_mE_m^−1/3. Lucas and Tomiyama [25] investigated the applicability of the above correlation to various bubbly flows and confirmed that it is applicable not only to air–water bubbly flows but also to vapor–water bubbly flows at high pressures.

The virtual mass force is given by (23) where C_V is the virtual mass coefficient and is given as a function of the aspect ratio [13,26]: (24) The turbulence dispersion force is given by [27] (25) where C_T is the turbulence dispersion coefficient. Since there are no available correlations of C_T for bubbles in rod bundles, the value of C_T was optimized by trial and error and is given in Section 2.6. The turbulence kinetic energy, k_L, of the liquid phase is given as a superposition of the shear- and bubble-induced turbulence kinetic energies, k^SI and k^BI: (26) The shear-induced turbulence kinetic energy is calculated using an anisotropic k–ϵ model (SZL model: Shih–Zhu–Lumley model) proposed by Shih et al. [28]. The bubble-induced energy is given by [29] (27)

The Reynolds stress tensor is modeled as a superposition of the shear- and bubble-induced turbulence stresses [30]: (28) The shear-induced Reynolds stress is calculated by using the SZL model. The bubble-induced turbulence stress is evaluated as (29) where the first and second terms are based on the Sato's model [30] and the Lopez de Bertodano's anisotropic model [29], respectively, and A^BI is the anisotropy matrix given by (30)

2.6. Numerical condition

The computational domain is shown in Figure 3. Due to the symmetry of the flow [9], the one-fourth of the four by four rod bundle was simulated. The Cartesian coordinates were used and the cell widths, Δx, Δy and Δz, in the x, y and z directions were 0.25, 0.25 and 0.75 mm, respectively. The number of cells assigned to the x, y and z directions were 108, 108 and 10, respectively. The rod diameter and pitch were 10 and 12.5 mm. The number of cells assigned to the rod diameter was therefore 40. As shown in the Appendix, we confirmed in our preliminary study that a good enough prediction for a single-phase flow in a rod bundle is obtained by using the Cartesian coordinates and this spatial resolution. The rod and channel box were simulated by imposing the no-slip condition to cells locating at their walls. The top and bottom boundaries were periodic. The time step, Δt, was 7.5×10⁻⁵ s. The solution algorithm can be found in [11]. The fluid properties of gas and liquid phases were those of air and water at room temperature and atmospheric pressure as in the experiments [9].

Figure 3. Computational domain.

Four cases, S1, S2, T1 and T2, were tested. Cases S1 and S2 were single-phase flows of J_L = 0.9 and 1.5 m/s, respectively. Cases T1 and T2 were two-phase bubbly flows at (J_L, J_G) = (0.9, 0.06 m/s) and (1.5, 0.06 m/s), respectively, where J_L and J_G are the liquid and gas volume fluxes.

In view of the experimental data [9], bubbles were classified into four classes (N = 4), whose diameters were 1.4, 2.2, 3.2 and 4.2 mm, respectively. The values of C_T for these bubble classes are given in Table 1.

Table 1. Turbulent dispersion coefficient C_T.

	CT in n-based	CT in α-based
d (mm)	method	method
1.4	0.50	1.0
2.2	0.40	0.50
3.2	0.02	0.50
4.2	0.01	0.50

The bubble shapes in the four by four rod bundle were also measured in [9]. Bubbles near the rod wall were apt to be spherical, E = 1, whereas E > 1 for bubbles in the gap due to the strong wall effect. Tomiyama et al. [14] measured the aspect ratios of bubbles in a two by two rod bundle and showed that E ≈ 0.75 for bubbles in the center of the subchannel. The aspect ratio was therefore given by the following interpolation function: (31)

3. Results and discussion

3.1. Liquid velocity in single phase flow

The numerical method was first validated for single-phase flows in the rod bundle. Predicted liquid velocities along the y-axis at x = 27 mm in Cases S1 and S2 are compared with the measured data [9] in Figure 4. Although some errors appear in the gap regions, the predictions are reasonable. The root-mean-square (rms) errors defined by (32)

Figure 4. Streamwise liquid velocity in single-phase flows (x = 27 mm).

were 0.12 and 0.13 m/s in Cases S1 and S2, respectively, where N_m is the number of data and the superscripts m and p denote the measured and predicted data. Baglietto and Ninokata [3] reported that a secondary flow in a single-phase turbulent flow in a tight-lattice rod bundle can be well predicted by increasing a model parameter controlling the intensity of anisotropy in the SZL model. Better predictions especially in the gap regions could be obtained by tuning the parameter in the present simulations. Therefore the errors in the gap regions can be attributed to the weakness of secondary flow due to the anisotropy. However, the main purpose of this study is to examine the effects of the numerical treatment of bubble motion and drag model. The present SZL model is acceptable for this purpose, and therefore, we do not make further discussion on the optimization of the parameters in the SZL model.

3.2. Effects of numerical treatment of phase distribution

The measured void fraction in Case T1 is shown in Figure 5(a). As described in Section 1, bubbles are apt to migrate from the subchannel center to the wall due to the lift force. The void fraction therefore increases from the center of an inner subchannel (y ≈ 14 and 27 mm) toward the rod. Since bubbles rarely enter into the gap region and they do not stay in the gap due to cross flows between subchannels, the void fractions in the gap regions (y ≈ 8 and 21 mm) are very low. Predicted α with the α-based method is compared with the measured data in the figure. Equation (17) was used for C_D in this simulation. The predicted void fraction in the inner subchannel is too low and the low void fraction in the gap is not reproduced. This is because the wall repulsive force cannot prevent the migration of bubbles toward the rod gap due to the lift force. To the contrary, the prediction with the n-based method qualitatively agrees with the measured data and the low void fraction in the gap is well reproduced, which implies that the geometrical restriction of the bubble lateral motion works well.

Figure 5. Void fraction and liquid and gas velocities in Case T1 (x = 27 mm). (a) Void fraction and (b) Liquid and bubble-relative velocities.

The measured and predicted V_L and the bubble-relative velocity, V_R ( = V_G − V_L), are compared in Figure 5(b). The prediction of V_L with the α-based method seems better than that with the n-based method. However the α-based method gives large errors in V_R in the whole domain, so that the prediction of V_R with the n-based method is much more reasonable.

As shown in Figure 6, the differences between the predictions of the α- and n-based methods are not so remarkable in Case T2 compared to those in Case T1. The volume fractions in the near wall and the gap regions are however better predicted with the n-based method than with the α-based method. The root-mean-square errors in the above predictions of α, V_L and V_R are summarized in Table 2.

Table 2. Root mean square errors e_rms in predicted α, V_L and V_R.

			erms in VL	erms in VR
Case	Scheme	erms in α	(m/s)	(m/s)
T1	n-based	0.037	0.137	0.035
	α-based	0.053	0.104	0.11
T2	n-based	0.016	0.191	0.17
	α-based	0.022	0.181	0.20

Figure 6. Void fraction and liquid and gas velocities in Case T2 (x = 27 mm). (a) Void fraction and (b) Liquid and bubble-relative velocities.

The measured area-averaged void fraction in the corner subchannel was 0.030 in Case T1 and those predicted using the n- and α-based methods were 0.030 and 0.105, respectively. These results imply that the geometrical constraint in the n-based method works well in the corner subchannel. However, in Case T2, both predictions of the n- and α-based methods, i.e. 0.0551 and 0.0696, largely differed from the measured value of 0.016. More precise treatment for the wall effect on the bubble aspect ratio may be required for more accurate predictions.

3.3. Effects of drag model

Simulations of Cases T1 and T2 are carried out using the n-based method with the drag models, Equations (16) and (17), to examine the effects of C_D models. The differences between the two models are not large in the predicted α in Case T1. However the bubble-relative velocity obtained using Equation (16) shows large errors as shown in Figure 7(a). Bubbles in the rod bundle are apt to be more spherical or prolate spheroidal compared to those in infinite liquids. According to Tomiyama et al. [12,13], the bubble rise velocity decreases with increasing E. Therefore the bubbles in the rod bundle are slower than those in infinite liquids. Since the C_D model without E cannot account for this effect, it gives small C_D for bubbles in the rod bundle. This is the main reason why V_R predicted with this model is higher than the measured data. The C_D model, Equation (17), with E gives better predictions.

Figure 7. Effects of C_D model on V_R (n-based) (x = 27 mm). (a) Case T1 and (b) Case T2.

As shown in Figure 7(b), the predictions with Equation (16) in Case T2 also show larger errors in V_R than those with Equation (17) due to the same reason.

The above discussion has made clear that the transport equation of the bubble number density and the drag model accounting for bubble deformation should be used to obtain reasonable predictions for bubbly flows in a rod bundle; otherwise the simulation fails in predicting the bubble motion especially in the near wall regions. Although the present experimental and numerical conditions differed from those in actual rod bundles, i.e. the present flows were the adiabatic air–water bubbly flows at relatively low liquid velocities, the above-mentioned fact is of great importance when predicting bubbly flows in rod bundles.

4. Conclusion

Hosokawa et al. [9] measured distributions of void fraction and gas and liquid velocities in a four by four rod bundle to provide a database for validation of numerical methods. The experiments showed that the presence of rod walls restricts the motion of bubbles near the wall, which strongly affects the distribution of void fraction. Bubble shapes are apt to change from oblate to distorted prolate spheroidal as bubbles approach the rod wall or the rod gap, which results in the reduction of bubble-relative velocity. It was expected that the former and latter effects can be predicted using a multi-fluid model based on the transport equation of the bubble number density [11] and a drag model accounting for the bubble deformation effect [12,13], respectively. In this study, numerical simulations of bubbly flows in the four by four rod bundle were therefore carried out using a multi-fluid model to numerically examine the above expectation. The main conclusions obtained from comparisons between the numerical and measured data are as follows.

1. The geometrical restriction on bubble motion due to the rods can be taken into account by using the transport equation of the bubble number density, whereas that of the void fraction cannot account for the restriction and causes large errors in a predicted distribution of void fraction.

2. The reduction in the bubble-relative velocity near the wall is predictable by using the drag model accounting for the bubble deformation effect.

It is therefore recommended to use the transport equation of the bubble number density as the continuity equation of the gas phase and the drag model accounting for the bubble deformation effect in simulations based on multi-fluid models. A question, however, arises as to how we should model the aspect ratio of bubbles in a rod bundle, especially for largely deformed bubbles in the gap region. Modeling the wall effect on bubble shape would be an important issue to obtain more accurate predictions of bubbly flows in rod bundles.

Nomenclature

ABI:	=	anisotropy matrix
CD:	=	drag coefficient
CL:	=	lift coefficient
CT:	=	turbulence dispersion coefficient
CV:	=	virtual mass coefficient
CW:	=	wall force coefficient
d:	=	sphere-volume equivalent bubble diameter [m]
De:	=	hydraulic equivalent diameter of subchannel [m]
erms:	=	root mean square error
E:	=	aspect ratio
Eo:	=	Eötvös number
EoH:	=	Eötvös number defined using bubble major axis
f:	=	shape factor
F, G:	=	functions of E
g:	=	acceleration of gravity (m/s2)
g:	=	magnitude of g (m/s2)
h:	=	function of EoH
JG:	=	gas volume flux (m/s)
JL:	=	liquid volume flux (m/s)
k:	=	turbulence kinetic energy (m2/s2)
L:	=	distance from channel wall (m)
M:	=	interfacial momentum transfer (N/m3)
MD:	=	drag force (N/m3)
ML:	=	lift force (N/m3)
MV:	=	virtual mass force (N/m3)
MT:	=	turbulent dispersion force (N/m3)
MW:	=	wall force (N/m3)
n:	=	bubble number density (m−3)
N:	=	Total number of bubble classes
nW:	=	unit normal to wall
P:	=	pressure (Pa)
Re:	=	bubble Reynolds number
Sr:	=	dimensionless shear rate
t:	=	time (s)
tW:	=	unit tangent to wall
V:	=	velocity (m/s)
x, y, z:	=	Cartesian coordinates (m)
α:	=	void fraction
Δt:	=	time [s]
Δρ:	=	density difference (Δρ = ρL − ρG) (kg/m3)
Δx, Δy, Δz:	=	cell size in the x, y and z directions (m)
θ:	=	cell volume (m3)
θB:	=	bubble volume (m3)
μ:	=	viscosity (Pa·s)
τ:	=	viscous stress tensor (Pa)
τRe:	=	Reynolds stress tensor (Pa)
ρ:	=	density (kg/m3)
σ:	=	surface tension (N/m)
φ:	=	model constant for bubble swarm effect on drag force
ω:	=	magnitude of liquid velocity gradient (s−1)
Subscripts	=
G:	=	gas phase
L:	=	liquid phase
m:	=	mth bubble class
Superscripts	=
SI:	=	shear-induced
BI:	=	bubble-induced
m:	=	measured data
p:	=	predicted data

Acknowledgements

The authors would like to express their sincere thanks to Mr Yuya Muramatsu for his effort in the experiments and numerical simulations.

Additional information

Funding

This work has been supported by the Japan Society for the Promotion of Science (JSPS) (grants-in-aid for scientific research (B)[grant number 24360070])

Appendix

Curvilinear coordinates can precisely deal with the rod wall. They however often suffer from the skewness and elongation of computational cells and the numerical treatments for the closure relations of interfacial momentum transfer are to be complicated in curvilinear coordinates. Cartesian coordinates could therefore be a better choice, provided that the jaggy boundary representation for rod walls does not deteriorate predictions.

Effects of coordinate systems and spatial resolutions on the predicted motion of single drops in vertical round tubes were investigated using three different coordinate systems, i.e. boundary-fitted, cylindrical and Cartesian coordinate systems, in our previous study [31]. The results confirmed that the Cartesian coordinate can give good enough predictions for Taylor drops that form thin liquid films between the drop and the pipe wall even with only 12 cubic cells along the pipe diameter.

Effects of Cartesian coordinates on a prediction of a flow in a rod bundle were also preliminary examined. Hosokawa et al. [32] simulated a single-phase flow in a two by two rod bundle by using a boundary-fitted coordinate method implemented in OpenFOAM. The horizontal cross-section of the rod bundle is shown in Figure A1(a). Since the flow was symmetric, the computational domain was one-fourth of the whole flow domain. The area-averaged liquid velocity was 2.0 m/s. The boundary-fitted coordinate system used in [32] is shown in Figure A1(b), in which the number of computational cells in the horizontal plane is 8000.

Figure A1. Simulation of single-phase flow in two by two rod bundle. (a) two by two rod bundle (flow is symmetric with respect to dotted-lines), (b) boundary-fitted cells and (c) predicted streamwise liquid velocities at x = 25 mm.

The domain was discretized using uniform cubic cells in our preliminary study. The spatial resolution was the same as in the simulation of the four by four rod bundle, i.e. 40 computational cells along the rod diameter. The number of cells in the fluid region was only about 3/20 of that in the boundary-fitted coordinate simulation. Predicted streamwise liquid velocities are shown in Figure A1(c). In spite of the lower spatial resolution, the prediction with the Cartesian coordinates agrees well with the prediction of boundary-fitted coordinates. Therefore, with the spatial resolution of 40 cells along the rod diameter, the effects of jaggy boundary representation for rod walls are sufficiently small for the prediction to converge to that with the boundary-fitted coordinates.