Textural equilibrium melt geometries around tetrakaidecahedral grains.

In textural equilibrium, partially molten materials minimize the total surface energy bound up in grain boundaries and grain-melt interfaces. Here, numerical calculations of such textural equilibrium geometries are presented for a space-filling tessellation of grains with a tetrakaidecahedral (truncated octahedral) unit cell. Two parameters determine the nature of the geometries: the porosity and the dihedral angle. A variety of distinct melt topologies occur for different combinations of these two parameters, and the boundaries between different topologies have been determined. For small dihedral angles, wetting of grain boundaries occurs once the porosity has exceeded 11%. An exhaustive account is given of the main properties of the geometries: their energy, pressure, mean curvature, contiguity and areas on cross sections and faces. Their effective permeabilities have been calculated, and demonstrate a transition between a quadratic variation with porosity at low porosities to a cubic variation at high porosities.

Introduction

The physical properties of partially molten materials depend crucially on the geometry of melt at the scale of individual grains. Properties like permeability or electrical conductivity can be radically different depending on whether melt forms a connected network or not. The aim of this contribution is to better understand the controls on the geometry of melt networks in order to ultimately better understand the physical properties of partially molten materials.

Surface energy plays a key role in determining melt network geometry. In the absence of external forcing, partially molten materials tend to a state of textural equilibrium in which the surface energies bound up in grain boundaries and grain–melt interfaces are minimized. While in many situations a state of textural equilibrium is not achieved (due to the action of additional mechanical and chemical processes), it provides an important reference model for understanding melt geometry.

It is straightforward to write down a mathematical statement of textural equilibrium, but solving the resulting equations is much less straightforward. Early work [1-4] provided analytical solutions in some simple special cases, or used fairly rough approximations to the geometry. Only in the late 1980s was the simplest three-dimensional problem—four grains meeting at a junction with tetrahedral symmetry—solved fully numerically by von Bargen & Waff [5], Cheadle [6] and Nye [7]. This simple model with tetrahedral symmetry provides important insights into when a melt network is expected to be connected, and also provides constraints on the expected permeability [5,6,8] and electrical conductivity [6,9] of such networks. However, the tetrahedral-symmetry junctions have an important drawback: there is no space-filling solid phase compatible with such junctions. It is thus difficult to use the results from these studies in models which need to describe processes occurring within individual solid grains (e.g. as needed in modelling diffusion creep).

This article provides an exhaustive account of textural equilibrium melt geometries around a particular choice of solid grains which do fill space. In the absence of melt, these chosen grains take the shape of plane-faced tetrakaidecahedrons (truncated octahedrons). The problem has a high degree of symmetry, and is only marginally more complex than the problem tackled by von Bargen & Waff [5], Cheadle [6] and Nye [7]. Indeed, this problem has already been tackled in a series of recent contributions by Ghanbarzadeh et al. [10-14], using a novel level-set method that allows for greater flexibility in grain shape than the approach taken here. However, there are some important differences between the results of that study and this study, which will be discussed in detail in §8. The main advantage of the approach taken here over that of Ghanbarzadeh et al. [10] is that this study exploits all the symmetries of the problem, which makes it easier to resolve the fine details of the melt geometries.

This article is organized as follows. First the model is described, along with a recap of some well-known mathematical results about textural equilibrium. Accounts are then given of the main properties of the melt topologies: energy, pressure, contiguity, and areas on contacts and cross sections. Permeability calculations are then discussed, and are followed by a more general discussion relating these calculations to previous work. Appendices provide more detail on the numerical methods, and provide some analytical solutions for special cases.

The model

The model geometry consists of an infinite tessellation of tetrakaidecahedral unit cells as depicted in figure 1. In the absence of melt, a single grain occupies the whole of a unit cell, and the grain takes the shape of a plane-faced tetrakaidecahedron (truncated octahedron). The faces of the tetrakaidecahedron form the grain boundaries (solid–solid contacts). When melt is present, it is assumed that grain boundaries continue to lie along the faces of the tetrakaidecahedral cell, although the area and shape of these contacts are allowed to vary. Textural equilibrium involves the minimization of surface energy E [1,4], where γ and γ are the surface energies per unit of area of grain–grain (solid–solid) contacts and grain–melt (solid–liquid) contacts, respectively. A and A are the corresponding areas of solid–solid and solid–liquid contacts per unit cell. In this work, it will be assumed that γ and γ are isotropic and constant. The factor of in (2.1) arises from the fact that the grain boundaries are on the faces of the unit cell, so per unit cell each solid–solid contact only counts for half in the total surface energy [4].

