Principle of virtual power and drilling degrees of freedom for dynamic modeling of the behavior of liquid crystal elastomer films

In this work, we aim to model the reorientation process of mesogens in nematic liquid crystal elastomers within the context of dynamics. We consider a continuum model with separate mappings for the deformation of the monolithic material and the orientation of the nematic director, where the latter describes the inclination of the mesogens. We achieve the inextensibility of the nematic director through the introduction of drilling degrees of freedom. We combine this approach with the application of the principle of virtual power and a mixed finite element formulation, in order to formulate distinct momentum and angular momentum balance laws for the two separate mappings. Furthermore, we include in our continuum model a volume load and a surface load associated only with the orientation mapping. We show in the presented three numerical examples that our formulation enables the fulfillment of all momentum and angular momentum balance laws.


Introduction
Liquid crystal elastomers (LCEs) are under investigation since a few decades due to their capability to respond with large deformations to external stimuli, such as heating and the application of electric fields. Because of this feature, a LCE can perform motion actuation and is thus evaluated as a suitable material for the development of artificial muscles. The behavior of LCEs is related to their instrinsic structure, as the characteristics of rubber are integrated with the orientational properties of liquid crystals [1][2][3]. Based on the grade of ordering of molecules, LCEs can be categorized into three classes: smectic, nematic and cholesteric [4]. The isotropic state of LCEs is achieved by heating above the transition temperature. In this work, we focus on nematic LCEs, in which rodlike parts of the molecule called mesogens are linked to the polymer backbone. In case of a monodomain microstructure, the mesogenic rods are aligned toward a unique nematic director [5] and LCEs are thus manufactured as thin films in order to keep the alignment of the nematic director through the whole thickness of the specimen [6]. Due to their unique features, nematic LCEs are currently the subject of intensive research and several works deal with the modeling of their behavior. In the works of Leslie [7,8], a continuum theory for liquid crystals was formulated from balance laws by including volume and surface Communicated by Andreas Öchsner.
forces related to the director in the former, and by introducing a viscous dissipation in the latter. As one of the main features of LCEs, the modeling of soft elastic response was the subject of several works, e.g., [9][10][11][12]. More recently, Agostiniani and DeSimone [13] developed a model for thin films at large strains by considering different nematic textures. Another crucial issue concerns the modeling of LCEs subjected to external stimuli, along with the experimental validation, e.g., [14] which leads to the model for the prediction of the actuation [15][16][17]. In this work, we present a variational-based approach for model the dynamic behavior of nematic LCEs. We have introduced two separate mappings for describing the bulk material and the nematic directors [18], whose reorientation is modeled through drilling degrees of freedom [19]. A list of the used symbols is reported in Table 1 2 The finite element formulation

