Computational homogenization of transient chemo-mechanical processes based on a variational minimization principle

We present a variational framework for the computational homogenization of chemo-mechanical processes of soft porous materials. The multiscale variational framework is based on a minimization principle with deformation map and solvent flux acting as independent variables. At the microscopic scale we assume the existence of periodic representative volume elements (RVEs) that are linked to the macroscopic scale via first-order scale transition. In this context, the macroscopic problem is considered to be homogeneous in nature and is thus solved at a single macroscopic material point. The microscopic problem is however assumed to be heterogeneous in nature and thus calls for spatial discretization of the underlying RVE. Here, we employ Raviart–Thomas finite elements and thus arrive at a conforming finite-element formulation of the problem. We present a sequence of numerical examples to demonstrate the capabilities of the multiscale formulation and to discuss a number of fundamental effects.


Introduction
The continuous advancement of technological innovations in the fields of digitalization and automation leads to increased demands of smart and multifunctional materials. In that context, the design of associated materials with tailor-made properties is paramount. Classical examples cover solids with coupled electro-and magneto-mechanical response such as ferroelectrics or magnetostrictives. More recently, several advanced materials exhibiting chemo-mechanical coupling have entered the scene and experience increased attention since then. Related materials show interactions between fluid flow and deformations and play an important role in emerging applications such as lithium-ion batteries (Wang et al. [40]), heterogeneous concretes (Wang and Ueda [41]), engineered biological tissues (Truskey et al. [38]), or fiber-reinforced superabsorbent hydrogels (Chen and Park [7]).
In order to allow for the theoretical development of related materials and structures, reliable continuum-mechanical models are needed. Related formulations have to take into account the mutual couplings between the fluidic diffusion of some solvent phase and the mechanical deformation of some solid phase in a single multiphase material. In fact, the development of continuum-mechanical models and associated numerical implementations have seen a lot of advancements recently. We refer to the works of Hong et al. [19,20], Chester and Anand [8], Nilenius et al. [28], Miehe et al. [25], Ehlers and Wagner [11] and Böger et al. [4] among many others.
In addition to that-when targeting the design of materials-the need for multiscale continuum approaches becomes evident. Associated homogenization techniques involving continuum microstructures trace back to the early works of Voigt [39], Reuss [33], Hill [18]. For the realization of multiscale numerical simulations we refer to Miehe et al. [26], Feyel and Chaboche [12]. A common feature of the above approaches is the identification of a representative volume element (RVE) that reflects the heterogeneity of a given microstructure. Here, one often assumes separation of length scales by postulating that typical length scales of the RVE are much smaller than typical length scales of a corresponding macroscopic problem. This then renders the notion of so-called first-order homogenization schemes. 1 Extensions of first-order homogenization towards the incorporation of heat conduction have been proposed by Özdemir et al. [29], Temizer and Wriggers [37], Temizer [36] and Chatzigeorgiou et al. [6]. Here, the transient behavior of the material is considered only at the macroscopic level, while at the microscale the problem is assumed to be stationary. In contrast to that, Larsson et al. [24] have proposed the variationally consistent homogenization of heat conduction that takes into account the transient behavior also at the microscopic level. 2 The latter approach leads to size-dependent macroscopic response. As could be shown by the authors, the size dependence vanishes when the size of the RVE approaches zero. While such effect seems in accordance with the idea of scale separation, it does not occur in the first-order homogenization of stationary problems. Related effects and consequences have been discussed by Kaessmair and Steinmann [21] who proposed a transient homogenization framework for coupled chemo-mechanical problems.
Close to the formulations of Larsson et al. [24] and Kaessmair and Steinmann [21] we propose a variationally consistent homogenization approach to the coupled chemomechanics of transient diffusion-deformation processes in soft, porous solids. In doing so, we critically revise the influence of the RVE size on the effective macroscopic response of a transient microscopic problem. In contrast to the above-mentioned homogenization schemes of thermo-and chemo-mechanics, we consider the computational homogenization of diffusion-deformation processes within a minimization-based variational formulation. The latter is adopted from the ideas of Miehe et al. [25] and Böger et al. [4].
In the minimization-based formulation, the deformation map and the solvent-volume flux act as the independent fields. Associated formulations have several advantages compared with direct formulations (that are based on balance equations) and standard saddle-point formulations (that usually exploit the chemical potential as independent field next to the displacement): 1. In contrast to direct approaches, variational structures lead to compact representations involving symmetric system matrices [4,25]. 2. For two-dimensional problems it could be shown that the minimization formulation is computationally more efficient than both direct approaches and saddle-point formulations [25,35]. 3. A minimization formulation is not restricted by the inf-sup condition, which is an inherent problem of saddle-point formulations [23,35]. 4. A minimization formulation is favorable for the investigation of structural instabilities since the coupled stiffness matrix is positive definite unless an instability occurs [14,31]. 5. In the context of homogenization methods, minimization structures can be embedded into Hashin-Shtrikman variational principles and thus allow the computation of bounds of effective material properties [17,27].
Next to that, we should however note that the above advantages come at the cost of a more involved spatial discretization of the solvent-volume flux based on H(Div)conforming finite elements. For the latter we refer to Miehe et al. [25], Böger et al. [4], Teichtmeister et al. [35].
The structure of the present contribution is as follows. In the section "Variational homogenization of diffusion-deformation processes", we start with the definition of a rate-type variational principle of homogenization based on the deformation map and the solvent-volume flux as independent field variables. The rate-type principle is then translated into an incremental variational principle of homogenization. The latter comes along with discretization in space and time. For time integration, we employ an implicit Euler scheme. The spatial discretization is realized by means of conforming Raviart-Thomastype finite elements. Based on the discrete setting, we are able to determine all relevant macroscopic quantities including the algorithmically consistent material moduli. In the section "Size effects in the homogenization of transient diffusion", we briefly comment on the size dependence of the homogenization problem involving the transient diffusion at microscale. This discussion serves as motivation for the set of numerical examples to be covered in the section "Representative numerical examples". There we present a sequence of numerical studies of the effective response of two-phase porous composites. A conclusion of the present work will finally be given in the section "Summary".

