Applying diffusion coupled deformation theory, we investigate how the elastic properties of a solid body are modified when forced to keep its chemical potential aligned with that of its melt. The theory is implemented at the classical level of continuum mechanics, treating materials as simple continua defined by uniform constitutive relations. A phase boundary is a sharp dividing surface separating two continua in mechanical and chemical equilibrium. We closely follow the continuum theory of the swelling of elastomers (gels) but now applied to a simple two phase one-component system. The liquid is modeled by a local free energy density defining a chemical potential and hydrostatic pressure as usual. The model is extended to a solid by adding a non-linear shear elastic energy term with an effective modulus depending on density. Imposing chemomechanical equilibrium with the liquid reservoir reduces the bulk modulus of the solid to zero. The shear modulus remains finite. The stability of the hyper-compressible solid is investigated in a thought experiment. A mechanical load is applied to a rectangular bar under the constraint of fixed lateral dimensions. The linear elastic modulus for axial loading is evaluated and found to be larger than zero, implying that the bar, despite the zero bulk modulus, can support a weight placed on its upper surface. The weight is stabilized by the induced shear stress. The density dependence of the shear modulus is found to be a second order effect reducing the density of the stressed solid (chemostriction).
I. INTRODUCTION
The thermodynamic state of a simple homogeneous liquid in thermal equilibrium is determined by the temperature and one more thermodynamic field variable, which can be density, pressure, or chemical potential. The value of the total free energy can be found by simply scaling the corresponding values of a reference system in volume or mass. Separate specification of mass and volume is not necessary, the two quantities being essentially equivalent. This equivalence is also reflected in the density. An increase in mass can always be compensated by a proportional increase in volume, leaving the density invariant. This invariance is a key principle of the thermodynamics of uniform liquids, as formalized by Gibbs. (Most of the material relevant to this paper can be found in the work of Hansen and McDonald.1) Density is a natural state variable for the description of an isothermal mechanical process. Alternatively, pressure can also be used.
However, complications arise for solids. An isochoric deformation of a homogeneous solid body, preserving volume and density, can still change its energy. The reason is the shear rigidity of solids. This is the subject of the classical theory of elasticity, which predates the Gibbsian thermodynamics of liquids. The state variables in this theory are strain and temperature. A variation in density is accounted for by an isotropic stretch of a reference state. (Out of the large number of textbooks on classical elasticity, we have selected for citation the recent work by Lubarda and Lubarda.2) Density is not treated as an independent degree of freedom and certainly does not contain sufficient information for the specification of the state of a solid (indeed, strain can be inhomogeneous even if density is homogeneous). There is therefore a disparity in the description of simple uniform liquids and solids. This is not a problem at interfaces between materially closed systems. The boundary conditions of continuum mechanics can take care of this. The question is how to match these different descriptions at the interface between liquid and solid phases if the liquid and solid are in chemical equilibrium, such as is the case if the liquid is the melt of the solid. This is the question addressed in the present contribution.
The study of chemical equilibrium between continuum models of liquid phases is an established field with a long history. This was the main driving force behind the development of the density functional theory (DFT) of inhomogeneous liquids.1,3 The discussion in the present paper is, however, strictly limited to boundaries between uniform phases as described in the framework of Gibbsian thermodynamics. Gibbs extended his scheme by representing the interface as a sharp dividing surface with its own energy scaling with the surface area of the interface.1,4 However, while these surface properties cannot be missing from a realistic description of phase boundaries, interface energy and surface tension will not be considered here. The focus is on how to generalize the dividing surface scheme for an interface between liquid and solid, bridging the difference in thermodynamic continuum theory for the two aggregation states.
This first of all requires an expression for the chemical potential of the solid accounting for the dependence on shear strain. The treatment of the balance of forces at the interface will also have to be refined now, involving the Cauchy stress tensor rather than only the hydrostatic pressure. Equal chemical potential no longer implies equal pressure (even in the absence of surface tension). Instead, chemical and mechanical equilibrium must be treated as separate but coupled boundary conditions. The theory to accomplish this generalization already exists. It is the diffusion coupled deformation (DCD) theory originally developed for the description of diffusion induced stress in solid solutions. The pioneers in this field are Larché and Cahn.5 Another influential theoretical contribution is the 1996 paper by Gurtin.6 For an update on the experimental evidence for Larché–Cahn theory, see Ref. 7. The theory was further extended for the study of the swelling of polymeric gels.8–14 A more recent application is in lithium ion battery research, in particular, for the understanding of the mechanics of charging and discharging of silicon cathodes.15–18
The investigation in this paper follows the diffusion coupled elasticity theory as presented in the book by Gurtin, Fried, and Anand (GFA)19 almost to the letter (including the notation). There is, however, one difference. The diffusive motion of a metal species in alloys is defined relative to a lattice of sites, which can be occupied by atoms of different species or even vacancies. Similarly, in a gel, the network of elastic polymers in the dry state serves as the reference to quantify the migratory displacement of the solvent. However, the solid in our application is a simple one-component system. There is no external lattice providing a reference. This is not a fundamental objection. As pointed out by Baek and Srinivasa, an equivalent separation can be achieved by a decomposition of the mass current in a convective and diffusive flow.8 The continuum mechanical basis of the theory remains valid. Differentiation between a migratory component (solvent) and an elastic but immobile component (elastomer) is made by imposing certain “swelling” constraints on the balance laws.8,13,14 Here, we carry on without constraints proceeding to the next stage, the design of a constitutive model of the simple solid.
The approach taken in this paper is somewhat unusual from the perspective of a physical chemist. DCD theory is an application of the methods of continuum thermodynamics.19,20 The principle equations in continuum thermodynamics are balance equations involving time derivatives of densities, entropy, and fluxes. Dissipation is built in at a fundamental level (see also Ref. 21). It may seem that using non-equilibrium equations to derive thermal equilibrium values is an unnecessary detour. Variational methods based on free energy minimization are more direct. However, the dynamical route has certain advantages even under equilibrium conditions. System specific constitutive relations are kept separate from generally valid field equations. This guards against inconsistencies, which are more of a risk for “multiphysics” theories, which introduce further state variables in addition to strain.19,22 Mass density in DCD theory can be regarded as such an extended field variable. Moreover, the distinction between density and deformation, as proposed here, is based on a decomposition of the mass current, which is a dynamical quantity. Here, it is perhaps helpful to draw a parallel to the Lorentz continuum theory of polarization for which the foundation is also a dynamical theory, the Maxwell–Lorentz equations.20,23 In fact, as suggested by Baek and Srinivasa, mass density can be treated as the divergence of mass polarization, replacing density as the fundamental degree of freedom.
Mass exchange across a liquid–solid interface also requires flexibility in constitutive modeling, allowing for a natural transition between liquid and solid. The free energy density of simple liquids is due to the short range steric repulsion and long range cohesive attraction characteristic of atomic fluids. These interactions, while treated in a local approximation, are poorly described by linearization of isotropic expansion as used in elastic theory. This suggests to retain the liquid-like local free energy density in the solid, with a different parameterization, and add a shear energy term depending on isochoric deformation. For an isotropic solid, such a model contains just one parameter more determining the strength of the shear elastic energy term. The continuum solid in our model is missing a reference lattice. Setting the shear modulus to zero, therefore, directly turns the solid into a liquid. Unfortunately, this makes the model intrinsically non-linear, forcing us to use the demanding machinery of non-linear continuum mechanics.19,24 Non-linearity is also the price to pay if we want to allow for a density dependence of the shear elastic interaction. This leads, as we will show, to chemostriction, creating a pressure gap between solid and liquid. While small, this effect is of interest as a manifestation of nonlinear elasticity.
Finally, a word about the motivation for the work reported here: While chemomechanical equilibrium at liquid–solid interfaces is of physical interest, what prompted the author to look into diffusion coupled deformation was a theoretical problem encountered in a previous study on electromechanics of dielectric fluids.25 The aim of that investigation was a systematic derivation of the stress tensor of a simple dielectric fluid polarized by an external electric field and to compare the corresponding force density to the force density obtained from density functional treatment of the same free energy functional. The two expressions for this observable quantity, generally known as the Korteweg–Helmholtz force, were found to be identical. More specifically, it was shown that the density obtained by optimization of the deformation for a fixed but arbitrary material density in reference space satisfies the variational DFT equation for density in the current (deformed) space. It is not at all clear whether this is a general principle. The results of the present study suggest that at least for inhomogeneous systems material density and deformation are independent degrees of freedom, which must be computed separately by solving a set of coupled equations. We will return to this issue in the conclusion.
With the fundamental questions mentioned above in mind, the present paper is a joint theoretical and applied investigation. The emphasis is on the thermomechanical methodology. While almost completely gathered from the textbook of GFA, which has become a standard reference in the continuum mechanics literature, this approach may not be familiar to everybody in the physical chemistry community. The derivations in the theoretical sections are therefore more explicit and lengthy than strictly necessary claiming the larger part of this paper. Only Sec. VII is devoted to the application of a supported bar of solid in a bath of its melt. This application is meant more as a thought experiment. The simple continuum model developed here is too primitive for the investigation of realistic solid–liquid phase boundaries. We conclude with a discussion in Sec. VIII, where we also comment on the connection to the DFT of fluids and solids.
II. KINEMATICS AND BALANCE EQUATIONS
A. Separating diffusive and convective currents
The diffusion coupled deformation methodology used in this work is patterned after the theory for the absorption (or release) of solvents by polymeric gels developed by the continuum mechanics community.8–14 (A recent review from a more physical chemistry perspective can be found in Ref. 21.) The solvent and the polymeric network are clearly distinct interacting components of a composite system, partly liquid and partly solid. The dry polymer is usually taken as the reference for both the elastic deformation of the polymer and the migration of the solvent. However, it was pointed out by Baek and Srinivasa that the balance equations can be formally developed, initially assuming that the gel is a one-component system with certain constraints imposed at a later stage8 (see also Refs. 13 and 14). We have adopted their approach leaving out the constraints. The argument presented in this section is a justification for treating the deformation of a one-component system and the density in the reference frame as separate independent degrees of freedom.
Separation of the total current j = ρv + h in convective mass flow ρv (solid arrows) and diffusive current h (dashed arrows), as written in Eq. (13). On the right is a deformed convecting volume (solid contour) with, on the left, its shape in material (reference) space (dashed contour). The curved arrows indicate the map χt of Eq. (5), moving points X from reference space to the positions x in the deformed body at time t. hR is the material diffusive current obtained by pulling the spatial current h back according to Eq. (24). While convective motion cannot cross the boundaries of the body, diffusive current h can. There is therefore residual diffusive current in reference space but no convective flow.
Separation of the total current j = ρv + h in convective mass flow ρv (solid arrows) and diffusive current h (dashed arrows), as written in Eq. (13). On the right is a deformed convecting volume (solid contour) with, on the left, its shape in material (reference) space (dashed contour). The curved arrows indicate the map χt of Eq. (5), moving points X from reference space to the positions x in the deformed body at time t. hR is the material diffusive current obtained by pulling the spatial current h back according to Eq. (24). While convective motion cannot cross the boundaries of the body, diffusive current h can. There is therefore residual diffusive current in reference space but no convective flow.
B. Mass balance in referential form
C. Rigidity and local Gibbs–Duhem relation: Overview
The decomposition [Eq. (13)] of the mass current j in hydrodynamic flow ρv and diffusive flux h is rather suspect from the perspective of Newtonian hydrodynamics. There is no mass diffusion in one-component fluids (h = 0). The only diffusion process is momentum diffusion, which is manifested as viscosity. This changes when the system is solidified. The mobility of atoms is now hindered. Rigidity sets in, and solids respond to the applied shear stress not by viscous flow but by elastic strain. However, the mobility of atoms in not fully arrested. The thermal concentration of vacancies still allows atoms to change places and migrate over finite distances, given enough time. Under these conditions, the distinction between convective motion (the deformation of the lattice) and particle diffusion is physical even in a one-component system. This was formalized in the seminal paper by Martin, Parodi, and Pershan (MPP) on crystal hydrodynamics.26 MPP argued that vacancy diffusion must be considered the third Goldstone mode complementing the two shear modes to the three modes required by the Goldstone theorem applied to the translational symmetry breaking in crystals. This implies that vacancy fluctuations are not only long time but also long range. This analysis was further elaborated in Ref. 27 and is the subject of an important chapter in the textbook by Chaikin and Lubensky.28
In addition to point defects, crystal mechanics must account for further defects not displayed by liquids. These are dislocations and other topological defects.28 Topological defects are also observed in mesophases such as liquid crystals and are characteristics of rigidity. Their dynamics can again differentiate between diffusive motion and convected motion. Restructuring of the near singular distortion inside dislocation cores is diffusive and is the rate limiting process for dislocation mobility. The relaxation of the long range elastic deformation induced by the dislocation is faster and can proceed by convective motion. This issue has been addressed in detail in recent studies of dislocation motion using the phase field crystal (pfc) method.29,30 The pfc method is a microscopic approach based on a modulation of the periodic structure of atomic density in crystals and is also used for the modeling of solidification and grain boundaries.31–35 These results suggest that separation between diffusion and convective motion is linked to the microstructure of density and could, therefore, also occur in systems without long range order such as fluids confined by hard walls.
Can the local Gibbs–Duhem equation [Eq. (28)] be given a similar interpretation in terms of some kind of scaling relation? To answer this question, we first need a local equivalent of the thermodynamic derivative relations [Eqs. (33) and (34)]. This is the purpose of the Lagrangian density factorization [Eq. (20)]. It enables us to treat the local particle number represented by ρR and volume represented by J as independent field variables. The flux decomposition [Eq. (13)] is needed to define the corresponding adjoint thermodynamic forces using the thermomechanical scheme for diffusion coupled deformation, as developed by GFA. This procedure is outlined in Sec. III. Application to the simple liquid in Sec. IV reproduces, as expected, Eq. (28), which in the process is elevated to a constitutive relation for a simple fluid. This result is valid whether the partition of current in a diffusive and convective component is physical or not.
Transverse stress is evaluated by solving the Cauchy equation taking account of the boundary conditions imposed by mechanical surface loading. In fact, body forces can be ignored for typical elasticity problems.2 The transverse stress is all that matters. It will also turn out the key point in the application presented in Sec. VII. Moreover, substituting for p in Eq. (29), the mechanical pressure p = −trT/3 is, in general, not correct either. The divergence of a (symmetric) tensor is not equal to the gradient of its trace. Ultimately, this is all again a consequence of rigidity. However, note that these complications not only apply to solids but also can affect inhomogeneous liquids. While liquids yield to shear stress, the equilibrium stress tensor at liquid interfaces can still have a deviatoric component. This has led to some controversy in the DFT of liquids.4 For a far reaching discussion, we refer to the papers by Percus and co-workers36,37 (see also the work of Lovett and Baus38,39).
The analysis in this paper remains strictly within the confines of continuum mechanics. Continuum thermodynamics is a self-contained theory, assuming a set of axioms (see the book of GFA). Derivations can, in principle, be carried out independently of an atomistic picture. A microscopic theory is needed at the application stage to fix constitutive relations for specific systems. The alternative is to obtain these relations from the experiment. Of course, a foundation in statistical mechanics is the ultimate goal in condensed matter science. This applies, in particular, to the separation in a convective and diffusive current for one-component systems introduced in Eq. (13). Fuchs and co-workers made significant progress toward this end for the case of crystal hydrodynamics.40,41 This concerns the problem raised by MPP macroscopic theory mentioned earlier.26,27 The question for microscopic theory is how to discriminate between the density in a crystal with vacancies and lattice deformation. Häring et al. designed a deterministic coarse graining procedure expressing these quantities in terms of instantaneous atom position. The advantage of defining hydrodynamic field variables by using a deterministic coarse graining scheme is that the fluctuations of these quantities can be analyzed by the methods of statistical mechanics and related to microscopic correlation functions.41 The coarse graining scheme of Ref. 40 is based on reciprocal space methods. An up-to-date review on real space deterministic coarse methods can be found in Ref. 42.
III. FREE ENERGY AND STATE EQUATIONS
A. Free energy imbalance
Continuum thermodynamics is a non-equilibrium theory in which fluxes play a central role. This is how we justified in Sec. II the inclusion of material density in the set of state variables describing the mechanical evolution of a single component system. The dynamical nature of continuum thermodynamics will again be used to derive equilibrium equations for these variables. Alternatively, the chemomechanical equilibrium conditions used for the application in Sec. VII could also be obtained using (free) energy-based variational methods. This is, in fact, the procedure followed by Baek and Srinivasa.8,13 However, the DCD theory of GFA is deeply rooted in non-equilibrium continuum thermodynamics, exploiting the constraints set by entropy production (see below). We will adhere to this philosophy and practice mainly because the constitutive theory for the solid proposed in Sec. V B is using another idea from the book of GFA, which is again formulated in terms of non-equilibrium continuum thermodynamics concepts.
A decisive step in the transformation of continuum thermodynamics into a practical tool for calculations was made by Coleman and Noll who showed how the second law can be used to derive formal state equations, given a set of state variables.43 This powerful procedure has become very popular in the continuum mechanics literature since it was proposed in 1963. However, departing in spirit and implementation from established methods in DFT, the Coleman–Noll (CN) procedure has received little attention from the physical chemistry community. Since the present paper is largely aimed at this audience, some of the key equations underlying the CN procedure are summarized below. Our main source is again GFA. A compact and accessible introduction can be found in the book by Kovetz20 who applied the CN procedure in a derivation of the Maxwell stress tensor.
The free energy imbalance [Eq. (38)], also known as the Clausius–Duhem equation, is the result of the combination of relations for part-wise energy balance and entropy production, assuming that the temperature is homogeneous and constant. Free energy imbalance plays an absolutely central role in the thermodynamics of deformable systems. The inequality is applied together with force balance equations such as the Cauchy equation [Eq. (40)], which are assumed, given, or are derived using the principle of virtual power.19,22 Here, the force balance equation is taken as a proven premise. A derivation using the virtual power balance can be found in Ref. 10.
B. Stress and chemical potential
Equations (62) and (63) are what one would expect from a direct variational calculation carried out in reference space using ψR as free energy density. Equation (64) goes beyond what free energy minimization can produce. It is a kinetic relation reminiscent of Fick’s law for diffusion. In the present investigation, only the thermodynamic derivative relations [Eqs. (62) and (63)] are needed. For further discussion about kinetics, the reader is referred to the book of GFA.
IV. SIMPLE CONTINUUM LIQUID
A. Constitutive relations for the simple liquid
Liquids have no memory. A Lagrangian detour via reference space should be an unnecessary exercise. Moreover, if we decide to take this route, we expect the choice of material density to be arbitrary. In this section, we verify that this is, indeed, the case. This derivation has been included not only for pedagogical reasons but also because this issue is going to be relevant later on for the constitutive definition of a simple solid at the melting line. We want the solid to resemble a liquid as much as possible in order to limit the number of order parameters to a minimum.
B. Recovering the local DFT of liquids
It may seem that this elaborate continuum thermodynamics argument takes a lot of effort to reestablish a relation that was essentially derived in a couple of lines in Sec. II C. Note, however, that the treatment of the chemical potential μ is more general. μ was defined in Eq. (41) as the energy transferred in diffusive transport. This eventually led via Eq. (55) to the local free energy imbalance [Eq. (57)]. Particle transport remains active even in inhomogeneous closed systems without outside supply of particles [h = 0 in Eq. (41)]. The continuum thermodynamics definition of chemical potential is, therefore, more general than the role assigned to it in DFT. We regard this as a significant advantage of the continuum thermodynamics view, ultimately justifying the definition of chemical potential for elastic solids in Sec. VI.
V. SIMPLE CONTINUUM SOLID
A. State equation for diffusion coupled deformation
B. Separating volumetric and isochoric deformation
VI. ISOTROPIC ELASTIC SOLID
A. Stress tensor and chemical potential
Shear stress is the signature of the breakdown of mobility at the microscopic level. Failing a detailed representation of atomic position in our model, the elastic coefficient Λ plays the role of a continuum order parameter for the transition from a liquid to a simple solid. This is why Λ has been made dependent on the spatial density ρ and not on the material density ρR in the reference frame. In this way, the shear elastic response can be triggered by the modest difference in density between the liquid and solid phases. The free energy density of Eqs. (108) and (109) can therefore be thought of as defining a Landau-type free energy function for solidification. However, in a continuum mechanics framework, this distinctly non-standard construction raises a number of technical questions. In particular, shear strain is now an isochoric deformation of a reference space attached to the current deformed density. Further discussion will be deferred until after the presentation of the application in Sec. VII.
Comparing to Eq. (74) for the liquid, we see that, in addition to a modification of the isotropic pressure term, now there is also an anisotropic term proportional to the shear Lamé modulus of the elastic model. This is, of course, a familiar result in non-linear elasticity theory (see the work of Ogden24). However, this expression is commonly given in terms of the full left Cauchy–Green tensor B. The fact that we can replace B by Bi is a convenient feature of the volumetric–isochoric decomposition of Eq. (92). The other important point to note is that the deviatoric stress of Eq. (123) and corresponding isotropic (total) pressure of Eq. (120) were derived in conjunction with a consistent expression for the chemical potential. Equation (113) makes it possible to define chemomechanical coupling coefficients in the form of the chemistry strain tensor introduced in Sec. VI B.
B. Chemomechanical response
C. Linear elastic modulus parameters
VII. PERMEABLE LIQUID SOLID INTERFACE
A. Mechanochemical equilibrium at the dividing surface
The liquid plays the role of a reservoir with fixed chemical potential and pressure . These quantities are not degrees of freedom but parameters. Equations (172) and (174) are therefore a set of coupled chemomechanical equations for the equilibrium reference density ρR and elastic degrees of freedom J and Ci of the solid.
B. Linearization: Solid bar immersed in its melt
Constrained axial loading of a solid bar (dark gray) immersed in a bath of its melt (light gray). Solid and liquid phases are in chemomechanical equilibrium (double half arrows). The bar is supported by a rigid substrate (black) inhibiting deformation in the two and three direction (stretches λ2, λ3 are unity). A weight is placed on the upper surface, applying a stress force pa compressing the bar in the one direction leading to a strain ϵ < 0. Experiments of this type are standard for the study of the swelling of polymeric gels. (A detailed discussion can be found in the textbook by Doi.44). Here, this setup is transferred in a thought experiment to a solid–liquid interface of a simple one-component system.
Constrained axial loading of a solid bar (dark gray) immersed in a bath of its melt (light gray). Solid and liquid phases are in chemomechanical equilibrium (double half arrows). The bar is supported by a rigid substrate (black) inhibiting deformation in the two and three direction (stretches λ2, λ3 are unity). A weight is placed on the upper surface, applying a stress force pa compressing the bar in the one direction leading to a strain ϵ < 0. Experiments of this type are standard for the study of the swelling of polymeric gels. (A detailed discussion can be found in the textbook by Doi.44). Here, this setup is transferred in a thought experiment to a solid–liquid interface of a simple one-component system.
Constrained axial contraction of a rectangular bar decomposed in volumetric compression followed by isochoric shear. The stretches λi, i = 1, 2, 3 for each operation are given below. Indicated are also the corresponding linear elastic moduli: M is the so-called P-wave modulus, K is the bulk modulus, and G is the shear modulus. In the linear approximation, M = K + 4G/3 [cf. Eq. (180)]. For an open system exchanging mass with a liquid reservoir, as realized in the system of Fig. 2, the bulk modulus vanishes (K = 0), eliminating the free energy cost of the volumetric compression (see Sec. VII B). The shear modulus G > 0 remains finite and therefore also M.
Constrained axial contraction of a rectangular bar decomposed in volumetric compression followed by isochoric shear. The stretches λi, i = 1, 2, 3 for each operation are given below. Indicated are also the corresponding linear elastic moduli: M is the so-called P-wave modulus, K is the bulk modulus, and G is the shear modulus. In the linear approximation, M = K + 4G/3 [cf. Eq. (180)]. For an open system exchanging mass with a liquid reservoir, as realized in the system of Fig. 2, the bulk modulus vanishes (K = 0), eliminating the free energy cost of the volumetric compression (see Sec. VII B). The shear modulus G > 0 remains finite and therefore also M.
C. Chemostriction and stability
Having determined the changes in geometry, we are ready for a closer investigation of the energetics and stability. First of all, there is the issue of the zero bulk modulus. We were led to this conclusion by the linear analysis of Sec. VII B. However, the lack of resistance to volumetric expansion is an inherent property of the original non-linear constitutive model, as formulated in Eqs. (108) and (109). With equal principle stretch in all directions, the isochoric Cauchy–Green tensor Bi is unity. The trace should be 3, as confirmed by evaluating Eq. (155) for λ1 = λ2 = λ3. Similarly, Eq. (157) for the deviatoric stress is the zero tensor. The shear elastic energy is eΛ = 0, and the Cauchy stress tensor becomes . There is no difference with the hydrostatic stress tensor of the liquid other than the actual density dependence of the LDA pressure .
This is what the volumetric–isochoric factorization of deformation built into the constitutive model of the solid was designed for. Under load and strain-free chemomechanical equilibrium, and . The liquid/solid dividing surface can be displaced, expanding or shrinking the volume of the solid, at no free energy cost, provided that the shape is not changed. This picture is consistent with classical nucleation theory. When there is no net decrease in free energy to be gained by enlarging the volume occupied by the growing phase, the nucleus will collapse under the stress exerted by the surface tension. Surface tension is not accounted for in our model, and hence, the volume of the solid body is arbitrary. Clearly, the introduction of surface tension is badly needed to make the diffusion coupled deformation scheme proposed here more realistic. Further discussion of this point is deferred to the concluding remarks in Sec. VIII.
VIII. SUMMARY AND DISCUSSION
This paper reports on a continuum mechanics investigation of the mechanical response of a solid in equilibrium with its liquid phase at the melting line. The liquid is modeled by the usual free energy density depending on local density only (no gradients). Steric repulsion and cohesion in the solid are represented by a similar local free energy density. Allowing for a different parameterization, this constitutive design makes the solid resemble the liquid as closely as possible, which is appropriate for a solid in equilibrium with its melt. The crucial difference is that the local density functional in the solid is extended with an elastic energy term accounting for shear deformation. The magnitude of the elastic term is controlled by a non-linear Lamé shear modulus varying with density. This elastic coefficient therefore effectively acts as an order parameter. A finite value turns a liquid into a solid without explicit modeling of positional ordering of atoms. Interface energy and surface tension have not been included in the present minimal model.
Chemo-mechanical equilibrium between solid and liquid is imposed by requiring that the chemical potential and normal component of the Cauchy stress tensor are continuous across the interface. The expressions for the chemical potential and stress tensor of the simple continuum solid are derived by applying diffusion coupled deformation theory. This theory was originally developed for the study of diffusion induced stress in multi-component systems, but following a suggestion by Baek and Srinivasa,8 the formalism is applied here to a one-component two phase system. Consistent with classical nucleation theory, we find that the bulk modulus of a solid body immersed in a bath of its melt vanishes. However, the solid keeps its shear rigidity. This is demonstrated by evaluating the elastic modulus for constrained axial loading of a rectangular bar. This example also illustrates how density dependence of the shear modulus gives rise to chemostriction. This somewhat unconventional name has been chosen because of the similarity to electrostriction in a dielectric fluid polarized by an external electric field. Chemostriction is manifested in a decrease in equilibrium density and pressure of the solid under mechanical loading. Similar to electrostriction, the effect is quadratic in the applied force and therefore the same for tensile dilation and compression.
The application is the subject of the final section with the larger part of this paper taken up by a detailed and inevitably technical explanation of the methodology. The reason that this was felt appropriate is that the diffusion coupled solid mechanics and constitutive theory as applied here have a number of non-standard aspects. There is first of all the question of the background against which diffusive motion is measured. In the theory of solvent permeation of gels, the reference is provided by the polymeric network. In Sec. II, we proposed that in a one-component system, the convective flow as defined in continuum mechanics can be used as a reference for migration of mass. Changing over to a Lagrangian picture, this argument leads us to treat material mass density in the reference space as an independent degree of freedom in addition to the deformation. The diffusion coupled deformation theory as presented by Gurtin, Fried, and Anand (GFA)19 is then applied without modification, leading to expressions for the chemical potential and stress tensor in terms of derivatives of the free energy density.
The constitutive model of the solid also makes non-standard assumptions about reference. This concerns the reference for shear deformation. The elastic coefficients in the model, including the non-linear shear modulus, are functions of the spatial (deformed) density. Similar to a liquid, specification of a material space is not needed. Current space is the (constitutive) material space. To implement this idea for the shear motion, the deformation gradient tensor was factorized in a volumetric expansion and an isochoric deformation following a mathematical scheme again taken over from the book of GFA. Shear deformation is effectively defined as displacive motion of covectors tangent to surfaces of constant volume density as quantified by J. However, isochoric deformation, while preserving local volume, can change inhomogeneous density. The example in the present application is Eq. (152), linking the change in density of the solid at constant chemical potential to the energy stored in the shear strain. As a result, shear deformation can modify the shear modulus, creating a self-consistency problem. This complication was evaded in the present application by linearization of the final mechanical equilibrium equations. The coupling to density was shown to be a second order effect and can be ignored for in the determination of the first order shear. However, in future fully non-linear applications (see below), these inconsistencies can no longer be ignored. Hopefully, we can learn from anelasticity theory where coupling of deformation to an evolving structure in material space is a major topic of research. The relatively “accessible” example is thermo-elasticity.19 Thermal expansion leads to another stress-free reference state.
The thermomechanical formalism and non-linear constitutive model developed in this paper are applied in a thought experiment, probing the mechanical response of a block of solids in chemomechanical equilibrium with its liquid phase. We first verify that enforcing equilibrium mass exchange across the dividing surface sets the bulk modulus of the solid body to zero. The application then continues showing how the mechanical response is stabilized by the remaining shear rigidity. This is illustrated by evaluating the elastic modulus for constrained axial loading of a solid bar in the bath of its melt. This analysis is carried out by linearizing the constitutive elastic model of the solid. The resulting strain is then used to estimate the chemostriction due to loading. The chemostriction is a consequence of the density dependence of the shear modulus and is a quadratic effect.
Clearly, our elementary model system is far removed from a phase equilibrium of a realistic solid–liquid system, such as ice cubes in a beaker of water. First of all, what is missing dearly is surface energy and surface stress. Moreover, for the description of the shape of a crystal, the surface energy will also have to be orientation dependent. This level of modeling is required to address questions such as the generalization of the Young–Laplace equation to solid interfaces.50 Furthermore, the migration of realistic solid–liquid phase boundaries can be expected to be hindered by kinetic barriers, such as activated accretion. However, we believe that the key message of the thermomechanical exercise reported in this paper still stands, namely that chemical and mechanical equilibrium must be treated separately at interfaces.
A simple elastic model, extended with surface tension, could already be useful for the investigation of soft capillarity. Similar to externally applied surface loads, surface tension will equally deform the interior of a soft solid. This is the subject of the new and fascinating field of elasto-capillarity.45–47 Solid nucleation will be a more ambitious goal. Again, the question can be asked how the capillary stress is modified by chemical equilibrium at the dividing surface. What will be the effect on the nucleation barriers? Elastocapillarity is not part of the analysis in most microscopic numerical studies of solid nucleation48–51 and freezing of confined nanophases.52 It could be important because at the melting line, solids are soft. The study of nucleation will also force us to reconsider another aspect of sharp interface modeling applied in this paper. The boundary between solid and liquid, while flexible in the deformed system, is fixed in material space, defining the reference volume and shape of the solid body. This issue touches on the fundamental question about the relation between thermodynamic tension and elastic stress.53
We conclude with a very brief comment on the relation to the DFT of inhomogeneous fluids. Density gradients can generate deviatoric stress. The well-known example is the Korteweg capillary stress tensor.54 This tensor can be derived from the square gradient (Cahn–Hilliard) free energy density by a Lagrangian procedure similar to the one applied in this paper. The free energy density is pulled back to a reference space and differentiated with respect to the deformation gradient, yielding a Piola stress tensor. Pushed forward to current space, the result is identical to the Korteweg stress tensor.55 There is finite stress but no strain (F = 1) as is also consistent with the freedom of choice of reference space. This raises the question that can stress in an inhomogeneous fluid, in principle, generate spontaneous (“eigen”) strain, given a suitably accurate density functional? In other words, can a variational minimization of a pure DFT functional extended with a search for the optimal shear deformation tensor break symmetry finding F ≠ 1? If it does, is this point where the fluid is no longer a fluid but a solid? This possibility is suggested by the observation that the higher order gradient functionals used in phase field crystal theory31–35 exhibit an elastic response. The same is true for the non-local density functionals used in the DFT of crystals.32,56,57 Unfortunately, evaluation of the stress tensor of a pattern forming (Hohenberg–Swift) density functional33 is very demanding and probably not very practical, at least not without further coarse graining.29 Here, a hybrid strain–density functional could be an option. The strain could play a similar role to the one-electron orbitals in the Kohn–Sham approach to the DFT of the electronic structure.
Similar questions concerning spontaneous strain in formally liquid systems arise in electric double layers modeled by Poisson–Boltzmann theory. The Maxwell stress tensor, generated by the electrostatic potential gradient, also comes, in general, with a deviatoric component. Can highly charged double layers develop shear elasticity and show chemo-mechanical response resembling the solid bar investigated in the work presented here? Another matter is possible coupling between Maxwell and steric stress. This requires a distinction between charge and mass density, which is usually not explicitly made in modified Poisson–Boltzmann theory.58 How are these electromechanical effects reflected in the differential capacitance profile? As mentioned in the Introduction, these questions, which are beginning to be explored,59,60 were the original motivation for the present investigation and are already the subject of ongoing research.23 Furthermore, a continuum mechanics based analysis of stress may also hopefully be useful for the interpretation of the results of surface force balance experiments carried out in electrically active soft matter.61
ACKNOWLEDGMENTS
Daan Frenkel is acknowledged for his interest and critical comments. No external funding was received for the research reported in this paper.
AUTHOR DECLARATIONS
Conflict of Interest
There are no conflicts of interest to disclose.
DATA AVAILABILITY
Data sharing is not applicable to this article as no new data were created or analyzed in this study.