Numerical Simulation of a Domain-Tessellation Pattern on a Spherical Surface Using a Phase Field Model

Show more

Received 15 February 2016; accepted 16 April 2016; published 19 April 2016

1. Introduction

The phase field model has been developed to simulate microstructures in materials, and various complicated patterns such as dendrite, cellular, lamellar, and polycrystalline structures have been successfully regenerated [1] - [4] . In addition, the model has been widely used not only in materials science but also in various fields such as topological optimization [5] [6] . We also applied the phase-field method to analyze the space-filling structure in two- and three dimensions [7] [8] . In the 2-D problem, an infinite plane can be divided into regularly repetitive patterns with typical symmetric unit shapes such as regular triangles, squares, and regular hexagons. In 3-D space, the repetitive space division is achieved using a regular truncated octahedron or a rhombic dodecahedron, which are also defined as the Voronoi tessellations of body-and face-centered cubic lattice points, respectively. Especially, the former is known as the Kelvin cell and is considered to have the best geometrical properties, i.e., the highest symmetry and smallest surface area. The author has further demonstrated 3-D space division by numerical simulation using the phase field model [9] , and successfully obtained these patterns. Moreover, it has been found that a significant advantage of the truncated octahedron is its stability. Microstructure modeling and simulation of foam or porous materials are some of the expected applications of this model.

Now, our next target is region division on a curved surface. First, we focus on a spherical surface. In a flat plane, the most preferred unit is a regular hexagon because of its symmetry. However, region division on a spherical surface is not straightforward; it is impossible to divide the spherical surface into equally shaped hexagons. This is revealed by reviewing polyhedra or platonic solids. For instance, a regular dodecahedron consists of 12 pentagons, and a regular icosahedron consists of 20 regular triangles. This geometrical feature means that the surface division with regular pentagons or triangles is possible if slight distortion along the curved surface is allowed. However, these are only the exceptions, and it is well known that regular polyhedra are realized only for the tetra-, hexa-, and octahedron in addition to the dodeca- and the icosahedron. Similar types of problem have been investigated in geophysics and meteorology and developed as the discrete global grid system [10] [11] , where the surface of the earth is divided into domains. Regular domain tessellation of curved surface has also been an intensive research subject in computer graphics [12] . As geometrical schemes, Voronoi tessellation and Delaunay triangulation are quite useful, and the computational algorithm has been improved [13] [14] . However, these techniques are intrinsically based only on geometrical features. Other candidates can be found in the molecular structures of fullerene, though they consist of two types of unit cells. For instance, a C60 consists of 12 regular pentagonal and 20 regular hexagonal faces. This type of semi-regular polyhedron, or Archimedean solid, helps to identify a better division; however, these cases also provide a limited solution. Therefore, we are motivated to find the space-dividing pattern by computer simulation, and apply the above mentioned phase-field method to this end, since the model can capture not only geometrical but also various kinds of physical and chemical features. Phase-field equations on the 2-D space are extended to the spherical surface. The fundamental equations to be solved are differential equations, and hence, numerical integration is required. The computational scheme as well as the phase-field equations is explained in Section 2, and the simulation results are presented in Sections 3 and 4.

2. Fundamental Equations and Numerical Scheme

2.1. Phase-Field Equations

The multi-phase-field model is applied in this study. This model was originally proposed by Steinbach et al. to simulate phase transformation among multiple phases [15] and subsequently adopted to polycrystalline modeling [16] . In the present study, area division is achieved by setting the initial domain nuclei and growing the nuclei using the phase-field equations. Each domain is represented by a multi-phase-field variable (i = 1 ~ N), where i indicates the i-th domain, i.e., the value of is 1 in domain i and 0 in the other locations. This variable changes continuously, but steeply, from 1 to 0 in the domain boundary, and the division pattern is represented by the boundary lines. The phase-field equations used in this paper are the following:

(1)

where

(2)

(3)

(4)

and

(5)

Here, , , and are based on Steinbach’s original model, though the notation is modified and the equations and parameters are standardized to dimensionless values. m_{ij}, a_{ij}, w_{ij}, and Δf_{ij} are parameters depending on the combination of i and j, and n is the number of phases on site. The last term, , is an additional term to control the area to be uniform, where S_{i} is the area of the i-th domain and is an adjusting coefficient. Note that the area S_{i} has a relative value with respect to the total surface area of the large sphere. This term is introduced to prevent domain overgrowth and to keep all domain area equal, because domain coarsening occurs due to the nature of the phase field model [8] .