Continuum model
We take into account the LCE material as a continuum (see Fig. 1); therefore, we consider the reference configuration B 0 at initial time t = 0 with the corresponding nematic director n 0 , which is a function of the material point X. We introduce two distinct mappings: the deformation mapping ϕ (X, t) : B 0 × T → B t and the orientation mapping χ (X, t) : B 0 × T → R n dim give back the spatial point x and the nematic director n t in the current configuration B t , respectively. Furthermore, the two mappings fulfill the conditions related to the reference configuration, i.e., ϕ (X, 0) = X and χ (X, 0) = n 0 . In order to deal with dynamic problems, two further variables must be included: the velocity vector v (X, t) :=φ (X, t) =ẋ as time derivative of the deformation mapping and the linear momentum vector p := ρ 0 v. Since we aim to formulate a distinct momentum balance for each mapping function [18], an orientational velocity vector v χ (X, t) :=χ (X, t) =ṅ t as time derivative of the orientation mapping and an orientational momentum vector p χ := ρ 0 l 2 χ v χ have been introduced, where l χ is the radius of gyration. Anderson et al. [18] defined the term ρ 0 l 2 χ as referential orientational inertia density and it keeps consistency with the units of the parameter, since the nematic director is a unit vector [18]. For a detailed description of the radius of gyration, we refer to the book of Warner and Terentjev [1]. Another approach for taking into account the inertia of the LCE continuum has been presented by Groß et al. [20].
From the deformation mapping ϕ, we introduce the deformation gradient F := Grad [ϕ], the right Cauchy-Green strain tensor C := F T F and the Jacobian determinant J := det[F], which are needed for formulating the Neo-Hookean material model. Anderson et al. [18] divided the strain energy into three parts: the elastic energy density, which is a function of C and of the nematic director in the reference configuration n 0 , the nematic energy density as a function of the orientation gradient G := Grad [χ] and the interactive energy density, depending on the invariant F T n t . In this work, we let out the terms depending on the orientation gradient, consistently with the strain energy density formulations which were described in other works, e.g., [2,21]. Therefore, the strain energy is described as the sum of elastic energy density ψ e and interactive energy density ψ i with the following parameters μ is the shear modulus and r is the ratio between lengths of steps parallel and perpendicular ⊥ to the nematic director [1,18]. The terms with prefactors c 3 , c 9 and c 10 in Eqs. (2)(3) couple the orientation mapping χ with the deformation mapping ϕ, i.e., the nematic directors both in the reference n 0 and in the current configuration n t are coupled with the deformation gradient F and both right and left Cauchy-Green strain tensors C and B, where the latter can be easily obtained from the term of the interactive strain energy density  with prefactor c 9 just by using few simple rules of tensor algebra. In contrast with the work of Anderson et al. [18], we assume that the ratio r is the same for both reference and current configurations, i.e., This assumption is motivated by the work of Warner and Terentjev [1]: they expressed the ratio r as a function of the orientational order parameter as well as the square of the ratio between the radii of gyration of the ellipsoid, which represents the chain distribution. Since our model aims to describe only isothermal loading cases, the changes in the orientational order are of minor importance [10]. Hence, it can be deduced that the ellipsoid keeps its shape with the same ratio between the radii of gyration parallel and perpendicular to the nematic director. Furthermore, the sum of the first two terms of the elastic strain energy density with the prefactor c 1 in Eq. (2) with the further terms having c 3 , c 9 and c 10 as prefactors corresponds to the well-known trace formula of Warner and Terentjev [1], if the assumptions of Eq. (5) are made and some rules of tensor algebra are applied.
In the scientific literature, the LCE material has been modeled either as incompressible, e.g., [4], or as nearly incompressible, e.g., [16]. We apply the Neo-Hookean model for compressible elastomers [22] and the nearly incompressibility is fulfilled by using suitable values for the shear modulus μ and the Poisson ratio ν, which define the second Lamé constant λ. 1 The rotation of the nematic director, i.e., its reorientation, is described by means of the drilling degrees of freedom, which are depicted in Fig. 2 for a single linear hexahedral element. The time derivative of the continuum rotationα is defined as the axial vector of the current velocity gradienṫ in order to model the kinematic reorientation process of the nematic director, where is the third-order Levi-Civita tensor andḞ is the time derivative of the deformation gradient. In the finite element formulation, the introduction of the time derivative of the continuum rotationα as a variable offers the advantages related to the potential applications of related Dirichlet boundary conditions; nonetheless, no boundary conditions forα have been imposed in this work. Moreover,α is needed for describing the time derivative of the orientation mapping as follows:χ =α × χ.
Equation (7) is the widely known formula for describing rigid-body dynamic rotations [23], which we impose as a constraint in our variational formulation (see Eq. (21)) for the purpose of forcing the orientation mapping and thus the nematic director to rotate without any stretching. The spin tensor of such rigid rotation is expressed of the skew-symmetric Kirchhoff stress tensor τ skw . According to Eq. (8), the Lagrange multiplier w τ contains the third-order Levi-Civita tensor, which must be single-contracted with the third-order Levi-Civita tensor in the right-hand side of Eq. (6), and with the time rate of the continuum rotationα in the left-hand side. Hence, the spin tensor, i.e., the skew-symmetric part of the spatial velocity gradientḞ F −1 , is constrained to be equal to the skew-symmetric second-order tensor associated with the time rate of the continuum rotation − ·α by means of the Lagrange multiplier 2 w τ . The Lagrange multiplier for Eq. (7) is the reorientation stress vector τ χ associated with the nematic director.

