{'fraction_non_alphanumeric': 0.03214935273349738, 'fraction_numerical': 0.009414785492671445, 'mean_word_length': 4.529429162969536, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Total Lagrangian Finite Element Formulation of the Flory-Rehner Free Energy Function

We address the total Lagrangian finite element implementation of the Flory-Rehner free-energy function in the framework of a hyperelastic material model. We explicitly give all the equations required to implement this material model in an implicit nonlinear finite element analysis, particularly, we show how to derive the so-called algorithmic or consistent linearized tangent modulus in the Lagrangian description. Some analytical and numerical results for different boundary-value problems are presented to validate the implementation.


INTRODUCTION
Flory-Rehner gels are a kind of soft matter that swell/shrink spontaneously due to a change of the chemical potential in the environment [1]. These gels are composed of a three dimensional polymeric network and small solvent molecules that can migrate in and out of the polymeric chains. The bonds between solvent and polymeric molecules can be thought as physical bonds on the order of Van der Waals forces [2]. Flory-Rehner gels are highly compressible elastic solids with a non linear material behavior that undergoes large deformations when they are in contact with a solvent at a different chemical potential. Therefore, due to these nonlinearities, it is not an easy task to develop simple theories to analyze them for general cases. Analysis of swelling of inhomogeous gels has been an active topic of study that dates back, at least, since the pioneering work of Tanaka and Fillmore [3] in the late 70's, when they tried to estimate, analytically, the characteristic time of swelling. Of particular interest is the mathematical work of Onuki [4,5].
More recently, Dolbow et al. [6,7] considered a continuum model for chemically induced volume phase transitions in hydrogels that allows for a sharp interface separating swelled and collapsed phases for the underlying polymer network. They used a hybrid eXtended-Finite Element/Level-Set method for obtaining approximate solutions to the governing equations for a two-dimensional model and get an estimate of the response time for swelling. Most notable is the recent work done by the group lead by Zhigang Suo [8,9,10], they presented a simple theory of coupled molecular migration and large deformation within the framework of non equilibrium thermodynamics where all the fields are given as a function of the material coordinates. They implemented a hyperelastic material model using the user supplied routine UHYPER in the commercial package ABAQUS. Kang and Huang [11,12] have reformulated and implemented the previous theory using the user supplied routine UMAT of ABAQUS to analyze surface instabilities of confined gels layers. The intent of this paper is to address the continuum mechanics theory and finite element implementation of a Flory-Rehner gel in a total Lagrangian framework. Although the theory presented in this paper is not new, still it has not been fully studied and to the knowledge of the authors, there is not a reference where explicitly a linearization of the internal nodal forces have been done in a total Lagrangian formulation for a Flory-Rehner gel. However, in [11] a tangent modulus tensor for the same material studied here was explicitly presented for current configuration. Someone could argue that a series of pull-back operation [13] can be applied to the tangent modulus in the deformed representation such that we could transform it into the total Lagrangian formulation. Nevertheless, we have decided to present the entire structure of the total Lagrangian implementation. Following same lines as in [8,14], we presented a purely mechanical theory of coupled species transport and large elastic deformation in a continuum mechanics approach.

GENERAL SETUP
Let introduce a fixed right-handed Cartesian coordinate system as presented in figure 1, with orthogonal unit base vectors, E i , e i , i=1,2,3 for the reference and current configuration, respectively. Select  0 as the volume occupy by the body in the reference state and  as the volume occupy by the body in the current state. We will assume that there exist a one-to-one and smooth mapping, X= −1  x , t , between the reference and current configuration for all time, t≥0 . The tensor field F=∇  is referred as the deformation gradient and should satisfy the requirement that J =det F0 for all t . The latter guarantee that the inverse motion, X= −1  x , t exist and it is unique. The boundary of the body is divided into three parts as follows:  h ∪ t ∪ u =∅ and with the restriction that = h ∪ t ∪ u . In the following subsections we will review the balance principles.

Molecular incompressibility for gels
The gel is a continuum body that could have a dramatic change in volume, therefore it is considered a compressible solid. If we look independently at the long chain network of the polymer and the solvent molecules that constitute the gel we will see that they are incompressible. Then, taking into account the incompressibility of the polymer and solvent molecules, and according to [14], the molecular incompressibility is defined in equation (1).
Where vC is the volume of the small molecules in the gel divided by the volume of the dry network. In figure 2 we present a pictorial definition of the hydrogel. Note that, although, the system is incompressible, the hydrogel is compressible, because the solvent particles could move in and out of the network. The kinetics of the swelling gels is presented in a pictorial form in figure 3. In the theory presented in this section, we have studied the state when the gel is in equilibrium with the solvent. In this steady state, there is no change in volume, unless a perturbation in the system is applied.

Mathematical derivation of the molecular incompressibility
We will assume that there is no void space or discontinuities in the gel. Therefore, the differential volume of the gel is given by equation (2).
Expressing the volume of the gel and the network with respect to the reference state (Lagrangian point of view) and noticing that V 0 refers to the volume at time t=0 and V refers to the volume at time t≥0 , we get Also, note that dV network is the volume of the small molecules (solvent) in the gel divided by the volume of the dry network. We will approximate the latter, using vC . Where v is the volume of the solvent molecules and C is a scalar field that describes the distribution of the solvent molecules in the gel. With these assumptions, we get the molecular incompressibility condition given in equation (1).

FLORY-REHNER FREE ENERGY FUNCTION
For the purpose of this analysis we will work with the Flory-Rehner free energy function [9], defined by equation (4) Where N is the number of polymeric chains per reference volume, kT is the absolute temperature in units of energy,  is a dimensionless measure of the enthalpy of mixing. In equation (4) we have two independent variables, F and C . Therefore, if we want to use equation (4) in a variational principle, then we will need to use a two-field variational principle [13]. Other way is to use a Lagrange transform and the molecular incompressibility to transform equation (4) into a function of the deformation gradient only. When the gel equilibrates with the solvent and the mechanical load, the chemical potential of the solvent molecules, μ, is homogeneous in the external solvent and inside the gel, as is presented in figure 4. In thermodynamics, the chemical potential is defined by equation (5), keeping entropy, volume and number of particles constant.
Using equation (5), the Legendre transform and the molecular incompressibilty, we get equation (6). Later we will replace C using the molecular incompressibility condition.
A combination of equation (4), equation (6) and equation (1) gives the desired free energy function.
Where I =F :F is the first invariant of the right Cauchy-Green deformation tensor. Seems like there is an inconsistency in equation (7) when J =1 , that is, solvent free, but it is not true because equation (7) is bounded, in other words, lim F=I  W  F , =0 .

HYPERELASTIC CONSTITUTIVE EQUATION FOR GELS
We will consider the hydrogel as a hyperelastic material, therefore the constitutive equation is given by equation (8), where P is known as the first Piola-Kirchhoff stress tensor [13] or nominal stress.
Next step is to replace equation (7) into equation (8) and take the derivative. In the following, we are going to use equation (9) ∂ For simplicity we will break the derivative of the free energy function in three parts.
Replacing part 1, part 2 and part 3 into equation (8), we get the constitutive equation for this particular problem.
Note that the first Piola-Kirchhoff stress tensor is singular when the network is solvent free.
Therefore, we can not use equation (10) in the reference state the dry network. In the following section is proposed to use a free-swelling as a reference state.

Free swelling as a reference state
To avoid the singularity in the free-energy function, we choose as our reference state the free swelling state [9]. We will need to specify equation (10) in terms of this state. To do it, we will use a multiplicative product of the deformation gradient, as in figure 5.
With reference to figure 5 and with the definition of the deformation gradient we get where in the last step we have replaced d x ' by F 0 d X and d x ' is a differential vector in the free swelling state, d x is a differential vector in the swelling state and d X is a differential vector in the reference state (dry network). Therefore, the multiplicative product of the deformation gradient is given by equation (11) F=F ' F 0 (11) The motion between the dry network state and the free swelling state is given by equation (12) x '= 0 X (12) where  0 is the free swelling stretch. Using equation (12) and the definition of deformation gradient, we get the deformation gradient between the free swelling state and the dry network state, equation (13), where I is the identity matrix.

Free swelling stretch
For the evaluation of the free swelling stretch, we will use equation (10) and the fact that in the free swelling state P=0 , where 0 is the zero matrix and F=F 0 .
with the above in mind, we get

It reduce to one equation with one unknowns
Nv Equation (14) is a nonlinear equation and could be easily solve using a Newton-Rahpson algorithm. Note that we have used the property: log u /v =−log v /u .

Free energy as a function of the free swelling state
We desire to reformulate equation (7) as a function of the free swelling state, that is, F' and  . Before proceeds, it is important to say that the free energy function is defined per unit reference volume. To developed the relation we will use a dimensional analysis. In the following, U is the strain energy in the swelling state, V 0 is the reference volume (dry network) and V ' is the volume in the free swelling state. Then, Note that the strain energy should be the same for both cases, the difference in the free energy function is due to the choice of the reference state. But, we can relate V 0 with V ' by mean of the Jacobian. Also, the Jacobian for the free swelling state is constant.

Free energy
The free energy as a function of the free swelling state is obtained using equation (7), equation (15) and the expression obtained in previous sections.
In the free swelling state, J '=1 , but it does not cause problem, unless,  0 =1 , that is, the gel is incompressible, therefore, violating our assumption of compressibility.

Constitutive equation in terms of invariant and elasticity tensor
In this section is introduced the most general form of a stress relation and elasticity tensor in terms of the strain invariant for the Flory-Rehner free energy function as is applied to hyperelastic material models. The most general expression for the evaluation of stresses using equation (16) is where S is the second Piola-Kirchhoff stress tensor I 1 and I 3 are the invariants of C and I is the identity matrix. For the description of isotropic hyperelastic materials at finite strain considering equation (16) in the coupled form, with invariant, I 1 and , the I 3 most general elasticity tensor is

Notes about derivatives of the Flory-Rehner free energy function respect to invariants of Green
Lagrange strain tensor.
All the equation on this section assumed that the Flory-Rehner free-energy function is a function of the first and thirds invariants of C . However, equation (16) is expressed as a function of I 1 and J =det F=I 3 1/ 2 . Therefore, to differentiate W respect to I 3 , we should use chain rule.

VARIATIONAL FORM OF THE EQUATION OF MOTION AND FINITE ELEMENT EQUATIONS
The weak form of the boundary-value problem is also known as the principle of virtual work in mechanics of materials. In this method we want to express the boundary-value problem in a variational form (integral form) instead of the original differential form. Therefore, the requirements of the continuity of the trial solution is lower than in the strong form, it suggest the name weak form. For a more detail explanation and a proof that both formulations are equivalent, see [15].
A formal statement of the weak form of the boundary-value goes as follows: Given b :   ℜ n ,  u :  u  ℜ n and  t :  t  ℜ n , find u∈V such that for all ∈W , where  is a second order tensor defined in terms of u and V is the space for all the test functions (solution space) and W is the space for all the trial functions (variation), such that the following is true, Figure 6 is a pictorial representation of the trial function. Note that  is arbitrary and the only condition is that must be zero when evaluated at the boundary where the displacement has been specified. We could look upon  as the virtual displacement field  u , defined on the current configuration. Taking this into account we want to recast equation (28) as a function of the variation of the displacement and the Cauchy stress tensor (  ). Considering the symmetry of the Cauchy stress tensor and Where  e is the variation of the Euler-Almansi strain tensor [13] defined in the current state.
Substituting equation (31) into equation (28) we get, This is because, in general, we do not know the domain of integration at the current time and the material model used in this paper is nonlinear.

Total Lagrangian formulation
In this paper we used the Total Lagrangian (TL) formulation to deal with the unknown integration domain. A TL formulation is based on the assumption that there exist an inverse mapping between the current and a reference state, that allows us, to refers everything to a known reference state at the beginning of the motion. To express the principle of virtual work in the reference state we will use the second Piola-Kirchhoff stress tensor and the Green-Lagrange strain tensor ( E ).
Equation (33) is the basic expression of the TL formulation used in this paper. Note that the integration is respect to the known reference configuration at time t =0 . It is equivalent as what we do in finite element analysis with the use of a parent element.

Integration in a finite element solution based on Total Lagrangian approach
Integrals on the reference configuration are related to integrals over the parent element by where the integration is carried out on the parent element domain (  ) and J  is the Jacobian between the reference configuration and the parent element, as shown in figure 7 for a two dimensional quadrilateral element. Note that, J  does not depends o the current configuration and therefore just need to be evaluated at the beginning of the calculations.

Finite element approximation
To solve approximately the equations on previous section, we will use an isoparametric finite element discretization. The problem domain is defined by a mesh. Geometry and independent variables are interpolated using compact piecewise linear functions. In general, the geometry and variables will be interpolated using equation (35) for each element in the mesh.
where  corresponds to the parent element coordinate, N  are the standard finite element shape functions [19] and  x and  u are the vectors of coordinates nodal values and nodal variables, respectively. Then, using the previous interpolations and Voigt notation [16] we can express equation (33) in matrix form and we will call it, the internal nodal force vector.
where the B 0I express the variation of the Green strain [16,17] and is given as follows for a three dimensional problem, and the second Piola-Kirchhoff stress "vector" is given by, Applying a linearization process to equation (36) and using finite element interpolation, we arrived to the following system of equations, Where K mat and K geo are the well known material and geometric stiffness matrix, respectively and R is the vector of external point forces. This matrices are defined by, where B 0I =∂ N I /∂ X , D is the algorithmic tangent modulus ( equation (19) ) in Voigt notation and I is the identity matrix.

Free swelling
Free swelling refers to a gel configuration in which there are no constraints for the motion. To simulate it into a finite element solution we removed the rigid body motion of the body, applying proper boundary conditions. In figure 8 is shown a comparison between the analytical and finite element results. The gel was simulated using one element as shown in figure 9.
Note also that the displacement field shown in figure 9 was measured respect to the free swelling state, because of our implementation of the Flory-Rehner free energy function. According to that, the relation used to calculate the free swelling stretch is given by equation (42) =1 where  0 is the initial free swelling stretch,  is the displacement, measure from the initial swelling state, and L is the original free swelling cube length.

A bar of a gel in contact with a solvent and subjected to a uniform axial stress
In this example we put in contact a bar of gel with a solvent and an axial force was applied. The longitudinal stretch is given by  1 and the traverse stretches are given by  2 = 3 . The analytical solution is given in equation (43) and it is used to calculate (using Newton's method)  2 for every known value  1 . The bar was simulated using one finite element as shown in figure 10. The comparison between the analytical and numerical solution is presented in figure   11, and it shows that our material model implementation is correct.

CONCLUSIONS
We have shown some preliminaries results that involves the use of a Flory-Rehner free energy