2.2. Numerical Integration

The multi-phase-field Equation (1) is solved numerically using the finite difference method (FDM). Grids are usually arranged on a rectangular lattice in 2-D FDM, but such gridding is not applicable on a curved face. Especially for a spherical surface, rectangles are unable to fill the surface completely. For example, the grid arrangement based on latitude and longitude division unavoidably creates distortion around the polar area. Triangular division has greater flexibility, and it is possible to fill the surface with mostly equally spaced grid lines, although precisely equal spacing is impossible. Even though the difference is negligible, grid symmetry cannot be realized. In a planar problem, each grid in the regular triangular division has six neighboring grids, but some of the grids on a spherical surface must have five neighbors, according to Euler’s theorem on the polyhedron. In this study, we create a triangular lattice by accepting this irregularity, and solve the conventional diffusion equation using the grids to check their validity before applying them to the phase-field equation.

Figure 1 shows the grid arrangement; the sphere under consideration, with radius R, is approximated by the regular dodecahedron (Figure 1(a)) inscribed in the sphere. The 20 large triangle faces are divided into small regular triangles, and then the sub-grids are projected onto the spherical surface along the radial lines from the center of the sphere. Figure 1(b) represents the grid disposition when each dodecahedron edge is divided into 40 grids; each grid is depicted by a small sphere and green color represents those originally disposed on the triangle edges. The projection involves distortion of the grid interval, and the distortion is adjusted to be approximately the same for all intervals. These operations do not produce precisely equal spacing, and deviation is unavoidable; Figure 1(c) represents the distribution of the average distance between the neighboring grids. Blue spots correspond to the vertices of the dodecahedron; around them, the points have only five neighbors, whereas other points have six neighbors. The distribution is actually gradual, and the color range is so small that the anticipated computational error is small enough for the current problem.

The discretized equation for the simple diffusion equation is represented below:

(6)

(a) (b) (c)

Figure 1. Grid arrangement on a spherical surface based on regular dodecahedron, and distribution of the local grid density expressed by the average distance between the neighbor grids. (a) Regular dodecahedron, (b) grid arrangement, (c) local grid density.

where n is the number of neighboring grids, Δt is the time increment, and Δr_{ij} is the distance between the central i-th grid and the neighboring j-th grid. This is a type of isotropic integration and is the case of forward difference with respect to the time increment. This form is reduced for usual FDM calculations with a regular square lattice when n = 4 and Δr_{ij} = Δx. In the current stage of our study, the influence of the curvature on the numerical calculation is not imposed.

2.3. Validation of Current Method

The presented numerical integration scheme is validated by solving a conventional diffusion problem. The non- dimensional heat conduction equation is solved for an initial temperature T = 0.0 in entire domain. Then, the temperature in one of the polar domains, i.e., the area where the normal vector is inclined up to five degrees in the positive direction of the z axis, is fixed at T = 1.0. Figure 2(a) represents the temperature distribution at the 3000-th and 10,000-th time steps. Heat conduction occurs radially from the polar domain and the isothermal lines form concentric circles; Figure 2(a) fairly represents this tendency. Figure 2(b) shows the temperature variation at selected points. The plot points A, B, and C follow the meridian lines inclining 10˚, 45˚, and 60˚ from the x axis, and the numbers following the alphabets represents the angle from the polar direction, as represented in the inset in Figure 2(b). The curves at the same latitude show mostly identical variation, which indicates that isotropic radial heat conduction is reasonably obtained by the present method. Actually, these curves are not exactly identical, but the difference is considered to be in the allowable range for the present qualitative study.

3. Simulation Results without Size-Controlling Term

3.1. Model and Conditions

A simple situation is first considered to validate the present method. Grids shown in Figure 1(b) is used for numerical integration, where the total number of grids is approximately 16,000. As shown in Figure 1(c), the interval between neighboring grids shows a slight deviation. The maximum interval is set to be Δx_{max} = 0.50, and the minimum interval is then Δx_{min} = 0.33. The sphere radius R is set as R = 50Δx_{max}.