Variational-based weak formulation
The weak forms are obtained starting form the principle of virtual power [19] t n+1 where fields with continuous time evolution, i.e., the time derivativesφ,ṗ,v,χ,ṗ χ ,v χ ,α and discontinuous time evolution w τ , τ χ , R appear. Equation (9) is valid for any time step T n := t n , t n+1 and the functional H is given by the sum of the internal energy int , both kinetic energies T and T n , as well as both external energies ext and ext n , associated with the deformation mapping ϕ and the orientation mapping χ Hence, the fulfillment of Eq. (9) is associated with The distinction between kinetic energies T and T n and external energies ext and ext n for the two mappings is made according to Anderson et al. [18]. We start by presenting each functional and calculating the corresponding variation. The functional of the deformational kinetic power depending on the fieldsφ,v andṗ readṡ The variations with respect to the fieldsφ,v andṗ are calculated by means of the Gateaux derivative The functional of the orientational kinetic power, which depends on the fieldsχ,v χ andṗ χ is given bẏ and the corresponding virtual orientational kinetic power Consistently, again with the work of Anderson et al. [18], we define as deformational and orientational the energies that are associated with the deformation mapping ϕ and the orientation mapping χ, respectively. The deformational and orientational kinetic powers in Eqs. (12) and (14) are obtained by summing up the integral of the time derivative of the kinetic energy density on B 0 with the integrals of the following constraints: v =φ andv =φ for the deformational kinetic power and v χ =χ andv χ =χ for the orientational kinetic power. Constraints are satisfied by means of the Lagrange multiplier method withṗ, p,ṗ χ and p χ as Lagrange multipliers. This procedure has been taken from Groß et al. [24].
In the functional of the deformational external power, the external loads such as the volume load B and the surface load T occur, which act upon the whole body B 0 and the Neumann boundary ∂ T B 0 , respectively. The reaction force R is a Lagrange multiplier, which acts on the Dirichlet boundary ∂ ϕ B 0 = ∂B 0 \∂ T B 0 for enforcing the constraint related to the Dirichlet boundary conditionφ =φ, which is expressed as time derivative [24] The virtual deformational external power with respect to the variations of the time derivative of the deformation mappingφ and the reaction force R reads We introduce the functional of the external power and the virtual orientational external power in Eq. (19) associated with the time derivative of the orientation mappingχ with the relating volume load B χ and surface load W , with the latter acting on the Neumann boundary ∂ W B 0 . In this work, we do not introduce any Dirichlet boundary condition for the orientation mapping; therefore, a corresponding orientational reaction force as Lagrange multiplier is not needed.
The time derivative of the internal energy˙ int , as it appears in Eq. (10), is the internal power In the internal power, we include the strain energy densities of Eqs. (1-3), the constraints for the time rates of the continuum rotationα and the orientation mappingχ and the two stress vectors w τ and τ χ as Lagrange multipliers. The internal power reads The virtual internal power is obtained by calculating the variations with respect to the five fieldsφ,χ,α, τ χ and w τ , considering that the time derivative of the deformation mappingφ occurs in the internal power as time derivative of the deformation gradientḞ : By summarizing all variations from the virtual powers, we obtain five weak forms which are solved by a monolithic approach. The deformation mapping is obtained from the weak balance of the linear momentum whereas the orientation mapping χ is given by the weak balance of the orientational momentum The weak balance of the reorientation stress expresses the reorientation stress τ χ as a function of the orientation mapping χ and the stress vector w τ of Eq. (8) The weak continuum rotation and the weak balance of the orientation rate correspond to the two constraints of Eqs. (6) and (7), respectively.

Balance laws
The weak forms of Eqs. (23-27) fulfill balance laws, which are obtained by using particular types of test functions. By substituting the test function δ * φ = c = const. in the weak balance of the linear momentum, we obtain where the derivative of the deformation gradient with respect to the deformation mapping occurring in Eq. (23) can be expressed as Since the gradient of the constant term c is equal to zero, we determine the balance of the linear momentum by canceling the constant c from the other integrals Analogously, we formulate the balance law of the orientational momentum by inserting the test function δ * χ = c = const. in Eq. (24) By summing up both balances from Eqs. (30-31), we obtain the following total momentum balance depending on the external loads, the reorientation stress vector τ χ and the partial derivative of the strain energy with respect to the orientation mapping In order to determine the balance law for the moment of linear momentum, we insert the test functions δ * φ = c × ϕ and δ * α = c in the weak forms of Eqs. (23) and (25), respectively, with c as constant term The balance law of the moment of orientational momentum results from employing the test function The space-time integration of the cross-product between the orientation mapping and the reorientation stress vector appears in both balance laws with opposite signs. Furthermore, based on the strain energy functions in Eqs. (3)(4), the following identity is valid: which is demonstrated in the Appendix. Consequently, by summing up Eqs. (33-34) we obtain the total moment of momentum, which is depending only on the external loads A further balance related to the orientation mapping results by inserting the test function δ * τ χ = χ in the weak form of Eq. (27) since the cross-product of a vector with itself is equal to zero, i.e.,α · (χ × χ) = 0. It leads to a further aspect of the reorientation process for the nematic director: given that Eq. (37) is valid, the two vectors χ andχ are perpendicular to each other and the nematic director is thus forced to rotate without undergoing any elongation or shortening. The inextensibility of the nematic director is linked to its length, which must be equal to one; therefore, we introduce the following weak balance of the reorientation function:

Space and time discretization
In order to implement our model, we discretize the five variablesφ,α,χ, w τ and τ χ in the weak forms of Eqs.  According to the continuous Galerkin method, all test functions δ * φ , δ * α , δ * χ , δ * w τ and δ * τ χ must be discretized with Lagrange polynomialsM I (ξ J ) of degree k − 1 as in Eqs. (42-43). A complete and concise description of the continuous Galerkin method can be found in the work of Erler and Groß [25]. Rules for time discretization must be of course strictly followed for the implementation of residuals and tangents, and the corresponding discretized equations are not reported in this work for the sake of brevity, but our approach has been partially described again in the work of Groß et al. [24].

Numerical results
We describe the results of the simulations on a LCE film. We consider the same geometry of LCE film for all simulations as a film of size 0.0125 × 0.075 × 0.0003 m [20] with 1500 hexahedral elements (H1). A quadrilateral element with drilling degrees of freedom could be considered as an alternative to the hexahedral element (see, e.g., [26]), but we refrain from using plane finite elements, since we would like to include three-dimensional motion in a future work. Moreover, we employ for all simulations a finite element method and the continuous Galerkin method with a degree k = 2 for the Lagrange polynomial of the time shape functions, which have been implemented in our in-house finite element software. For a detailed description on the time discretization for fields having either a continuous time evolution or a discontinuous time evolution, we refer to the work of Groß et al. [19]. The time step size is 0.0005 s from 0 up to 0.1 s, which corresponds to 200 time steps. We have considered the infinity norm of the global residual as convergence criterion, i.e., the highest value in the residuals of the weak forms from Eqs. (23-27) must be lower than the prescribed tolerance of 10 −8 . We have taken into account the data reported from de Luca et al. [2] for the shear modulus μ = 3.43 × 10 4 N/m 2 and the ratio between polymer step lengths r = 1.88. Although these data are referred to a temperature of 60 • C, in this work the temperature does not influence in any other way the mechanical response of the LCE. We have taken the radius of gyration l χ = 28 × 10 −10 m from the book of Warner and Terentjev [1] in the direction parallel to the nematic director for the temperature of 60 • C. Since the radius of gyration perpendicular to the nematic director is equal to 23 × 10 −10 m, the difference between the two radii of gyration in both directions is negligible. The density is ρ 0 = 1760 kg/m 3 [20] and the Poisson ratio is set to ν = 0.49 in order to lend a quasi-incompressible behavior to the material.