Variational homogenization of diffusion-deformation processes
In the present section, we discuss the computational homogenization of diffusiondeformation processes based on a rate-type variational formulation. We consider a variationally consistent approach close to the ideas of Larsson et al. [24]. However, in contrast to the latter work, our goal is to develop a homogenization scheme that is based on a minimization principle. The spatial discretization of the coupled problem will be realized within a conforming finite-element formulation based on Raviart-Thomas-type finite elements. In order not to overload the presentation, we will keep the motivation for the use of certain functionals as concise as possible. For more information on related topics the interested reader is referred to Hong et al. [19,20] and Böger et al. [4].

Constitutive state variables at the microscopic level
In the context of a minimization formulation, a diffusion-deformation process can be described by two independent fields given by the deformation map ϕ and the solventvolume flux H having the following properties (Gurtin et al. [16], Miehe et al. [25], Böger et al. [4]) ϕ : and H : where X denotes a material point at the microscopic reference configuration B 0 , which is mapped to the current configuration B t by the deformation map ϕ. The solvent-volume flux H describes the flow of the solvent volume relative to the solid skeleton. Such interpretation follows Biot's approach to the modeling of porous media [3]. For a detailed derivation of the equations used in the present contribution we refer to Coussy et al. [9].
Furthermore, we introduce the concentration of solvent volume describing the amount of solvent molecules in a referential infinitesimal volume element s : with s = νc. Here, ν denotes the volume of a single solvent molecule and c is the number of these molecules. The solvent-volume concentration s can be related to the solvent flux by considering the conservation of solvent-volume concentration (Gurtin et al. [16]). It states that the change of solvent-volume concentration in a given arbitrary domain of a body is equal to the flux of solvent volume across the surface of that domain. Locally, this can be expressed aṡ Considering a first-order approach of homogenization, the independent fields can be decomposed into microscopic and macroscopic contributions where ϕ and H denote microscopic fluctuations of the deformation map and the solventvolume flux, respectively. We refer to Fig. 1 for a graphical illustration. Further, considering a Taylor expansion of the macroscopic contributions up to the linear terms yields 3 l macro l micro π(F , , s) Graphical illustration of the computational homogenization of chemo-mechanical processes. The macroscopic dual fields and moduli at a macroscopic material point X are determined by solving a microscopic boundary value problem. The domain of the microscopic problem is chosen such that it represents the heterogeneity of the microstructure in a representative manner and is thus referred to as representative volume element (RVE). In a minimization-based variational formulation the RVE is driven by the deformation gradient F as well as the solvent-volume flux H and its divergence Div H where F = Grad ϕ is the macroscopic deformation gradient andṡ = −Div H is the rate of the macroscopic solvent-volume concentration. In what follows, we assume that the origin of the coordinate system is located at the geometric center of the RVE, i.e., B 0 X dV = 0. From (4) we obtain the deformation gradient F and the divergence of the solvent flux Div H as F = F + F and Div H = Div H + Div H .
The latter equation implies that the solvent-volume concentration can be decomposed into fluctuative microscopic and constant macroscopic parts as s = s + s. Integrating the equations in (6) over the domain of the RVE with |B 0 | = vol(B 0 ) yields where we have considered the following relations The former equation is satisfied for continuous fluctuations of the deformation map along the whole microstructure, where continuity also relates to the transition across the (arbitrary) boundary of a periodic RVE. The restriction w.r.t. the continuity across an At the respective material points X + and X − on these boundaries, the deformation map is considered to be periodic (ϕ(X + ) = ϕ(X − )) and the normal solvent-volume flux H is considered to be anti-periodic (H(X + ) · N + = −H(X − ) · N − ) RVE's boundary is usually expressed as ϕ = 0 on ∂B 0 , where (·) = (·) + − (·) − denotes the jump of a quantity (·) across the boundary of an RVE. Analogously equation (8) 2 is satisfied for continuous normal projections of the solvent-volume flux, i.e., anti-periodic normal projections H · N = 0 on ∂B 0 . We refer to Fig. 2 for a graphical illustration.