The original phase and growing phases are expressed by and (i = 1 ~ N), respectively, and the parameters are listed in Table 1. The parameters with the subscript ij generally consist of independent values be-

(a) (b)

Figure 2. Temperature distribution and variation at selected points according to the preliminary calculations. (a) Temperature distribution, (b) variation of temperature.

Table 1. Simulation parameters applied in this paper.

tween the i-th and j-th phases; however, in this study, identical values are applied for each combination. Parameter Δf_{ij} represents the difference in free energy between two phases, and this value is assumed non-zero only for the interaction between the original phase and growing phases. First, the effect of the size-controlling term is excluded to observe the stability of the formed pattern. Then, to investigate the instable cases, this term is turned on by applying a non-zero value to. This value is given only between the growing phases, so that the growths against the original phase are not affected.

The nuclei of the growing phases are set as follows: Case 1: intersection points of the spherical surface with the x, y, and z axes (6 points). Case 2: intersection with the equivalent eight directions denoted by vector (±1, ±1, ±1) (8 points). Case 3: Case 1 plus Case 2 (14 points). The phase-field value is assigned for the grids within radius r_{n} = 0.01 from the i-th nuclei point and is assigned for the other grids.

3.2. Results for Case 1

The simulation results for Case 1 are shown in Figure 3. The color indicates the value of. This value

equals 0 in the original phase (the region where) since i = 0 is not included in the equation; it equals 1 in the growing phase, where one of the values of equals 1; and it takes value between 0 and 1 at the boundary domain. These values are colored in blue for 0, red for 1, and continuous gradation between 0 and 1. Each illustration shows a snapshot projected on the x-y plane, unless otherwise noted, at the time steps noted below each figure.

Figure 3(a) represents the initial condition. The nuclei are set at the positive and negative directions along the x, y, and z axes. Theoretically, these phases should grow isotropically, according to the characteristics of the phase-field equations, because anisotropic properties are not included. Figure 3(b) and Figure 3(c) show the concentric growth from all nuclei, and hence, it can be concluded that the numerical scheme introduced here is valid. The growth rates from each nucleus are also the same, and the domains collide at four points, in the visible range, at the same time. Phase transformation is completed in the entire surface by the 1500-th time step, and the domain boundary is formed, as shown in Figure 3(d). Continuing the computation, no change is observed by the 5000-th time step, as shown in Figure 3(e) and Figure 3(f), and this shape is considered stable. Figure 3(g) shows the projection on the x-z plane at the 5000-th time step, which clearly indicates that each domain has the same morphology both in size and shape, and each domain has a barrel shape with a curved boundary and four triple junctions at the corners. Figure 3(h) is the view from another angle, and the color range is modified to emphasize the triple junctions. In a planar calculation of the current phasefield model, it is proved that the hexagonal space division is stable. Each hexagon consists of six triple junctions with 120˚ angle at the corner and six straight equal-length edges, i.e., 120˚ junction and straight boundary is the stable state for the present phase

(a) (b) (c) (d)(e) (f) (g) (h)

Figure 3. Variation of the domain shape for Case 1. Color indicates the value. (a) 0 step, (b) 500 steps, (c) 1000 steps, (d) 1500 steps, (e) 2000 steps, (f) 5000 steps, (g) x-z steps, (h) rotated angle with enhanced color.

field model. In this simulation, as shown in Figure 3(h), the triple junction locally has a 120˚ angle. This stability is considered to yield total equilibrium in the macroscopic morphology over the energy loss due to the curved boundaries.

3.3. Results for Case 2