LCE film under initial rotation
In the first simulation, the LCE film is subjected to an initial rotation ω 0 = 25 rad/s about an axis parallel to the z-axis (direction +e z ), passing through the center of mass of the considered geometry, as shown in Fig. 3a, b. In the reference configuration, the nematic director is oriented in the direction e y . Neither Dirichlet nor Neumann boundary conditions are applied in this case; therefore, all momenta and moments of momenta (both linear and orientational) must be conserved. In Fig. 3c, d, the distribution of the magnitude of the reorientation stress vector on the LCE film surface is shown for the time step in the middle, i.e., t = 0.05 s and the final time step, i.e., t = 0.1 s, along with only one nematic director every 20 nodes, whereas a single nematic director is depicted every 40 nodes in figures representing the reference configuration for the sake of understandability. As Lagrange multiplier of the constraint in Eq. (7), the reorientation stress vector corresponds to the stress which makes the nematic director rotate about the direction describing the time rate of the continuum rotation and thus impedes the nematic director from stretching. The nematic directors rotate uniformly keeping their orientation with respect to the specimen. The magnitude of the reorientation stress vector increases from both ends of the LCE film up to the end of the simulation, when a spot-like pattern for the magnitude of the reorientation stress arises and the highest magnitude is finally reached over the whole specimen. Figure 4a refers to the total momentum in the left y-axis and the total moment of momentum in the right y-axis over time. The total momentum in the current time step L T (t n+1 ) = L (t n+1 )+ L n (t n+1 ) is given by the balance in Eq. (32) as a sum between momenta in the previous time step L (t n ) + L n (t n ) and the terms in the right-hand side of Eq. (32). Analogously, the total moment of momentum in the current time step J T (t n+1 ) = J (t n+1 ) + J n (t n+1 ) is given by the balance in Eq. (36) as sum between the total moment of momentum in the previous time step J (t n ) + J n (t n ) and the terms in the right-hand side of Eq. (36). Figure 4a shows that the total momentum and total moment of momentum are conserved, as expected. In Fig. 4b, the orientational momentum and the moment of orientational momentum are shown, given by the balances in Eqs. (31) and (34), respectively. The orders of magnitude of the orientational momentum and its moment are much lower by comparison with the total moment and total angular momentum, which includes also the linear momentum and its moment. This is due to the orientational momentum p χ being determined by the square of the radius of gyration, as shown in Eq. (15). The orientational momentum and the moment of orientational momentum are conserved as well, but an abnormal trend can be observed by approaching the end of the numerical test. The sum of the orientation stress vector with the partial derivative of the strain energy with respect to the orientation mapping mainly affects the trends of the orientational momentum and its moment. The trend might be linked to energy consistency issues, which have been not addressed in this work. Figure 4c shows the error of the balances for the linear momentum (Eq. (30)) and its moment (Eq. (33)), whereas Fig. 4d shows the error of the balances for the orientational momentum (Eq. (31)) and its moment (Eq. (34)). All errors are lower than the prescribed tolerance and all momentum balances are thus fulfilled. In Fig. 4d, each error exhibits a trend which is symmetrical to the associated orientational momentum and its moment with respect to the x-axis. Since the error is given by the left-hand side of Eq. (34) subtracted by its right-hand side and then divided by the prescribed tolerance, Fig. 4 Total momentum and total moment of momentum (a), orientational momentum and moment of orientational momentum (b), error of the linear momentum balance and error of the moment of linear momentum balance (c), error of the orientational momentum balance and error of the moment of orientational momentum balance (d) over time for a LCE film subject to the prescribed initial rotation. In all plots, the left y-axis is referred to the momenta and the right y-axis to the moments of momentum the terms have opposite signs. The orders of magnitude for the orientational momentum and its moment are much lower than those of the correlated errors because the latter are divided by the tolerance of 10 −8 .