Constitutive functions at the microscopic level
As diffusion-deformation processes are dissipative in nature, they can be modeled by a set of two constitutive functions given by an energy-storage functionψ and a dissipationpotential functionφ. The functionψ(F , s) models the energy storage due to the deformation of the material and the influx of the solvent molecules. It is assumed to have the following additive form (Hong et al. [19,20], Böger et al. [4]) whereψ mech ,ψ chem andψ coup denote the mechanical, chemical and coupled chemomechanical contributions. For the modeling of hydrogels, a neo-Hookean material model can be considered for the mechanical energyψ mech . The chemical energyψ chem is usually governed by a Flory-Rehner-type constitutive function (Flory and Rehner [13]). The term ψ coup (J, s) in the energy-storage function models the coupling between deformation and solvent volume and by relating the volume change of the material to the solvent-volume concentration s. In the variational minimization formulation to be developed below, the dissipation potentialφ(H ; F , s) is a convex function of the solvent flux H . 4 The forms of these functions and their relevance for the specific problems at hand will be discussed in the following sections.

Rate-type minimization principle of computational homogenization
The macroscopic potential density π(˙ ϕ, H ;Ḟ , H ) of the homogenization problem at a given macroscopic constitutive state {Ḟ ,ṡ = −Div[H ]} can be described by the mini- This variational principle can be considered as a generalized Hill-Mandel condition and involves the minimization of the volume-averaged microscopic potential density in the admissible function spaceṡ

Incremental variational principle at the microscale
In the numerical setting, the above given rate-type variational principle (10) is solved incrementally at discrete times 0 < t n+1 ≤ T , where T denotes the ending time of a process. It is assumed that at the beginning of a time interval [t n , t n+1 ] all variables are known. In the following, we consider the notation (·) n for discrete variables at time t n and drop the subscript for the variables at time t n+1 , i.e., (·) := (·) n+1 . Considering an implicit Euler time integration with a time step τ = t n+1 − t n , we replace (10) with an alternative incremental macroscopic potential density (13) in terms of the incremental microscopic potential density Here, the solvent-volume concentration is determined via implicit Euler integration such that s = s n −τ Div H . Analogously, we determine s = s n −τ Div H at the macroscale. Note that in (14) the dissipation-potential function takes into account known state variables {F n , s n } at time t n in order to preserve a variationally consistent structure of the problem (Miehe et al. [25], Böger et al. [4]). The admissible spaces in the incremental setting are defined as

Euler-Lagrange equations and linearization of the variational formulation
The necessary condition of the incremental variational formulation (13) yields where δ ϕ ∈ W ϕ and δ H ∈ W H . Using integral theorems as well as employing the periodicity conditions (15), we arrive at the incremental Euler-Lagrange equations In (17) we identify the constitutive equations for the microscopic dual fields, i.e., the first Piola-Kirchoff stress tensor P := ∂ Fψ and the chemical potential of the solvent-volume concentration μ := ∂ sψ . The equation (17)  Linearization of (16) leads to the following expression including the incremental microscopic moduli tensor In the following sections, we will take advantage of (16) and (18) when we consider the finite-element implementation of the homogenization problem.