Figure 1.

Three-dimensional rendering of the tetrakaidecahedral unit cell. Grain–grain (solid–solid) contacts are shown in pink. The region where melt (liquid) is present is shown in yellow. Translucency has been used to allow a view through the grain.

Minimization of E is performed subject to constraints, which reflect assumptions about the geometry. Here, grain centres are assumed to reside on a body-centred cubic (bcc) lattice, all grain boundaries are assumed to be planar, and all grains are identical. The unit cell shown in figure 1 is the Wigner–Seitz cell of the bcc lattice of grain centres. As a consequence of the cubic symmetry, the basic computational domain can be reduced to a single region that is 1/48th of a grain or 1/8th of a quadruple junction (figure 2).

Figure 2.

The fundamental computational domain for a porosity ϕ=0.03 and dihedral angle θ=30°. (a) Melt quadruple junction showing the solid–liquid interface with the fundamental computational domain in blue with finite-element triangulation in black. The rest of the melt junction can be obtained from the fundamental domain by applying the eight elements of the point group D2 (tennis-ball symmetry). (b) The same fundamental computational domain in blue, but shown in pink are the solid–solid contacts (grain boundaries) associated with the fundamental domain. The full grain can be produced by application of the 48 elements of the point group O.

Grain centres are fixed in these calculations, and thus there is a constant distance d between opposing square faces in the unit tetrakaidecahedral cell. The volume of the unit cell is V =d3/2 and it has area . Each edge has length . The volume fraction of melt is prescribed, by introducing a Lagrange multiplier λ and minimizing the functional J, where V is the volume of liquid in the domain and V * is the target (prescribed) value.

The variational problem is discretized by representing the surface by a finite-element mesh of triangles, using the Surface Evolver software [15-17]. Numerical optimization is used to find the unknown mesh node coordinates that minimize the surface energy subject to the given constraints. Further details on the numerical methods can be found in appendix A.

At a minimum, δJ=0, and from variational calculus it follows that, at the minimum, λ=λ*, where Here, H is the mean curvature of the solid–liquid surface, defined by where R1 and R2 are the principal radii of curvatures, and n is the normal vector (chosen here to point outwards from solid/inwards to liquid); ∇⋅ represents the surface divergence operator. Thus, in textural equilibrium, the mean curvature of the solid–liquid surface is a constant. Moreover, the Lagrange multiplier enforcing the volume constraint has a physical interpretation: it represents the pressure difference between solid and liquid, λ*=ΔP=P−P, and (2.3) is the Young–Laplace equation relating pressure differences to mean curvature.

Also from δJ=0 it follows that the following force balance holds along triple lines where three surfaces meet [18]: where is the co-normal to each surface at the triple line and γ is the corresponding surface energy per unit area. This reduces to the familiar expression for the dihedral angle θ at which two solid–liquid surfaces meet a solid–solid surface (figure 3),

Figure 3.

(a) The dihedral angle θ as the angle at which two solid–liquid surfaces meet a solid–solid surface. (b) The corresponding force balance at the triple line.

It should be noted that the approach taken here, of discretizing the variational problem (minimize J in (2.2)), differs from the approach taken in the classic studies by von Bargen & Waff [5], Cheadle [6] and Nye [7], and also the more recent work by Ghanbarzadeh et al. [10]. The starting point for all of these studies is the statement that the solid–liquid surfaces are of constant mean curvature and meet the solid–solid surfaces at the dihedral angle. These studies are all based on a numerical discretization of the mean curvature of the solid–liquid surfaces, whereas here in discretizing J only areas and volumes are discretized. The approach taken here, working from the variational problem, is essentially identical to the approach pioneered by Beeré [1], although here a more refined discretization is used.

By scaling, the variational problem can be reduced to being a function of just two dimensionless parameters: the porosity ϕ (the volume fraction of melt) and the dihedral angle θ. For example, if E represents the surface energy contained in one tetrakaidecahedral cell, then a scaled energy can be defined by E′=E/(γA), where A is the area of the bounding tetrakaidecahedral cell (). All results in this article are presented using scaled variables as functions of ϕ and θ.

Melt topologies