The simulation results for Case 2 are presented in Figure 4. Free growth from the nuclei occurs properly. The first collisions are observed on the x-y, y-z, and z-x planes, and straight boundaries are formed. Phase transformation is completed by the 1500-th time step, and the cross-shaped boundaries are generated, as shown in Figure 4(c) and Figure 4(d). This state is stable; no change is observed by the 5000-th time step, as shown in Figure 4(e) and Figure 4(f). Moreover, a view on the x-z plane shows the same shape in Figure 4(g). Figure 4(h) shows another view from a rotated angle. Each phase region is surrounded by three boundaries, producing a triangular domain. The border lines are, of course, globally curved along the spherical surface, but here, we consider the line along the maximum circle of the sphere to be straight. In this state, four domains intersect perpendicularly (with 90˚ angels) at a point to make a quad-junction. As mentioned in Section 3.2, the phase field model yields a stable state for a 120˚ junction with a straight boundary. In Case 2, the junction is not stable, but the straight edge is preferable. In the planar calculation, a similar tendency is observed; the square division with a 90˚ quad-junction shows meta-stability. On a spherical surface, the triangle domain is considered equivalent to the square division on the flat plane. Since this state is not stable, but metastable, a slight fluctuation causes the state to collapse. Actually, it is confirmed that this pattern is destroyed when a slight deviation in the initial nuclei position is provided, and the results are shown in Figure 5. Almost the same pattern is generated by the 1500-th time step (Figure 5(a) and Figure 5(b)), but the quad-junction is soon segregated into two triple junctions, as shown in Figure 5(c). Then, the domain shape is globally modified and a barrel-shaped pattern similar to Case 1 is obtained, as shown in Figure 5(e), and no change was observed by continuing the calculation. It was also confirmed that the pattern obtained for Case 1 was stable under the same type of fluctuations.

(a) (b) (c) (d)(e) (f) (g) (h)

Figure 4. Variation of the domain shape for Case 2. Color indicates the value. (a) 0 step, (b) 500 steps, (c) 1000 steps, (d) 1500 steps, (e) 2000 steps, (f) 5000 steps, (g) x-z steps, (h) rotated angle with enhanced color.

(a) (b) (c) (d) (e)

Figure 5. Variation of the domain shape with a slight deviation causing instability of the domain shape. (a) 1000 step, (b) 1500 steps, (c) 3000 steps, (d) 5000 steps, (e) 10,000 steps.

3.4. Results for Case 3

The simulation results for Case 3 are shown in Figure 6. The nuclei are set at 14 points, and circular growth is observed, as shown in Figure 6(b). The phase transformation is completed by the 1000-th time step. Unlike the previous two cases, the 14 sites are not all equivalent, and collisions of the growing domain occur at two types of sites. As a result, two types of domains are formed, as shown in Figure 6(c): one is quadrangle, whose center falls in the <100> directions, and another is triangle, whose center falls in the <111> directions. Such inhomogeneity brings instability in the phasefield model, and the formed pattern does not last long. The triangle domains shrink gradually and the quadrangle domains grow larger, and, finally, the triangle domains vanish by the 4000-th time step (Figures 6(d)-(h)). The preserved quadrangle domains adjust their shapes, and settle into the same pattern as that obtained in Case 1 by the 4000-th time step, as shown in Figure 6(h). Since this state is stable, no change is observed by continuing the calculation.

4. Results and Discussion with Size-Controlling Term

4.1. Simulation Condition

In the Case-3 simulation, the small triangular domains disappear soon after formation. This result reflects the natural tendency of the phasefield model, but this feature sometimes limits its application. In this section, simulations are demonstrated by imposing a condition to prevent the domain size from reducing. The other conditions and parameters are set the same as those in the previous section, except that an extra term is activated with parameter. The nuclei are set as in Case 3 (Section 4.2) and at random sites (Section 4.3).

4.2. Results for Case 3 with Size-Controlling Term

The simulation results are presented in Figure 7. The initial domain growth is identical to that presented in Figure 6 since the original phase is excluded from the sizecontrol. Domain tessellation is completed by the 1000-th time step, and triangular and quadrangular domainsare formed. Unlike the case in Figure 6, the triangular domains do not shrink thereafter, and the triangular shape is maintained, as shown in Figure 7(c) to Figure 7(g). Instead, the quadrangle domains change their shape; they are initially nearly square (the central domain in Figure 7(c)), but change to rhomboid by the 6000-th time step, as shown in Figures 7(d)-(g). Figure 7(h) is a view from another angle at the 6000-th time step, where the color range is adjusted to enhance the position of the domain junctions. Indeed, every junction is not a single quad-junction, but two closely spaced triple junctions. Accordingly, the domains that appeared triangular are actually hexagons composed of three long edges and three very short ones. As for the rhombic domain, there are two types of directions, which are arranged alternately to make a unique domain-tessellation pattern.

(a) (b) (c) (d)(e) (f) (g) (h)

Figure 6. Variation of the domain shape for Case 3. Color indicates the value. (a) 0 step, (b) 500 steps, (c) 1000 steps, (d) 1500 steps, (e) 2500 steps, (f) 3000 steps, (g) 3500 steps, (h) 4000 steps.