Finite-element implementation of the homogenization procedure
For the numerical solution of the microscopic diffusion-deformation problem, conforming quadrilateral Q 1 RT 0 Raviart-Thomas-type finite elements are implemented [32]. This choice is due to the required function space H(Div, B 0 ) of the solvent-volume flux H . An alternative non-conforming discretization on the basis of standard Q 1 Q 1 finite elements usually yields non-physical behavior. For recent developments of equivalent minimization principles incorporating a node-based finite-element discretization combined with reduced integration we refer to [35].
In the framework of the Raviart-Thomas-type finite-element discretization, the fluctuations of the deformation map ϕ are approximated continuously using bilinear Q 1 shape functions, which are related to the nodal degrees of freedom d e ϕ of a finite element. The fluctuations of the solvent-volume flux H are approximated considering linear RT 0 shape functions, which are related to the fluctuations of the flux degrees of freedom d e H across the edges of the quadrilateral element. We refer to Fig. 3a for a qualitative illustration.
Note that only the normal component of the solvent-volume flux h K = E K H h · N e dA at an edge of the element E K is considered as degree of freedom. We refer to Raviart Raviart-Thomas finite element and its mapping from the parametric to the physical space. a The deformation map ϕ is continuously interpolated using node-based Q 1 shape functions. The solvent-volume flux H is approximated using edge-based RT 0 shape functions. b To obtain the RT 0 shape functions and their derivatives in the physical space, a Piola transformation from the parametric space according to (21) is considered The factor 1/2 is due to the length of the element edge in the parametric space such that d K and Thomas [32], Brezzi and Fortin [5] for theoretical aspects and to Schwarz et al. [34], Anjam and Valdman [1], Böger et al. [4], Teichtmeister et al. [35] for discussions of the numerical implementation. The interpolation of the independent fields arises as where the scalar shape functions N I ϕ are defined at the node I and the vectorial shape functions N K H correspond to the edge K . While the former can be directly defined in a so-called parametric space A = [−1, 1] × [−1, 1], the latter have to be computed via a Piola transformation of the corresponding shape functions N K H in that parametric space (from here on the notation (·) indicates quantities (·) formulated in the parametric space). Similarly, d I ϕ denotes the array of the displacement components at the node I and d K H corresponds to the normal component of the solvent-volume flux at the edge K . An illustration of the Raviart-Thomas shape functions in the parametric space A is given in Fig. 4.
The RT 0 shape functions considered in the present contribution are given in the parametric space as The degrees of freedom need to be invariant under the change of coordinates from the parametric space to the physical space, i.e., d K Thus, we apply a Piola transformation to map the vectorial shape functions and their derivatives from the parametric space A to the reference element B e 0 . We refer to Brezzi and Fortin [5] for more details and to Fig. 3b for a graphical illustration. It follows where X h : A → B e 0 refers to a map from the parametric space to a reference finite element and is determined by X h (ξ) = n node I=1 N I ϕ (ξ)X I in terms of the nodal coordinates X I of the considered finite element.
Furthermore, it is important to point out that, in order to have a consistent finite-element implementation for the flux degrees of freedom, we need to make sure that the degrees of freedom on a common edge of two adjacent finite elements are identical. In order to achieve this, we assume a positive orientation for an edge in a finite element B e whenever the node numbers of the vertices of that edge increase in counterclockwise direction. In contrast, when the node number decreases, we augment the RT 0 shape function of that edge with a negative sign. As a consequence, an outflux at an edge of a particular finite element will always be associated with influx across the same edge of a neighboring finite element. We refer to Schwarz et al. [34], Anjam and Valdman [1], Böger et al. [4], Teichtmeister et al. [35] for more details on the numerical implementation.
Based on (19) we obtain the fluctuations of the deformation gradient and of the divergence of the flux fields as Having the above discretization at hand, the necessary condition of the incremental minimization principle (16) in the finite-element setting can be written in form of the microscopic residuum vector where the first Piola-Kirchhoff stress tensor P is implemented as a vector, i.e., P = [P 11 , P 22 , P 12 , P 21 ] T in two dimensions. The corresponding mechanical and coupled moduli tensors follow a similar structure. The tangent matrix of the linearized system of equations reads A converged microscopic state is obtained by solving the system of equations (16) iteratively using, for example, a Newton-Raphson method