One of the key features of textural equilibrium is that different kinds of melt topology are possible for different values of the dihedral angle and at different porosities [2,3]. For example, there is the well-known result that, at small porosities, the melt network is only connected along the grain edges for dihedral angles less than 60°. The problem tackled here allows for a rich variety of melt topologies as a function of porosity and dihedral angle, which are summarized in the regime diagram of figure 4 and discussed below. For some combinations of parameters it is possible to find multiple topologies: each is a local minimum of the energy functional, but not necessarily a global minimum.

Figure 4.

Regime diagram showing regions where different melt topologies have been found to exist. ‘c’ denotes melt connected along grain edges (figure 5). ‘i’ denotes melt in isolated pockets at the grain corners (figure 6). ‘s’ denotes melt that wets the square faces of the grains but not the hexagonal faces (figure 7). ‘h’ denotes melt that wets the hexagonal faces but not the square faces (figure 8). ‘d’ denotes disaggregation, where grains are isolated spheres surrounded by melt. Pluses indicate regions where multiple topologies are found.

Figure 5.

Examples of melt topology ‘c’, connected along grain edges. On the left, in yellow,is a view of the quadruple junction. In the middle, in pink, is a view of an individual grain. On the right is a view of an individual grain + melt. Corresponding labels give the porosity ϕ and dihedral angle θ. The first two rows have positive mean curvature; the final row has negative mean curvature.

Figure 6.

Examples of melt topology ‘i’, where melt is isolated at the grain corners. Note that θ=65° has positive mean curvature, whereas θ=120° has negative mean curvature. Note that the bottom row has the same combination of porosity and dihedral angle as the bottom row of figure 5.

Figure 7.

Examples of melt topology ‘s’, where melt is connected and wets the square faces. Note that the bottom two rows have the same combinations of porosity and dihedral angle as the bottom two rows of figure 5.

Figure 8.

An example of melt topology ‘h’, where melt is connected and wets the hexagonal faces.

The first type of topology, marked ‘c’ in figure 4, has melt forming a connected network along the grain edges. Examples of such a topology are shown in figure 5, and are broadly similar to those depicted by von Bargen & Waff [5], Cheadle [6] and Nye [7], except that the melt junction does not have tetrahedral symmetry. The second type of topology, ‘i’ (figure 6), which consists of melt isolated at grain corners, is also broadly similar to previous calculations which assumed tetrahedral symmetry, for which analytical expressions are available [2]. At very large porosities the grain–grain contacts disappear and grains are completely surrounded by melt. This is marked as ‘d’ for disaggregated in figure 4, where the minimum energy configuration simply consists of spherical grains. ‘d’ topologies exist for all porosities greater than .

The problem considered here admits additional melt topologies because of the lower degree of symmetry of the quadruple junction. These are depicted in figures 7 and 8. The first of these represents the case ‘s’, where the square faces become wetted but the hexagonal faces do not. The second represents the case ‘h’, where the hexagonal faces are wetted but the square faces are not.

For low porosities, the dihedral angles for which topologies ‘c’ and ‘i’ are found are essentially the same as found for the problems with tetrahedral symmetry: isolated solutions exist only for dihedral angles greater than 60°. For dihedral angles greater than this, there is a region of overlap between the ‘c’ and ‘i’ topologies as porosity increases (the region between the ‘pinch-off’ and ‘wetting’ boundaries as described by von Bargen & Waff [5]).

Perhaps the most important new transition in the present study is the transition from ‘c’ to ‘s’ topologies as porosity increases (i.e. the preferential wetting of the smaller square faces). For low dihedral angles, this occurs close to ϕ=0.11, but for larger dihedral angles there is an overlap where both ‘c’ and ‘s’ topologies exist. The wetting of square faces at around ϕ=0.11 is also seen in calculations of wet Kelvin foams [19-24]. Wet Kelvin foams are essentially a limiting case of the problem considered here: the case where the dihedral angle approaches 0°. The network of melt considered here is termed a Plateau border network in the foam literature. The only difference is that a wet Kelvin foam allows for some curvature along what are the grain–grain contacts in the present model, but such curvature is extremely slight, and unlikely to significantly affect where the transition to wetted square faces occurs.

One curious feature of figure 4 is that for dihedral angles just below 100° the ‘s’ solutions exist down to very small porosities. If such a solution were realizable, there would be the potential for percolation of melt at small porosities even for some dihedral angles greater than 60°. However, as will be seen in the next section, such ‘s’ topologies have higher energy than the ‘i’ topologies and are thus less likely to be realized.