(a) (b) (c) (d)(e) (f) (g) (h)

Figure 7. Variation of the domain shape for Case 3 using the size-controlling term. Color indicates the value. (a) 0 step, (b) 500 steps, (c) 1000 steps, (d) 1500 steps, (e) 2500 steps, (f) 4000 steps, (g) 6000 steps, (h) rotated angle with enhanced color.

4.3. Results for Random Nuclei Arrangement

Finally, the pattern formation initiated with randomly arranged domain nuclei is simulated. Figure 8 shows a typical result obtained by 12 domain nuclei, whose initial state is shown in Figure 8(a). The domain growth is isotropic, and circular growth is observed in the initial stage (Figure 8(b)). The collision of the neighboring domains occurs intermittently, and an irregular shape is generated (Figure 8(c)). Quadrangles, pentagons, and hexagons are observed at this time. These domains vary in shape, so that balance is achieved between neighboring domains (Figures 8(d)-(g)). Finally, they settle into a steady state presented in Figure 8(h), where the color range is modified to enhance the domain boundaries. This pattern is constructed by 12 pentagons. The lengths of the edges of each pentagon are almost identical, and the five angles are also the same; hence, the spherical surface is divided into regular pentagons. This pattern resembles a regular dodecahedron, which has 12 regular pentagons as its surfaces. Notably, although the 3-D phase-field simulation generates a rhombic dodecahedron with 12 faces, its shape is not pentagonal but rhombic. This implies that the present model properly reproduces a stable state by not only the face area but also a totally favorable domain tessellation.

A variation in the domain size and area distribution is drawn in Figure 9. The domain size R denotes the distance between the furthest apexes in a pentagon, and the maximum and minimum values in the 12 pentagons are represented as R_{max} and R_{min}, respectively. The value length is expressed in Δx unit. The maximum and minimum areas S_{max} and S_{min} are also plotted in the same figure, where the area is expressed in the ratio to the average area of all domains. Variation of the domain area shows that each domain area increases by the 1000-th time step, when free domain growth is completed. Then, the maximum area starts decreasing, while the minimum area keeps increasing. This means that large domains shrink and small ones keep enlarging, and hence, the size adjustment works effectively. The domain size variation does not necessarily follow this trend if the shape is not the same even if the area is equivalent. In Figure 9, both the area and size become identical at around the 8000-th time step, revealing that the steady state is achieved. Actually, it was confirmed that no change in the pattern was observed by continuing the calculation further.

5. Conclusion

Numerical analyses were carried out to investigate the space-dividing pattern on a spherical surface. The phasefield model was applied to simulate the domain distinction and its growth. A simple method was introduced for the finite-difference calculation on a spherical surface, and its validity was confirmed by solving a conventional diffusion equation resulting in concentric isolines. Pattern formation simulation was subsequently demonstrated and unique patterns were obtained; a stable barrel-shaped pattern with six faces was formed, whereas triangular division with cross-boundaries was meta-stable. The size-controlling term allowed the reproduction of a wider

(a) (b) (c) (d)(e) (f) (g) (h)

Figure 8. Variation of the domain shape for random nuclei arrangement using the size-controlling term. Color indicates the value. (a) 0 step, (b) 500 steps, (c) 1000 steps, (d) 1500 steps, (e) 3000 steps, (f) 5000 steps, (g) 7500 steps, (h) 10,000 steps with enhanced color.

Figure 9. Variation of the maximum and minimum domain size and area.

variety of patterns, and a dodecahedral division with 12 pentagons was exhibited as an example. Consequently, the present method reveals to be valid and effective for pattern formation analysis. In this study, only a limited number of models were investigated, i.e., N = 6, 8, 14 for regular nuclei arrangement, and N = 12 for random arrangement. Nevertheless, many other patterns have been obtained for different values of N. These data will be analyzed in our next study based on geometrical theory in relation to the surface area and boundary or edge lengths. Furthermore, this method will be applicable for some practical purposes. For example, this method is expected to be utilized in materials design by coupling the geometrical features with interfacial energy and some other physical and chemical properties.

References

[1] Provatas, N. and Elder, K. (2010) Phase-Field Methods in Materials Science and Engineering. Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim. http://dx.doi.org/10.1002/9783527631520