LCE film under volume load
The second representative simulation consists in the LCE film specimen being clamped in its left side as Dirichlet boundary condition (See Fig. 5a). The Neumann boundary condition is a constant orientational volume load over time equal to 0.1 m 2 /s 2 in direction −e y . In the reference configuration, the nematic director is oriented in direction e x . In Fig.5c, d the distribution of the magnitude of the reorientation stress vector and the nematic directors are shown for the same time steps as in the previous example. The rotation of the nematic director is more pronounced by moving away from the clamping in the left end of the LCE film. The greatest rotations of the nematic director occur noticeably at the right end of the specimen. The rotation of the nematic directors is accompanied by the deformation of the embedding bulk elastomer. Figure 6a displays Total momentum and total moment of momentum (a), orientational momentum and moment of orientational momentum (b), error of the linear momentum balance and error of the moment of linear momentum balance (c), error of the orientational momentum balance and error of the moment of orientational momentum balance (d) over time for a LCE film undergoing the prescribed volume load. In all plots, the left y-axis is referred to the momenta and the right y-axis to the moments of momentum the total momentum and total angular momentum over time. A decreasing trend for the momenta in the xand y-direction is observed along with the decrease in the angular momentum in the z-direction. The decrease in the trends is due to the term associated with the reaction forces, the term L (t n ) in Eq. (32) and the term J (t n ) in Eq. (36), since no other Neumann boundary condition is active, e.g., gravity. In all other directions, the total momentum and total moment of momentum are constant. We notice an analogous response for the orientational momentum and moment of orientational momentum in the same directions, although the trend appears scattered with a significantly smaller order of magnitude. The balance laws for momenta and moments of momentum are fulfilled, as shown in Fig. 6c, d. The increasing trend for the error of the moment of linear momentum, which is noticeable in Fig. 6c, is related to the decreasing trend of the total moment of momentum in Fig. 6a.  As third example, we consider the same specimen of LCE film being clamped in its left side with an orientational surface load acting upon the very thin surface in the right end of the specimen. The orientational surface load is 3 × 10 −7 Nm −1 in direction −e y . We do not impose any other Neumann boundary condition. The nematic director in the reference configuration is oriented in direction e x in this case too. The Dirichlet boundary condition is shown in Fig. 7a, which is the same as considered for the LCE film under orientational volume load. In Fig. 7a, the blue arrows refer to the direction of the orientational surface load. Figure 7b, c depicts the reorientation stress distribution and the nematic directors for the two considered time steps. The response of the LCE film to the orientational surface load is similar to the one observed for the orientational volume load, but the rotation of the nematic director is much lower and the highest magnitude of the reorientation stress vector is detected near the narrow area subjected to the Neumann boundary condition. Accordingly, the resulting deformation of the LCE film is significantly lower as compared to the deformation of the LCE film subjected to an orientational volume load (cf. Figs. 5c, 7c). Concerning the total momentum and total moment of momentum of Fig. 8a, we observe a remarkable decrease in the y-direction of the total momentum and the z-direction of the moment of momentum. Such trends are influenced in the same way as the LCE film under the orientational volume load. In Fig. 8b, the trends of the orientational momentum and the moment of orientational momentum appear concentrated because the orientational surface load acts only in a restricted region of the specimen. Figure 8c, d shows that the balance laws for momenta and moments of momentum are fulfilled.
For all representative simulations, the balance law of the reorientation function of Eq. (38) has been proven as well. The error of the reorientation function over time is shown in Fig. 9 for the case of the LCE film under initial rotation, but all presented loading cases exhibit an identical response for this error.

Conclusion
The goal of our work was the modeling of the reorientation process of mesogens in a LCE film under dynamic loads or free motion. For this reason, we considered a continuum model with two separate mappings for describing the deformation and the nematic director, where the latter represents the unique direction of the mesogens. The introduction of two separate mappings led to distinct formulations for the energy, consistently with the work of Anderson et al. [18]. The reorientation of the nematic directors was accomplished through the drilling degrees of freedom and the related Lagrange multipliers as independent variables, for the purpose of ensuring the inextensibility of the nematic director. We simulated the response of a nematic LCE film subjected to orientational volume load and surface load [18], as well as to free rotation. The weak forms and the balance laws were determined by using the principle of virtual power, and all presented momentum balances were preserved thanks to the usage of the continuous Galerkin method. This work lays the foundations for the prediction of the motion actuation, which is the most remarkable characteristic of LCEs. Moreover, the formulation could be improved by considering further variables, which were let out in this work, such as the orientation gradient. A future work should also deal with the issue of energy consistency.

Appendices A
In order to demonstrate the identity in Eq. (35), the following rule for the cross-product between second-order tensors must be introduced [19] A × B := : AB T . (A1) In the elastic strain energy density, only the deformation gradient appears and hence only the left-hand side of Eq. (35) is involved. In this case, the cross-product between the deformation gradient and the derivative of all Neo-Hookean terms must be zero, in order to fulfill the identity in Eq. (35), as shown in Eqs. (A2)-(A7))) In Eq. (A8), the left Cauchy-Green strain tensor B = F F T appears analogously to the free energy density described by Conti et al. [21] F T n t 2 = F T n t · F T n t 1 2 2 = F F T : (n t ⊗ n t ) .
However, in Eq. (A8) the double contraction between the second-order tensors F and n t ⊗ n 0 [23] gives back a scalar, which is factored out from the cross-product. We calculate also the right-hand side of Eq. (35) for the interactive energy density by substituting the orientation mapping n t = χ (n 0 , t) as the nematic director in the current configuration χ × ∂ c 9 F T n t 2 ∂χ = 2c 9 n t × F F T n t (A10)