It should be noted that figure 4 does not provide a map of all possible melt topologies, only a key subset that have been chosen to be investigated. Firstly, the topologies investigated are only those consistent with the symmetries imposed. Topologies with lower degrees of symmetry are possible, e.g. isolated topologies with melt at some quadruple junctions but not others. Moreover, there are additional isolated topologies possible that are consistent with the imposed symmetry but that have not been calculated here. For example, melt can also be isolated in the middle of the grain edges, or at the centres of the grain faces. One can also have some combination of isolated melt in all three places—edges, faces and vertices of the unit cell. However, at the onset of melting it is expected that melt first forms at the quadruple junctions, so it is the case of melt at these junctions that is typically of most interest.

Energy and pressure

Figure 9 plots a scaled surface energy for each of the melt topologies as a function of porosity and dihedral angle. For those regions of ϕ−θ space that admit multiple topologies, figure 9 can be used to identify which topology has the overall lowest energy (and hence is more likely to be found). For example, the ‘h’ topologies (which wet the hexagonal faces) always have higher energy than the ‘s’ topologies (which wet the square faces). This is what one might intuitively expect—it is more likely for smaller faces to become wetted than larger faces.

Figure 9.

Scaled (dimensionless) energy, E′=E/(γA), plotted as a function of porosity ϕ for different values of the dihedral angle θ (key on the right). In this plot, and subsequent plots, thin solid lines represent connected topology ‘c’, thick solid lines represent square-wetted topology ‘s’, dotted lines represent isolated topology ‘i’ and dashed lines represent hex-wetted topology ‘h’. The black solid line represents disaggregated topology ‘d’ (here independent of the dihedral angle).

The situation is more complicated in the region of overlap between ‘c’ and ‘s’. For small dihedral angles, the ‘c’ topologies have lower energy than the ‘s’ topologies for most of the region of overlap, until close to the boundary where only ‘s’ exists. For larger dihedral angles the ‘s’ topologies are lower energy for most of the range of overlap. Similarly, in the region of overlap between ‘s’ and ‘i’, for lower dihedral angles the ‘i’ topologies are lower energy for much of the region of overlap, whereas at higher dihedral angles the ‘s’ topologies are lower energy for the whole region of overlap.

Figure 10 plots a scaled version of the pressure difference ΔP between solid and liquid for the various melt topologies. Owing to the Laplace–Young equation (2.3), this is also a plot of a scaled mean curvature. This pressure difference is calculated during the solution of the variational problem as ΔP=λ* is the Lagrange multiplier that enforces the fixed volume constraint. The curves in figure 10 are related to the slopes of the curves for energy in figure 9, because of the relationship which is a consequence of the variational problem. Here, E* represents the energy at equilibrium, and V * and V * the volumes of the liquid and unit cell, respectively. In the scaled (dimensionless) variables used in figures 9 and 10, this relationship becomes

Figure 10.

Scaled pressure (or scaled mean curvature), ΔP′=ΔPd/γ=2Hd, as a function of porosity.

Figure 10 shows a number of features that are consistent with previous studies. For example, for the isolated topologies ‘i’, the curvature is positive for θ<71° and negative for θ>71° [2]. Interestingly, topologies ‘i’, ‘c’ and ‘s’ have zero mean curvature for some combinations of porosity and dihedral angle, and are thus examples of ‘minimal surfaces’ (surfaces defined as having zero mean curvature).

The pressure difference between the two phases becomes singular as the porosity approaches zero, but the form of singularity varies with the dihedral angle. For isolated topologies the porosity dependence is of the form ΔP∝ϕ−1/3 (see (B.5) in appendix B). For small dihedral angles, the melt geometry can be closely approximated by tubes along the grain edges, for which ΔP∝ϕ−1/2 (see (B.16) in appendix B).

Effective pressure

