On the deformation dependency of the diffusion flux in solids at large deformations

In material modeling, when dealing with diffusion at large deformations, there are usually two different variants for the diffusion flux: an isotropic law in the current placement and an isotropic law in the reference placement. The first one causes diffusion behavior, which is independent from the initial shape of the body, i.e., it causes a deformation-independent behavior. The second one relates the diffusion solely to the initial shape of the body, which results in a deformation-dependent behavior in the current state. In most of the works in the literature, one of these two possible formulations is chosen arbitrarily. While the modern description of diffusion at large deformations mostly evolved in the last two decades, to our best knowledge, there are no works which discuss or motivate the choice for one of these two versions really in detail. In the present article, we approach the motivation for the choice of the two different types of diffusion flux formulations. We illustrate their characteristics and discuss their application under different circumstances. It is important to note that the deformation dependency which arises from choosing the isotropic reference placement formulation is quite specific and strongly differs from the actual behavior of many materials. We investigate such a case with a more individual deformation dependency based on a very simple artificial microstructure. We determine the properties on the macroscale using representative volume elements within numerical homogenization.


Introduction
Diffusion is quite important when trying to understand the behavior of materials which are exposed to different environmental conditions. Typically, gases diffuse into surface zones of machine parts and lead to chemical aging. In some cases, especially at high temperatures or when dealing with polymer materials, a large deformation framework can be necessary. When modeling a process with an isotropic diffusion law, one has to choose whether the material is supposed to behave isotropic in the reference placement or isotropic in the current placement. The diffusion can be deformation dependent or deformation independent as the result of this choice.
Diffusion or mass transport in the presence of large deformations is mostly known from polymer materials. Often, this is considered in the context of swelling due to the intrusion of small molecules. For polymeric solids, Rajagopal [1] created an overview on articles dealing with diffusion at large deformations. A quite early approach can be found in the work of Durning and Morman [2] in 1993. At the beginning of the century, further articles dealing with this were published [3][4][5][6][7]. Not all of the authors formulated the diffusion flux Communicated by Andreas Öchsner.
J. Voges (B) · F. Duvigneau · D. Juhre Institute of Mechanics, Otto von Guericke University Magdeburg, Universitätsplatz 2, 39106 Magdeburg, Germany E-mail: jannik.voges@ovgu.de based on the gradient of the chemical potential, but those who did, used either the isotropic current placement formulation or the isotropic reference placement formulation. In 2008, a quite well-known or often-cited article was published by Hong et al. [8] which again approached coupled diffusion and large deformation in a polymer material, or more specifically in polymeric gels. In contrast to the other works mentioned before, the authors explicitly motivated the choice of the diffusion flux formulation in their model. They assumed isotropic and deformation-independent diffusion behavior and implemented this with the isotropic current placement formulation. This was adopted, e.g., by Zhang et al. [9] and Bouklas and Huang [10].
Other works, e.g., one by Doi [11] on gel dynamics, defined the mass transport based on pressure gradients (Darcy's law). However, for the paper at hand, we stick to formulations based on the chemical potential.
In 2010, a very detailed derivation of the equations for mass transport under large deformations was done by Duda et al. [12]. They wrote that their purpose was to suit a theory more with the original ideas of Gibbs [13] on heterogeneous substances. They mentioned, e.g., the usage of the gradient of the chemical potential. Similarly as some aforementioned works, they chose the diffusion law to be isotropic in the current placement. Other important works were done by the group of Anand [14][15][16]. In one of their articles [14], they gave a very detailed derivation of a large deformation diffusion theory using the isotropic reference placement formulation. They chose this formulation as it directly follows from their derivation without making additional constitutive assumptions. They also included a comparison with the recent works of Hong et al. [8] and Duda et al. [12]. In two articles [15,16], Anand and colleagues also approached the diffusion and oxidation underlying elasticviscoplastic deformations in turbine blades or high-temperature applications in general. Similar to their other works, in these articles they used an isotropic law in reference placement. The same group also assumed such a formulation for modeling the swelling of silicon electrodes [17,18]. Later in 2015, in an article [19], they briefly stated that they adopt the diffusion flux formulation of Hong et al. [8] and Duda et al. [12], while in previous works, they had preferred the isotropic reference placement formulation.
Another work that includes the isotropic current placement version was done by Drozdov [20] on the swelling of pH-responsive cationic gels. In more recent works by other authors, the isotropic reference placement formulation can be found in a formulation on the swelling of rubber [21] and on the oxidation of silicon carbide fibers [22]. In a work on polyelectric gels [23], both formulations were shown in the theory; however, in the application part of the article, the isotropic reference placement formulation was implemented. Another recent isotropic current placement formulation was used for electro-chemo-mechanical modeling of ion transport [24].
The decision for a specific formulation of the diffusion flux often plays a comparatively small role in the cited works. Most authors do not motivate their choice. Even nowadays, it does not seem to be completely clear which formulation is the more appropriate one for specific applications.
While the cited works mostly used the isotropic diffusion flux formulations just as a small part of their application-oriented models, there are also works which focus on the investigation of the deformation dependency of the diffusion or mass transport. In one example, Masoud and Alexeev [25] determined the effective anisotropic diffusion properties for deformed polymer networks. In another example, Markert approached deformed high-porosity polyurethane foams [26].
In this present article, we compare the current placement formulation and the reference placement formulation in detail and point out the consequences from choosing one instead of the other. Based on this, we assign these formulations to different types of diffusion mechanisms like the movement of atoms in a crystal lattice or the movement of particles in between molecule chains. Furthermore, the choice for the reference placement formulation leads to a specific deformation-dependent diffusion behavior. We compare this with a complete different deformation dependency which directly arises from a microstructure. For this, we perform numerical homogenization for the diffusion flux within the finite element method (FEM).

Basics for diffusion under large deformations
In the following subsections, we introduce basic quantities like the diffusion flux or the chemical potential. Additionally, we briefly describe diffusion mechanisms and the deformation kinematics used in this article.

Diffusion mechanisms
Large deformations do not influence all diffusion mechanisms in the same way. In this section, we want to briefly introduce some specific cases which we later discuss in the presence of large deformations. On the left-hand side in Fig. 1, we illustrate an impurity as a substitutional atom in a crystal lattice as well as a much smaller interstitial atom. The diffusion of both impurity atoms is based on changing places in the crystal lattice. This plays a role for impurities which can influence the properties of, e.g., steel [27]. In the middle drawing, the single atoms (representing also small gas molecules) can move quite freely within long molecule chains which need to be bypassed. One example is diffusing oxygen in polymers, which results in chemical reactions as a part of aging processes. This leads to embrittlement and changed properties in general [28]. The third drawing illustrates a material with two solid phases which can provide different properties for a diffusing gas. The phase distribution determines the effective behavior on a larger scale.

Deformation
We use a typical continuum mechanics notation. It contains the symbol B R for a body in the reference placement. In the current placement, it is B which, in contrast to B R , can change over time. The deformation gradient Here, x are spatial coordinates, while X denote material coordinates. The capital letter in the gradient operator indicates the usage of partial derivatives w.r.t. material coordinates, i.e., using the Lagrangian description. Lower case letters, i.e., grad(. . . ), include partial derivatives w.r.t. spatial coordinates, using the Eulerian description. The same holds for divergence operators. To reduce the amount of indices, mostly capital letters are used for reference placement and lower case letters for current placement. This also holds for volume increments which are related via dv = det(F)dV . In some cases, the distinction is pointed out by the index 'R'. In the whole paper, an orthonormal basis is used with {e 1 , e 2 , e 3 } in the Euclidean space. Since we solely apply homogeneous deformation fields and, additionally, here, we are not interested in stress analysis, the balances of linear and angular momentum and a specific mechanical material law are not introduced.

Conservation of mass
A substance in the domain B has the mass M and the mass density ρ. The temporal change of the mass,Ṁ = dM dt , is equal to the change of the integral of ρ over B. This change is nonzero for nonzero mass fluxes h leaving or entering the domain over its boundary ∂B, Typically, a positive flux h quantifies mass leaving the domain, which, here, results in the minus sign for the case of a surface increment vector da which points outward. M is proportional to the amount of substance N (i.e., the number of atoms or molecules divided by the Avogadro's number). The mass density ρ and the mass flux h are proportional to the molar concentration c and the molar flux j, respectively. Since they all share the molar mass as the proportionality factor, we can divide Eq. (1) by it and geṫ The time derivative of the integral term is The overline is used to clarify that the dot is applied to the whole expression det(F). In Eq. (3), we first converted the domain to the reference placement having a volume which is not time dependent. Then, the product rule is applied and the conversion is reversed. Applying the divergence theorem to the right-hand side in Eq. (2) and combining it with the concentration terms of Eq. (3) give In reference placement, there is no need for conversion of the volume and no product rule needs to be applied. We get Later, we need a local version of Eq. (4). We can writė which simplifies for the equilibrium state to div(j) = 0.

Chemical potential and diffusion flux
The chemical potential μ can be derived from an energy imbalance [29], For convenience, we use the reference placement with the Helmholtz free energy density ψ R and the first Piola-Kirchhoff stress tensor P. Applying the chain rule on ψ R which depends on F and C, then rearranging terms yields For arbitraryḞ andĊ, we can deduce relations for the stress tensor and for the chemical potential. From the latter one, we get Based on the remaining third term in Eq. (9), the diffusion flux can be related to the gradient of the chemical potential in the specific form [29,30] The second-order tensor M is usually called mobility tensor. To satisfy the inequality, it must be positive semi-definite. Local differences of the chemical potentials in neighboring points drive the diffusion. Vectorial fluxes can be transformed with [31] Based on this, for the third term in Eq. (9), we can perform the transformation to analogously arrive at the relation for the current placement flux with another second-order mobility tensor m. For chemical potentials which solely depend on the concentration, we get using the chain rule. In this article, we also use the diffusion coefficient tensors d and D. Note that couplings like concentration-dependent elastic properties or concentration-dependent volume changes can lead to additional variables in μ which then generate additional gradient terms.

Different formulations of isotropy
Both formulations, Eqs. (11) and (14), can imply the same material behavior. The difference results from the choice of m and M. To clarify this, we discuss the isotropic case for both fluxes. Using the scalar mobilities m and M besides the identity tensor 1, we can write We use Eq. (12) for a pullback of the flux j and for a push-forward of the reference placement flux J. Additionally, we use Grad(. . . ) = F T grad(. . . ) and the inverse version of this and get New possibly anisotropic mobility tensorsm andM emerge. They correlate with the inverse of the right Cauchy-Green tensor, C −1 = (F T F) −1 and the left Cauchy-Green tensor, b = FF T . The outcome of equations (11) and (14) is as follows. The Eulerian formulation is based on the gradient using spatial points. The Lagrangian formulation is based on gradients using points which are fixed on the material. The choice for either Eq. (17) or Eq. (18) can cause deformation dependencies. In the following, we investigate the consequences of this choice in detail to establish a basis to clearly differentiate these effects from another effect which includes deformation dependencies as well.
With respect to Fig. 2, we consider two hexahedral blocks: One has its shape due to a deformation (block 'B') and the other one is undeformed, but has the same hexahedral shape (block 'C'). For simplicity, the deformation is homogeneous and isochoric. Both blocks are exposed to a gas.
The crucial difference in the two formulations is that when using Eq. (17), both the blocks 'B' and 'C' develop the same gas concentration profiles. The formulation results in deformation-independent behavior. At a specific distance from the surface, we have the same concentration value in 'B' and 'C', cf. the red point in Fig. 2. The point is approximately at 25% of the long edge length. When we now would remove the A B C Fig. 2 Illustration of different deformed and undeformed blocks with equidistant vertical dashed lines for orientation. The red circles illustrate deformation-independent diffusion behavior, and the blue crosses illustrate deformation-dependent diffusion behavior deformation (which gives block 'A') after the diffusion process, the red point would still be at 25% of the cube edge, but the absolute distance would be smaller. Using Eq. (18), in contrast, will generate a different concentration profile for the deformed block 'B'. It is not the same as in block 'C' anymore. The red points do not apply here; instead, we use the blue crosses. The behavior is deformation dependent. The concentration profile is completely based on the initial shape 'A', even though we still compare 'B' and 'C'. A particle which travels a distance of length L/2 from the right edge of block 'C' would have made double the distance in block 'B'. This is because a particle, which would travel a distance of length L/2 in the initial geometry 'A', would reach the material coordinate located at 50% of the edge length of the cube. This coordinate corresponds to 50% of the long edge in 'B'.
These differences result directly from the choice, either of isotropic mobility in the current placement or of isotropic mobility in the reference placement.

Assigning the formulations to different material types
In this section, we relate the current placement formulation and the reference placement formulation to different material types or to different diffusion mechanisms, respectively.
For all following considerations, we choose a body on the macroscale. First, we assume a crystalline metal or ceramic material and exemplify the diffusion behavior using a lattice of atoms like it is illustrated in Fig. 1 on the left. An elastic deformation changes the distances between the lattice points. Imagine the cube 'A' in Fig. 2 consisting of such a lattice of atoms and an atom travels from one side of the cube to the other side of the cube. For this, a specific amount of lattice point swaps is needed. When the same ion travels from one side to the other side in the elastically deformed body 'B', the amount of lattice point swaps does not change. The lattice points behave similar like material coordinates. It seems natural that the proper macroscopic continuum description of this diffusion process is Eq. (18). This holds for both the substitutional lattice diffusion and the interstitial lattice diffusion. Note that possible changes in the swapping mechanism due to a higher distance between centers of adjacent lattice points are ignored. However, any deviations from the simplified consideration can be covered either in the free energy function ψ R , or in a possibly anisotropic mobility tensor M using the more general Eq. (11) rather than Eq. (18).
Often, when modeling metals at large deformations, the elastic contribution is quite small and most of the deformation can be referred to plastic deformations. We thus now assume that an elastic deformation of a metal is sufficiently small and does not change the diffusion behavior significantly and focus on an additional large plastic deformation. In crystal plasticity, the crystal structure changes due to dislocation glides, dislocation climbs or other mechanisms [31]. However, as an example, having a strongly deformed, but mostly stressfree body implies that the crystal structure has been rearranged. In such a plastically deformed case, for a traveling atom, a larger distance would mean the need of more lattice point swaps. In this case, the distance in the current shape of the body controls the diffusion behavior. The initial shape and thus the deformation are not relevant. The other flux definition, i.e., Eq. (17), seems to be the better choice. Again, any deviations from this simplified consideration can be covered in the free energy function ψ R , or in a possibly anisotropic current placement mobility tensor m with Eq. (14) rather than Eq. (17). This refers to, e.g., newly emerged crystal defects or microcracks. The influence of microcracks in terms of a macroscopic damage variable on the diffusion behavior was investigated in several works [32,33]. Now, referring to the middle drawing of Fig. 1, for a second consideration, we leave the topic of crystal lattices to focus on free atoms or molecules moving in a medium which is occupied by obstacles. As an example, these obstacles can be long polymer chains. Here, the movement of the diffusing particles is not restricted in a way like in the crystal structure. For a particle which is moving through the voids and gaps of such an amorphous structure, an elastic deformation of the structure simply changes the distance which it has to travel to cross the whole block. Like in the plastic deformed crystal lattice, the current placement type of the diffusion flux seems to be the logical choice, i.e., Eq. (17). Deviations from the simplification, i.e., when expecting an influence of the rearranged molecule chains on the effective diffusion behavior, once again, can be covered in the mobility tensor m or in the free energy function ψ R . Particularly, when dealing with larger diffusing molecules in polymer materials, the process can be more complex. Just to mention one mechanism, a deformation can enable the movement of a molecule by opening up new space between molecule chains [34].

Comparison of deformation dependencies
In this section, we first quantify the deformation dependency which arises from choosing the isotropic reference placement formulation. In the second part of this section, we calculate a different type of deformation dependency which originates from the deformation of a microstructure. The same exemplary deformation gradient is used for both considerations.

Isotropic reference placement formulation
Using Eq. (20) and rewriting the expression in terms of the gradient of the concentration, cf. Eq. (15), we get whered is a diffusion coefficient tensor in current placement. It includes the deformation dependency which arises from choosing the isotropic reference placement formulation as the material law. We use a deformation which we illustrate in Fig. 2. For the orthonormal basis, the coefficient matrix of the corresponding deformation gradient is The coefficient matrix ofd then isd Without specifying M and ∂μ ∂c , we can observe thatd 11 = 8d 22 = 8d 33 which implies a strong anisotropy in the current state. However, this deformation dependency is very specific since it is fully determined by the push forward in Eq. (20). The anisotropy then is the same for all materials. In the following section, we investigate a case of a more individual deformation dependency which is significantly different.

Artificial microstructure
In the drawing on the right-hand side in Fig. 1, we show a material consisting of two different solid phases on the microscale. Imagine that on an even smaller scale, each of these two phases behaves like it consists of long molecule chains (middle drawing in Fig. 1). Based on our previous argumentations, the diffusion process of a gas which is moving in these two phases is thus assumed to be deformation independent, i.e., we use the Fig. 3 Comparison of a deformed specimen '2' and an undeformed specimen '3' with the same shape. The circles indicate detail views of an artificial microstructure which is deformed in '2' and undeformed in '3'. In '2b' and '3b', representative volume elements are shown. The microstructure is based on a red phase and a white phase isotropic current placement formulation for both. Even with this foundation, the effective macroscopic material properties can still be deformation dependent due to changes on the microscale when one phase allows faster diffusion than the other one. This is similar to shapes like voids and tunnels which open up or close. In the following, we investigate such a case by using a simple artificial microstructure.
In Fig. 3, we show three different blocks. Imagine a big sample of a material with an underlying microstructure and the block '1' and the block '3' are cut out of this. The block '2' is achieved by a deformation of block '1' and has the same outer geometric shape as '3'. The artificial microstructure is indicated by the red and the white regions in the detail views. Since only one is deformed, the microstructures look different. In this generic example, we assume that both the deformation on the macroscale and the deformation on the microscale are homogeneous. Furthermore, the white phase represents a material which provides an increased diffusivity compared to the rest of the block. Even though these regions with increased diffusivity are solid material, figuratively, we label them as 'tunnel'. The low-diffusivity red region is labeled as 'bulk'. We use as the flux in each phase. It is based on Eq. (15) with an isotropic diffusion coefficient tensor d = k1. The scalar diffusion coefficient k determines the diffusivity for the tunnel material or for the bulk material, depending on the coordinate x. For our quite arbitrary case, we use the relation k tunnel = 10 k bulk . In Fig. 3, we also illustrate a representative volume element (RVE) for each of the blocks '2' and '3' and label them as '2b' and '3b'. To compare the fluxes inside the RVEs, in the following paragraphs, we solve the partial differential equation Eq. (7) numerically. We insert Eq. (24) and get 0 = div(j) = − div(k grad(c)).
For periodic boundary conditions [35], we can write where the vectorial quantity L represents the length of periodicity. The quantity grad(c) is an effective gradient which we apply on the RVEs. Multiplying this with L yields a concentration change from one face to the opposite face. The minus in Eq. (27) results from the orientation of the surface normals which usually point outward.
To compute a numerical solution of the problem with FEM, we use the software FEAP [36]. For qualitative rather than quantitative results, we choose a quite rough discretization with 64 × 64 × 64 hexahedral linear elements in a uniform mesh. The material distribution is defined on the integration-point level [37][38][39], rather than by fitting the element boundaries to the material boundaries. This engages easy implementation, especially for rounded corners which we used to replace any sharp corners inside the RVE to avoid flux singularities.
In Fig. 4, we plot the resulting vectorial fluxes for different effective concentration gradients along the RVEs. For the gradients in this illustration, we choose the negative e 1 -direction and the negative e 2 -direction, each with a positive scalar value α. We omit a plot for the negative e 2 -direction for '3b', as in this case, due to symmetry, the results would be just a rotated version of the shown one. For '2b', we consider both directions. Having a look at the plot on the left-hand side in Fig. 4, we can observe that along the main diffusion direction, three segments, i.e., 'I', 'II' and 'III', emerge. The first one and the third one characteristically consist of a slender region which has a high diffusivity in between a low-diffusivity area and thus localizes huge diffusion fluxes. In segment two, the magnitudes of the diffusion fluxes are smaller, but arranged over a larger part of Fig. 4 Transparent material interfaces within uniform spatial distribution of diffusion flux vectors, colored and scaled with respect to their magnitude value. Shown are the effective concentration gradient in negative e 1 -direction for RVE '2b' (left) and the effective concentration gradient in negative e 2 -direction for both RVE '2b' (middle) and RVE '3b' (right). A second one for '3b' is dropped due to symmetry the cross section. The transition of the diffusion flux distributions from segment to segment is smooth. We also have such segments in the other diffusion direction, i.e., in the middle plot of Fig. 4. Here, in contrast, these segments are small, while the transitions are more dominant. Compared to the other two plots, the segments in the '3b' RVE on the right-hand side appear to be developed moderately, i.e., less than on the left-hand side plot and more than on the middle plot. When using effective diffusion properties to represent the diffusion behavior in the RVEs, we can observe an analogous order. In the following, we quantify these properties by performing homogenization.
Following a basic numerical homogenization procedure for vectorial fluxes [35], we can write the equation which includes the effective flux j and the effective diffusion coefficient tensor d besides the already introduced effective gradient. After performing an FEM simulation, we can calculate j by volume averaging the flux field. Using different effective gradients then allows to determine d for a specific RVE with Eq. (28). For this, we use the gradients αe 1 , αe 2 , αe 3 and combinations like α(e 1 + e 2 ). For α, we use 1 mol μm 4 which is concentration per length scale. Due to linearity, a ten times bigger gradient would lead to a ten times bigger flux which then gives the same d. The edge lengths of the cube '3b' are 1 μm. The deformed edge lengths, i.e., those of '2b', are 2 μm and 0.707 μm according to Eq. (22).
The resulting coefficient matrices for our chosen basis are Note that other relations than k tunnel = 10 k bulk will give different results. With respect to Eq. (29), we can observe that the diffusion flux coefficient tensor for the specific deformed microstructure in '2b' of Fig. 3  The deformation dependency, which results here from the microstructure, has a completely different degree compared to the one, which is implied by the isotropic reference placement formulation in Sect. 5.1 which has the huge anisotropyd 11 = 8d 22 = 8d 33 .
We only performed this procedure for one specific deformation gradient. Others can be used to determine the tensor d 2b as a function of F. Then, d 2b can be inserted in Eq. (15) to use the specific behavior in a macroscale computation.

Conclusion
As described in "Introduction", the current literature on diffusion at large deformations lacks in detailed motivations for choosing a specific diffusion flux formulation. Often one is interested in a simple and thus isotropic material law. For this, however, usually one has to choose between an isotropic current placement formulation or an isotropic reference placement formulation. In a deformed body, the latter one leads to a specific deformation dependency which is entirely determined by a push forward of the material law. The other one gives deformation-independent behavior. We discussed these two formulations and related them to different diffusion mechanisms under large elastic deformations. We think that for metals or ceramics, or more specifically, in the presence of lattice diffusion, the isotropic reference placement is the appropriate choice under large deformations. However, the elastic deformation is usually quite small in comparison with possible plastic deformations, e.g., in metals. For such cases with small elastic deformations and large plastic deformations, the isotropic current placement formulation can be more appropriate. Based on our considerations, we also state that this formulation is reasonable for the diffusion of particles in amorphous structures which coincides with the choice of Hong et al. [8]. However, these isotropic formulations are simplifications, and thus, sometimes, adjustments in the free energy function and in the mobility are necessary, which can lead to anisotropic formulations. This can involve a deformation dependency which is individual and thus usually different from the one which arises from choosing the isotropic reference placement formulation. We used an exemplary deformation gradient to compare these kinds of deformation dependencies. In our example, the more individual one originated from an underlying simple artificial microstructure. For this, we homogenized the diffusion properties using RVEs. The resulting properties differ strongly from the properties, which result from the isotropic reference placement formulation.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by Projekt DEAL. We gratefully acknowledge the financial support for this work by the Research Training Group 1554/2: Micro-Macro-Interactions in Structured Media and Particle Systems of the DFG (German Research Foundation).

Conflict of interest
The authors declare that they have no conflict of interest.