Effective macroscopic dual fields and moduli tensor
In the numerical setting, taking the variation of π τ at the macroscopic material point yields where R := π τ ,d = n elem A e=1 π τ ,e ,d e vanishes at the solution point of the microscopic boundaryvalue problem. As a consequence, the macroscopic first Piola-Kirchhoff stress tensor P and the negative gradient of the chemical potential M (the latter being the driving force of the macroscopic solvent-volume flux) are determined as averages of the corresponding microscopic analogues However, the macroscopic chemical potential follows from a different averaging process We observe that μ does not only depend on the average of the microscopic chemical potential but also on its gradient and the spatial dimensions of the RVE. As a consequence, the macroscopic effective response depends on the size of the RVE, which will be discussed in detail in the next section.
Considering a continuous setting and exploiting tensorial manipulations in (27) and (28), we can express the effective dual fields in terms of the surface integrals over the boundary of the RVE as Adopting the notation F := [F , Div H , H ] T , we can write the first-order terms in the linearization of (26) at the macroscale in a compact form as where Combining (23), (26) and (30), we obtain an elimination equation at the macroscopic level for the microscopic degrees of freedom, that is Considering the results of the latter equation in (30), we obtain wherein we have identified the macroscopic coupled moduli Now having determined the macroscopic dual fields (27) and (28) as well as the macroscopic moduli (35), we could solve a macroscopic boundary value problem under prescribed boundary and initial conditions by using a conforming finite-element discretization analogous to the section "Finite-element implementation of the homogenization procedure".

Size effects in the homogenization of transient diffusion
In equations (28), (31) and (32), we observe terms that are scaled with the coordinates of the RVE. This gives rise to a size-dependent effective response at the macroscale. To have more insight into the problem, we simplify the above described formulation (10) to the diffusion in rigid solids by neglecting any elastic effects. Then, the problem at hand reads To determine the macroscopic response of the material, we take the variation of (36) and arrive at the definition of the macroscopic chemical potential Considering the dissipation-potential functionφ(H ) = 1 M H 2 for the diffusion process in rigid solids with mobility parameter 6 M, we can reformulate the size-dependent term μ size of the chemical potential μ in (37) by considering (4) and (5) as Furthermore, after some simple manipulations, we obtain We observe that the macroscopic effective response depends on the size of the RVE and the mobility parameter M. Note that the first and second term on the right-hand side of (39) vanish for specific realizations of RVEs. 7 In general, however, we expect variations among the effective responses depending on the particular selection of RVEs from periodic microstructures. Due to the size dependence, the effective response of a single unit cell D differs from the effective response of an ensemble of identical unit cells nD, n ∈ N 3 . In the following, we provide some numerical examples that describe the mentioned size effect for RVEs with different sizes and material parameters. 6 The mobility parameter M can be computed via the Stokes-Einstein relation where ν is the volume of a single solvent molecule, D is the diffusivity, k B is the Boltzmann constant, T is the absolute temperature, R is the radius of the solvent molecules and η is the viscosity of the solvent [10,20]. 7 One such realization is given by a symmetric RVE that is parameterized with coordinate axes that originate from the spatial center of the RVE.

Representative numerical examples
In order to demonstrate the capabilities of the variational homogenization procedure discussed above, we present some numerical examples. The first numerical examples consider pure diffusion phenomena and investigate the influence of material parameters, geometry as well as the size of RVEs on the effective, macroscopic response of composite materials (see, the section "Investigation of the macroscopic properties of porous rigid solids"). In a second study, we extend the formulation by elastic effects of the solid matrix and couple it to the diffusion processes. In doing so, we will investigate swelling-induced pattern transformations of periodic hydrogels (see, the section "Coupled diffusion-deformation processes in periodic hydrogels").
In what follows, we assume that the relaxation time of the material is much longer than the solvent's diffusion time. We thus neglect any viscoelastic effects of the material and consider standard Fickian-type diffusion [8,19,20]. For an implementation that takes into account the viscoelasticity of the material we refer to Govindjee and Simo [15].

Investigation of the macroscopic properties of porous rigid solids
In this subsection, we analyze the influence of material parameters and the size of RVEs on the macroscopic response of two-phase materials. We will find that the macroscopic mobility parameter lies within the classical Voigt and Reuss bounds for vanishing RVE size.