The pressure plotted in figure 10 is a measure of the change in energy of the system with liquid content, assuming the cell size remains fixed (and thus as the volume of liquid increases, the volume of solid must decrease). As pointed out by Park & Yoon [4], in several physical situations what one would like to know is the change in energy with liquid content, but keeping the solid volume (V *) fixed. For example, in problems of compaction or sintering one may have a partially molten medium surrounded by a reservoir of liquid, and want to know whether liquid will be drawn in or expelled from the partially molten medium. Park & Yoon [4] quantify this in terms of an ‘effective pressure’, P, defined by Figures 11 and 12 show scaled energy and scaled effective pressure, analagous to figs 2 and 4 of Park & Yoon [4]. They are scaled such that, as porosity varies, the solid volume remains fixed (rather than the cell size). Scaling is made with respect to the length d0, which is defined as the distance between grain centres in the absence of melt with the same volume of solid grain. The area of the unit cell in the absence of melt is given as , and the volume as . The plot of effective pressure in figure 12 is related to the slope of figure 11 by

Figure 11.

Scaled energy, scaled in the same way as figure 2 of Park & Yoon [4], . Filled circles show the minimum values, which correspond to zero crossings in figure 12. Dash-dotted lines show the energy for the disaggregated topology ‘d’ (which varies with dihedral angle in this scaling, unlike the scaling in figure 9).

Figure 12.

Scaled effective pressure, scaled in the same way as figure 4 of Park & Yoon [4], .

If the sign of P is positive, it indicates that there is a driving force sucking liquid into the partially molten material. If negative, the driving force acts to expel liquid. Zero P represents a stable state (a minimum of the scaled energy plotted in figure 11). The key result of figures 11 and 12 is that topologies for dihedral angles greater than 60° are always unstable, and the driving force is such that the two phases try to separate (P is always negative). For dihedral angles less than 60° there exists a critical porosity at which the effective pressure is zero, and hence the topology is stable. The smaller the dihedral angle, the larger the critical porosity. If the partially molten medium has a porosity less than this critical porosity there will be a tendency of melt to be drawn in, whereas if it is above this critical porosity there will be a tendency for melt to be expelled.

In Park & Yoon’s original study, solutions with zero effective pressure could be found for dihedral angles up to 75°, rather than the 60° limit found here. One difference between Park & Yoon’s study and the present study is that Park & Yoon [4] consider a unit cell taking the shape of a rhombic dodecahedron. However, the discrepancy with the results of this study is likely to have arisen from Park & Yoon’s approximation of the grain–melt surfaces by circular arcs, which is not consistent with the surfaces being of constant mean curvature.

Geometrical properties

Figures 13–18 summarize a number of key geometrical properties. The first of these (figure 13) shows contiguity φ, which measures the relative area of grain–grain contact to total grain area. Figures 14 and 15 show the individual areas of solid–solid and solid–liquid contact, normalized to the area of the unit cell. Contiguity is used as a fundamental variable in a number of micromechanical models, e.g. in determining elastic properties [25] and effective viscosities due to creep [26]. For small porosities, figure 13 shows the same main trends as observed by von Bargen & Waff [5] and Cheadle [6], with contiguity being larger for larger dihedral angles. Contiguity shows a steady decrease with increasing porosity as more of the grain becomes wetted. There is a notable change in slope during the transition from ‘c’ to ‘s’ topologies, with a shallower slope for the ‘s’ topologies than the ‘c’ topologies.

Figure 13.

Contiguity. φ=A/(A+A).

Figure 18.

Scaled cross-sectional channel area (sectioned at the mid-grain edge), scaled as A/(d2ϕ). Isolated topologies ‘i’ have no melt at mid-grain edges and thus have zero area and are not shown. The black line labelled ‘u’ on the y-axis shows the expected channel area for melt tubes of uniform cross section along the grain edges.

Figure 14.

Total area of solid–solid contact, normalized to the area of the unit cell, A/A.

Figure 15.

Total area of solid–liquid contact, normalized to the area of the unit cell, A/A.

In this problem, there are two different kinds of grain–grain contact—those associated with the square faces and those associated with the hexagonal faces. Figure 16 plots the areas of the square faces, and figure 17 that of the hexagonal faces. For the ‘c’ topologies, the areas of the square faces decrease as porosity increases, but note that they do not vanish at the boundary between the ‘c+s’ and ‘s’ topologies: thus wetting of the square faces occurs discontinuously from a finite area of square grain–grain contact to a zero area as porosity increases. Only for a zero dihedral angle does the ‘c’ topology smoothly go to zero area of the square face as porosity increases. Such discontinuous jumps in behaviour are common in area minimization problems. The classic example of this is in the minimal surface of revolution problem: finding the soap film which minimizes the area between two parallel circular hoops [27]. Beyond a critical separation, the catenoid solution which joins the two hoops breaks down to the Goldschmidt discontinuous solution with two separate films spanning the two circles. For small dihedral angles, there is a close agreement between these results and those for wet Kelvin foams: hexagonal and square areas reported in figure 5 of [22] are very similar to what one would expect on extrapolating the results of figures 16 and 17 to zero dihedral angle.