[2] Wang, Y. and Li, J. (2010) Phase Field Modeling of Defects and Deformation. ActaMaterialia, 58, 1212-1235.
http://dx.doi.org/10.1016/j.actamat.2009.10.041

[3] Steinbach, I. and Shchyglo, O. (2011) Phase-Field Modelling of Microstructure Evolution in Solids: Perspectives and Challenges. Current Opinion in Solid State and Materials Science, 15, 87-92.

http://dx.doi.org/10.1016/j.cossms.2011.01.001

[4] Furrer, D.U. (2011) Application of Phase-Field Modeling to Industrial Materials and Manufacturing Processes. Current Opinion in Solid State and Materials Science, 15, 134-140.

http://dx.doi.org/10.1016/j.cossms.2011.03.001

[5] Zhou, S. and Wang, M.Y. (2007) Multimaterial Structural Topology Optimization with a Generalized Cahn-Hilliard Model of Multiphase Transition. Structural and Multidisciplinary Optimization, 33, 89-111.

http://dx.doi.org/10.1007/s00158-006-0035-9

[6] Takezawa, A., Nishiwaki, S. and Kitamura M. (2010) Shape and Topology Optimization Based on the Phase Field Method. Journal of Computational Physics, 229, 2697-2718.

http://dx.doi.org/10.1016/j.jcp.2009.12.017

[7] Uehara, T. (2014) Numerical Simulation of Foam Structure Formation and Destruction Process Using Phase-Field Model. Advanced Materials Research, 1042, 65-69.

http://dx.doi.org/10.4028/www.scientific.net/AMR.1042.65

[8] Uehara, T. and Suzuki, H. (2012) Numerical Simulation of Homogeneous Polycrystalline Grain Formation Using Multi-Phase-Field Model. Applied Mechanics and Materials, 197, 2610-2614.

http://dx.doi.org/10.4028/www.scientific.net/AMM.197.628

[9] Uehara, T. (2015) Phase-Field Modeling for the Three-Dimensional Space-Filling Structure of Metal Foam Materials. Open Journal of Modelling and Simulation, 3, 120-125.

http://dx.doi.org/10.4236/ojmsi.2015.33013

[10] White, D., Kimerling, A.J. and Overton, W.S. (1992) Cartographic and Geometric Components of a Global Sampling Design for Environmental Monitoring. Cartography and Geographic Information Systems, 19, 5-22.
http://dx.doi.org/10.1559/152304092783786636

[11] Gregory, M.J., Kimerling, A.J., White, D. and Sahr, K. (2008) A Comparison of Intercell Metrics on Discrete Global Grid Systems. Computers, Environment and Urban Systems, 32, 188-203.

http://dx.doi.org/10.1016/j.compenvurbsys.2007.11.003

[12] Mahdavi-Amiri, A. and Samavati, F. (2014) Atlas of Connectivity Maps. Computers & Graphics, 39, 1-11.
http://dx.doi.org/10.1016/j.cag.2013.09.003

[13] Du, Q., Gunzburger, M.D. and Ju, L. (2003) Voronoi-Based Finite Volume Methods, Optimal Voronoi Meshes, and PDEs on the Sphere. Computational Methods in Applied Mechanics and Engineering, 192, 3933-3957.
http://dx.doi.org/10.1016/S0045-7825(03)00394-3

[14] Jacobsen, D.W., Gunzburger, M. Ringler, T. Burkardt, J. and Peterson, J. (2013) Parallel Algorithms for Planar and Spherical Delaunay Construction with an Application to Centroidal Voronoi Tessellations. Geoscientific Model Development, 6, 1353-1365. http://dx.doi.org/10.5194/gmd-6-1353-2013

[15] Steinbach, I., Pezzolla, F., Nestler, B., Seeselberg, M., Prieler, R., Schmitz, G.J. and Rezende, J.L.L. (1996) A Phase Field Concept for Multiphase Systems. Physica D, 94, 135-147.

http://dx.doi.org/10.1016/0167-2789(95)00298-7

[16] Steinbach, I. and Pezzolla, F. (1999) A Generalized Field Method for Multiphase Transformations Using Interface Fields. Physica D, 134, 385-393. http://dx.doi.org/10.1016/S0167-2789(99)00129-3