Microscopic constitutive functions
In the section "Constitutive functions at the microscopic level", we have discussed the basic forms of the energy-storage and dissipation-potential function. In case of a pure diffusion processes (for ϕ = const. with F = 1), the energy-storage function reduces tô ψ(s) =ψ chem (s). As our main focus is to investigate fundamental effects, we select the simple quadratic function (Kaessmair and Steinmann [21]) where A is a chemical material parameter. As dissipation-potential function we again considerφ(H ) = 1 M H 2 with the mobility parameter M.

Description of the problem
In the following, we consider periodic microstructures with circular inclusions. For such microstructures it is straightforward to indentify a periodic unit cell (see Fig. 5a). In a typical homogenization framework, the effective response of the material could be computed uniquely by considering an RVE that resembles exactly that unit cell. However, as we could see in the section "Size effects in the homogenization of transient diffusion", the homogenization of transient diffusion processes comes along with well-known size effects in the sense that the physical size of the considered RVE influences its effective response. In order to analyze this size effect, we will alter the size of associated RVEs in the following way. On the one hand, we will change the size by adjusting the lateral length of the underlying unit cell. On the other hand, we will change the size by composing an RVE of a certain number of unit cells. In both cases the RVEs will be generated from a primitive square-shaped unit cell having the lateral length of 1 mm. We adopt the following name scaling of the dimensions of the primitive unit cell and mm indicates a pattern of m × m primitive unit cells. For example, 2RVE11 and RVE22 both represent an RVE with side length of 2 mm, but while the former is realized by scaling the primitive unit cell by a factor of 2, the latter is made up of 2 × 2 primitive unit cells. Please refer to Fig. 8 for a corresponding visualization. In all following simulations, we consider a volume fraction of the inclusions of f 0 = π/16. In a purely diffusive process, the RVEs can be driven by the macroscopic solvent-volume flux H and its divergence Div H . As the size effect is due to the linear terms appearing in the solvent flux vector in (5) 2 (see also (39)), we take into account the following macroscopic loading conditions where γ is a loading parameter. In the current case, we assume that the divergence of the macroscopic solvent-volume flux increases linearly with time, i.e, γ ∝ t, see

Influence of mobility parameters on the effective macroscopic response
We now investigate the influence of mobility parameters of the individual constituents on the macroscopic response of the two-phase material.
In a first step, we consider different kinds of square-shaped RVEs with constant mobility parameter M mat = 0.1 mm 4 /(Ns) of the matrix and altering mobility parameters of the inclusions.
In Fig. 6a we observe discrepancies in the results depending on the realization of the respective RVE. An increasing size of the RVE leads to a higher magnitude of the effective chemical potential. This effect is due to the last term in (39), which arises as the result of the transient nature of the problem. We refer to Larsson et al. [24] for a similar analysis. As expected from (39), the discrepancies of the effective responses of different RVEs decrease when we increase the mobility parameter of the inclusions. However, since we keep the mobility parameter of the matrix unchanged at a comparably low level (M mat = 0.1 mm 4 /(Ns)) the size effect is still pronounced. Further comparing the responses of 2RVE11 and RVE22, we observe that when the mobility parameter of the matrix and the inclusion have a similar magnitude, the corresponding effective responses nearly coincide (Fig. 6b). In fact, if we would assume M mat ≡ M incl , the responses would be identical since then the RVEs are homogeneous.
In Fig. 7, we show the influence of the mobility parameter of the matrix material on the effective response for fixed mobility parameter of the inclusions M incl = 10 −4 mm 4 /(Ns). We observe decreasing discrepancies in the effective responses with increasing mobility parameter of the matrix. This effect can be attributed to the higher volume fraction of the matrix.
Contour plots of the microscopic chemical potential μ for different realizations of RVEs are shown in Fig. 8. Remarkable differences in the microscopic response of the RVEs are evident. This confirms that the underlying periodicity cannot be used to reduce computations over an enlarged RVE to computations over a unit-cell RVE.
Above we have demonstrated that the effective chemical potential strongly depends on the mobility parameters of the individual phases as well as on the size and realization of the RVE. We now take a close look at the evolution of the chemical potential for various RVEs over time (Fig. 9). In analogy to the previous examples, we consider square RVEs that are either made up of a certain number of equal-sized unit cells or made up of a scaled-up unit cell. The macroscopic loading is again formulated in terms of the loading parameter γ , which will first be linearly decreased to a minimum value and then increased to zero. After that, it will be kept fixed. Please refer to Fig. 5c for a graphical illustration, in which the negative sign of the loading parameter is related with an influx of the solvent. For time incrementation we consider the time step τ = 0.1 s.
In Fig. 9a we observe that the different RVEs show the expected differences in their initial effective response. However, after a certain period of time the responses of all RVEs converge to the same stationary solution, where low mobility parameters of the inclusions lead to longer relaxation times (compare Fig. 9a-c). We further observe that an increase of the mobility parameter of the inclusions has negligible influence on the behavior under loading, but strongly influences the unloading phase (compare Fig. 9a,d with Fig. 9b,e respectively). In contrast to that, changing the mobility parameter of the matrix strongly influences the loading phase, but shows qualitatively similar trends in the unloading regime (compare Fig. 9a,d with Fig. 9c,f respectively). As a final observation we note that the stationary solutions of different RVEs converge to the same values, even when the material parameters of the individual phases are different. (compare Fig. 9b with Fig. 9e). This behavior could be expected from the size-dependent contribution of the effective chemical potential given in (39). There, the first term has been assumed zero per definition and the last term of the right-hand side is zero after full unloading. As a consequence, the second term fades out over time due to internal redistributions of RVE88 RVE88 Fig. 7 Effective chemical potential depending on the size of the RVE and the mobility parameter of the matrix. The evolution of the macroscopic chemical potential μ for different realizations of RVEs is depicted as a function of the divergence of the macroscopic solvent-volume flux Div H for selected mobility parameters of the matrix given by a M mat = 0.1 mm 4 /(Ns) and b M mat = 5.0 mm 4 /(Ns). Like in the previous study (Fig. 6) we observe that the effective chemical potential μ depends strongly on the size of the RVE. The observed discrepancies decrease with increasing mobility parameter of the matrix  the solvent such that finally H → 0 and thus μ size → 0. In the stationary state, the sizeindependent term of the macroscopic chemical potential (37) is identical for all considered RVEs since we assume the same chemical parameters for both phases in the microscopic chemical energy (40).
Furthermore, for vanishing size of the RVE, the macroscopic mobility parameter ∂ 2 HH π lies firmly within the Voigt and Reuss bound according to (42)