Figure 16.

Square solid–solid contact area, normalized to the area of the unit cell, A/A. The ‘s’ and ‘d’ topologies have zero square contact area and are not shown.

Figure 17.

Hexagonal solid–solid contact area, normalized to the area of the unit cell, A/A. The ‘h’ and ‘d’ topologies have zero square contact area and are not shown.

Figure 18 plots a scaled version of the channel area A along the middle of each channel edge (i.e. the area of the truncated faces of the quadruple junctions shown in figures 5–8). If the melt resided in tubes along the grain edges of uniform cross section, one would expect for small porosities (black line on the y-axis in figure 18). For small porosities, and small dihedral angles, figure 18 shows that this value is approached, reflecting the fact that for such porosities and dihedral angles the melt geometry is well approximated by tubes. At small porosities, the mid-edge channel area decreases with increasing dihedral angle, reflecting the fact that more melt resides near the grain vertices than resides near the grain edges. The opposite trend is seen for the ‘s’ topologies at large porosity, where the larger dihedral angles have larger channel areas.

Permeability

The effective permeabilities of the melt topologies are shown in figure 19, and details of their calculation can be found in appendix C. Permeabilities are scaled in figure 19 by d2ϕ2, based on the expected behaviour for low porosities (see below). Most of the main features of figure 19 can be understood in the context of simpler models. For small porosities, and small dihedral angles, the melt geometry for topology ‘c’ takes the form of tubes along the grain edges. A very simple model of permeability was derived by Frank [29], consisting of tubes of melt with circular cross section lying on the edges of a tetrakaidecahedron. For this model, the permeability is As discussed by von Bargen & Waff [5], this formula gives an overestimate of the permeability at low porosity for the melt networks given here, because the melt tubes do not have circular cross section. For a cross-sectional shape more appropriate for low dihedral angles (three circular arcs meeting at zero dihedral angle), the melt flux for a given cross-sectional area and pressure gradient is a factor of 0.505 lower. A modified version of Frank’s formula is then and this formula is marked as ‘u’ on the y-axis of figure 19. The equation above explains the limiting behaviour seen in figure 19 for low porosities and low dihedral angles, which tend to finite values of k/(ϕ2d2) as , and whose geometries can be well described by tubes along the grain edges with approximately uniform cross section. For larger dihedral angles the cross sections become much less uniform, with more melt residing at the vertices of grains rather than the edges. The effect of this is to reduce the permeability. von Bargen & Waff [5] provide an approximate formula, which is in good agreement with the results here for small porosities and dihedral angles around 40°. Such a formula gives the permeability to within 15% for all porosities less than 2%. The expression quoted by Cheadle [6] (k=ϕ2d2/3000) seems to underestimate the permeability, and the permeabilities calculated here are much closer to the results of von Bargen & Waff [5].

Figure 19.

(a) Scaled permeability k/(d2ϕ2). The solid black line in the plot is for topology ‘d’, a bcc array of isolated spheres (the permeabilities calculated in this case are in agreement with previous calculations, e.g. ch 9 of [28]). The short black line labelled ‘u’ on the y-axis shows the expected permeability for melt tubes of uniform cross section along the grain edges. (b) Same data as shown in panel a, but zoomed-in to show the behaviour of ‘c’ topologies at low porosities (similar to figure 14 of von Bargen & Waff [5]).

For larger porosities (greater than 10%), lines on figure 19 are roughly linear, which reflects the fact that at large porosities permeability scales as k∝d2ϕ3. The ϕ3 behaviour can also be understood from simpler models, and is the expected behaviour for melt lying as thin sheets along the grain boundaries. A reasonable approximation for 0.1<ϕ<0.3 and the ‘s’ topologies is which lies within about 25% of the calculated values for dihedral angles θ<60°. Indeed, it should be noted that the relative range in permeability is very small when θ varies for porosities above 10%: a factor of 2 at most. It is interesting to note that the behaviour of permeability with dihedral angle differs at higher porosities: for high porosity, the larger the dihedral angle, the greater is the permeability; although the magnitude of this effect is rather modest.