Influence of the size of the RVE on the macroscopic moduli
Next, we study the effect of the size of the RVE on the resulting macroscopic moduli. The influence of the size of the RVE for M mat = 10 −1 mm 4 /(Ns) and M incl = 10 −4 mm 4 /(Ns) on the effective response is shown in Fig. 10 We observe that as the size of the RVE increases the macroscopic chemical parameter ∂ 2 s s π as well as the macroscopic mobility parameter 1/(∂ 2 (35)), increase monotonically although the microscopic material parameters are not changed. However, for l micro → 0, the macroscopic moduli saturate. As becomes visible in Fig. 10a, the macroscopic chemical parameter approaches the average value of the microscopic chemical parameters, i.e., ∂ 2 ss ψ := 1 B 0 | B 0 ∂ 2 ss ψ dV = 10 N/mm 2 . Analogously, in Fig. 10b, the macroscopic mobility parameter for l micro → 0 is bounded by the classical Voigt and Reuss bounds (Zohdi and Wriggers [43]) To further study the influence of the size of the RVE depending on the mobility parameter of the microstructure, we provide three-dimensional plots for the different mobility parameters of the matrix M mat = {0.001, 0.1, 1.0} mm 4 /(Ns) in Fig. 11.
In Fig. 11a-c, we observe that the macroscopic chemical parameter is highly dependent on the RVE size l micro as well as the mobility parameter of the matrix material M mat . As M mat increases the macroscopic chemical parameter decreases. However, the influence of the mobility parameter of the inclusions is lower than the one of the matrix material. In Fig. 11c, we see that as M incl → 0 the effective chemical parameter increases in the given range. In all instances of M mat , M incl and l micro , the macroscopic chemical parameter approaches 10 N/mm 2 when the RVE size approaches zero. In Fig. 11d-f, we observe that the macroscopic mobility parameter is highly dependent on the size of the RVE and the mobility parameter of the matrix. We also see that for lower values of matrix mobility parameter, the influence of the mobility parameter of the inclusions becomes more pronounced. Additionally, the RVE size becomes more relevant for increasing magnitudes of M mat .
Finally, in Fig. 11g-i we illustrate regions of macroscopic mobility parameter that are bounded from above by the classical Voigt bound (42). It is evident that low magnitudes of matrix mobility parameter lead to larger regions bounded by the Voigt bound. The size of the RVE has a similar influence on this region.