Permeability data are often parameterized with a single law of the form k=ϕd2/C (e.g. [30,31]). Note that such a parametrization in terms of a single C and n value would not be able to describe all the subtleties seen in figure 19, and particularly the shift from k∝d2ϕ2 at low porosities to k∝d2ϕ3 at larger porosities.

Discussion

As remarked in the Introduction, the textural equilibrium geometries calculated here have also been calculated by Ghanbarzadeh et al. [10-14], although these authors consider a more limited range of porosities and dihedral angles. Ghanbarzadeh et al. [11] introduce a novel approach to calculating textural equilibrium geometries based on the level-set method. In this approach, the interface between the solid and liquid is represented by the level set of a function, and the interface evolves at a velocity proportional to the surface Laplacian of mean curvature. Over time, the surface should approach a state of constant mean curvature. A partial differential equation determines the evolution of the level-set function, and additional terms are added to the PDE to enforce the dihedral angle constraint. Ghanbarzadeh et al. [11] discretize the PDE using high-order finite differences.

One of the advantages of Ghanbarzadeh et al.’s approach is that it allows a great deal of geometric flexibility, including the ability to solve for many grains at once with different grain shapes. As such, they do not explicitly impose symmetry on the solution (as is done here). Their method is considerably more computationally expensive than the approach taken here (which can solve for the melt geometries in seconds).

In broad terms, many of the topologies calculated by Ghanbarzadeh et al. [11] are similar to those calculated here, but there are some notable discrepancies. One of the most important discrepancies concerns the wetting of the square faces. Ghanbarzadeh et al. [10] state that ‘In contrast with prevailing assumptions, the smaller square grain boundaries become fully wetted at θ=10° and ϕ≥5%.’ As can be seen in figure 4, this study places the transition to wetted square faces at closer to ϕ=11% for θ=10°: a factor of 2 difference. As remarked earlier, transitioning at ϕ=11% is consistent with results for wet Kelvin foams [19,22,23].

Another discrepancy is in the calculated mean curvatures of the melt topologies. Hesse et al. [32] presented a plot of mean curvature against porosity for isolated ‘i’ and connected ‘c’ topologies for a dihedral angle θ=90°. In theory, this plot should be similar to the curvature information shown in figure 10, but in fact bears little resemblance. First, Hesse et al. [32] showed ‘c’ topologies for porosities as low as 1.5%, whereas the suggested lower bound for ‘c’ topologies here is around 4% (very similar to the pinch-off boundary in fig. 7 of von Bargen & Waff [5]). Moreover, the behaviour of the mean curvature with porosity for the isolated ‘i’ topologies shown by Hesse et al. [32] is very different from that seen here. The shape of the isolated topologies is independent of porosity, apart from a scaling of the coordinates [2]. Thus mean curvature H should be proportional to ϕ−1/3 for the ‘i’ topologies, and be singular at zero porosity (as indeed seen in figure 10). The curve shown by Hesse et al. [32] does not show this behaviour, and instead appears to curve the wrong way, and seems to approach a finite value for small porosity.

These discrepancies matter, because they lead to different predictions about upscaled physical properties like permeability. Indeed, the permeabilities that have been reported by Ghanbarzadeh et al. [13,14] differ from those calculated here by a significant margin, in some cases by almost an order of magnitude. For example, for a dihedral angle of 10°, Ghanbarzadeh et al. [13,14] report that their permeability data can be fitted well by an expression of the form k=d2ϕ2.6/595.66. For a porosity of ϕ=0.01 this implies k/(d2ϕ2)≈1×10−4, whereas here a value of k/(d2ϕ2)≈7×10−4 has been estimated (figure 19). Thus, Ghanbarzadeh et al. [13,14] seem to underestimate the permeability by a factor of 7 for these parameter values. For a porosity of ϕ=0.1, Ghanbarzadeh et al.’s expression implies k/(d2ϕ2)≈4×10−4, whereas here k/(d2ϕ2)≈1.4×10−3 has been estimated, more than a factor of 3 greater. Perhaps the most probable reason for the discrepancies is that the simulations of Ghanbarzadeh et al. [11] are under-resolved, and struggle to accurately capture small-scale variations in the geometry, e.g. with cusps at a small dihedral angle or small radii of curvature at low porosities. The accuracy of the solutions obtained in this study is discussed in appendix A.

Another point of difference between this study and that of Ghanbarzadeh et al. is that Ghanbarzadeh et al.’s study essentially aims to produce surfaces of constant mean curvature compatible with the dihedral angle constraints. Constant mean curvature is a necessary, but not a sufficient, condition for minimum energy. Constant mean curvature guarantees that the energy is extremized, but does not say whether the surface is a minimum, maximum or a saddle point of the energy. Indeed, it is well known that some solutions of the Euler–Lagrange equations for area minimization problems can be unstable, as in the minimal surface area of revolution problem mentioned earlier. Hence it is possible that the scheme proposed by Ghanbarzadeh et al. produces melt topologies that are not minimum energy, and hence would not be found by the approach taken here.

The existence of multiple melt topologies has important consequences for upscaled physical properties (e.g. permeability, electrical conductivity, effective viscosities). First, transitions in such properties can be discontinuous on varying parameters, as the underlying transitions in topology are discontinuous. Moreover, as remarked by von Bargen & Waff [5], and more recently by Hesse and co-workers [14,32], there is the potential for hysteresis in topology and hence hysteresis in upscaled physical properties. For example, suppose the dihedral angle is 40° and porosity is increased slowly from zero such that textural equilibrium is maintained (figure 4). The topology would remain as type ‘c’ until a porosity of 12% when a discontinuous transition to the type ‘s’ topology would occur (wetting of the square faces), when the ‘c’ state is no longer stable. If porosity were then slowly reduced, the topology would stay as type ‘s’ until a porosity of 9% when a discontinuous transition back to ‘c’ would occur, as the grain–melt surfaces begin to touch on the square faces (de-wetting of the square faces). Similar hysteresis can occur in other parts of the regime diagram, for example the ‘i’ to ‘c’ transition discussed by von Bargen & Waff [5].

There are several natural ways in which the present work can be extended. In the same way that this study presents a model with a lower degree of symmetry than the studies by von Bargen & Waff [5], Cheadle [6] and Nye [7], it would be fruitful to examine the consequences of relaxing more of the symmetry assumptions. Symmetry plays a crucial role in constraining the solutions produced here, and topologies which are stable under the symmetry constraints imposed here may not be stable if some of these symmetry constraints are relaxed.

In this study, a tetrakaidecahedral unit cell was chosen, but it would be straightforward to repeat the same calculations for other choices of unit cell, such as the rhombic dodecahedral unit cell used by Park & Yoon [4]. In a tessellation of tetrakaidecahedrons, at each vertex four edges meet, and each vertex is thus said to have a coordination number of 4. In a tessellation of rhombic dodecahedrons, half of the vertices have a coordination number of 4, and half have a coordination number of 6. X-ray tomographic imaging suggests that the predominant coordination number in olivine–basalt systems is 4 [33], so the tetrakaidecahedral unit cell is the more appropriate geometry for the partially molten rocks of the Earth’s mantle.

To simplify the calculations, grain boundaries were forced to be planar, and forced to lie on certain symmetry planes. Grain boundaries are in general not perfectly planar, and these constraints could be relaxed in future work. Indeed it is the curvature of grain boundaries that allows for grain growth and coarsening, an effect that is suppressed in the present calculations by the assumptions that have been made. Another natural extension is to consider the anisotropy in surface energy, both between grains and between the grains and the melt. Anisotropy can lead to faceted solid–liquid and solid–solid interfaces (e.g. [34]) that will present challenges for numerical simulation. Another challenge for the future is to move beyond textural equilibrium, to study the interplay between deformation and melt network geometry.

Conclusion

The aim of this study has been to produce a simple reference model which describes the geometry of melt within a partially molten material. This model, based around a tessellation of tetrakaidecahedral unit cells, has an advantage over most previous models in that the melt geometry is compatible with a space-filling tessellation of grains. A lot of the fundamental behaviour of the melt topologies calculated here agrees with the classic models of von Bargen & Waff [5], Cheadle [6] and Nye [7], particularly at small porosities. The main differences occur at higher porosities, where topologies that wet the square faces of the tetrakaidecahedral grains exist.

These geometries form a starting point for the calculation of upscaled physical properties which depend on the geometry of melt at the grain scale. One example of such a property has been presented here, namely the permeability, but the intent is to describe a wider range of physical properties in future articles. In particular, calculations on the effective creep properties will be presented, building on the simpler grain-scale models of Takei & Holtzmann [26].