Coupled diffusion-deformation processes in periodic hydrogels
We now present a numerical example that demonstrates the coupled chemo-mechanical response of a two-phase periodic hydrogel microstructure. To be specific, we analyze a pattern transformation due to diffusion-induced swelling phenomena (we refer to Zhu et al. [42] for experimental evidence). In order to model such phenomenon, we need to consider appropriate constitutive equations that take into account the coupled response of hydrogels. The subsequent sections are thus first devoted to the discussion of the used microscopic constitutive functions. Thereafter we analyze the behavior of a periodic hydrogel in consideration of typical geometrical and physical properties under given macroscopic loading conditions.

Microscopic constitutive functions
We assume that the mechanical contribution of the energy-storage function in (9) can be described by a neo-Hookean function. The chemical energy function is considered to be of the Flory-Rehner-type. These functions are given aŝ where α is a mixing modulus and χ is a dimensionless Flory-Huggins interaction parameter. The coupling of the mechanical and the chemical response of the material is modeled with the function where is a penalty parameter enforcing the volumetric constraint J = 1 + s, see Böger et al. [4]. The dissipation-potential function is a convex function of the solvent flux φ(H ; F n , s n ) = 1 2Ms n C n : (H ⊗ H ), where C is the right Cauchy-Green tensor and M is a mobility parameter (Böger et al. [4]). We note that the swelling-volume concentration s ≥ 0 is a non-negative quantity and s = 0 corresponds to a dry hydrogel polymer network. However, due to the present singularity of (43) 2 at the dry state we employ a stress-free pre-swollen state as the reference configuration for the Flory-Rehner energy function. We refer to Hong et al. [19,20], Böger et al. [4] for detailed discussions. Consequently, we write the final forms of the energy-storage function where J 0 is a Jacobian that characterizes the volume change due to uniform pre-swelling. The initial solvent-volume concentration is determined from the assumption of a stressfree reference configuration as The macroscopic solvent concentration s 0 at the pre-swollen state is defined as the average of the microscopic solvent-volume concentration s 0 . Penalty parameter, N/mm 2 10 10 −3 6.

Description of the problem
To analyze the effective chemo-mechanical behavior of a hydrogel we consider a twophase square-shaped RVE with size l micro = 2 mm build out of four periodic unit cells (the unit cell is depicted in Fig. 5). The volume fraction of the inclusions is π/16. The inclusions are softer than the matrix and have a higher mobility parameter. The material parameters are listed in Table 1.
The RVE is loaded with the macroscopic fields

Swelling-induced pattern transformation of a periodic hydrogel
As a numerical example, we analyze a swelling-induced pattern transformation of a twophase periodic hydrogel. The periodic RVE is discretized by using 16, 000 Q 2 RT 0 finite elements. In Fig. 12, we illustrate contour plots of the P 11 -component of the first Piola-Kirchhoff stress at three different time instances. We observe that initially the soft inclusions having a volume fraction of π/16 shrink uniformly until a critical loading point is reached. At critical loading, the RVE undergoes a pattern transformation in the form of a buckling mode. A similar behavior has been observed in experiments by Zhu et al. [42]. When the microstructure buckles the shape of the inclusions becomes non-circular, a state which is typically referred to as the diamond-plate pattern [42]. Such an instability mode has also been observed by Bertoldi et al. [2] for periodic elastic structures with voids or soft inclusions under purely compressive loading. To trigger the buckling mode shown in Fig. 12, we have slightly perturbed the finite-element mesh in the given directions.

Summary
We provided a computational homogenization framework for diffusion-deformation processes within a variationally consistent minimization-based setting that takes into account the deformation map and solvent-volume flux as independent field variables. We have discussed the theoretical aspects as well as the finite-element implementation of the formulation. The latter was realized by means of a conforming Raviart- Thomas-type discretization. By doing so, we were able to compute the macroscopic dual fields and incremental moduli in the numerical setting. Consistent with previous works from the literature, we found that the formulation yields a size-dependent macroscopic response. We confirmed this effect in a number of numerical examples by consideration of different RVEs with different mobility parameters and sizes. We further confirmed that the macroscopic mobility properties are bounded by the classical Reuss and Voigt bounds when the RVE size approaches zero. We also presented a numerical example showing a swelling-induced buckling mode of a soft periodic hydrogel.