{'fraction_non_alphanumeric': 0.05710335707630812, 'fraction_numerical': 0.020346827758242418, 'mean_word_length': 3.88820332011165, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical Study of Liquid Crystal Elastomer Using Mixed Finite Element Method

We aimed to use finite element method to simulate the unique behaviors of liquid crystal elastomer, such as semi-soft elasticity, stripe domain instabilities etc. We started from an energy functional with the 2D Bladon-Warner-Terentjev stored energy of elastomer, the Oseen-Frank energy of liquid crystals, plus the penalty terms for the incompressibility constraint on the displacement, and the unity constraint on the director. Then we applied variational principles to get the differential equations. Next we used mixed finite element method to do the numerical simulation. The existence, uniqueness, well-posedness and convergence of the numerical methods were investigated. The semi-soft elasticity was observed, and can be related to the rotation of the directors. The stripe domain phenomenon, however, wasn't observed. This might due to the relative coarse mesh we have used.


Introduction
Nematic liquid crystal elastomer is a relatively new kind of elastic material which combines properties from both incompressible elasticity and rod-like liquid crystals. It has some unique behaviors such as the stripe-domain phenomenon and the semi-soft elasticity. These phenomena are usually observed in the "clamped-pulling" experiement, in which a piece of rectangular shape liquid crystal elastomer is clamped at two ends and pulled along the direction perpendicular to the initial alignment of the directors.
Mitchell et al. [17] did the "clamped-pulling" experiement for acrylate-based monodomain networks. They found that the directors rotate in a discontinuous way: after a critical strain, the directors switch to the direction perpendicular to the original alignment. On the other hand, Kundler and Finkelmann [13] did the same experiment for polysiloxane liquid crystal elastomer. They found that the directors rotate in a continuous manner. They also observed the interesting stripe domain phenomenon. Later, Zubarev and Finkelmann [20] re-investigated the "clampedpulling" experiment for polysiloxane liquid crystal elastomers, to investigate the role of aspect ratio of the rectangular shape monodomain in the formation of stripe domains. Figure 1 shows the process of the "clamped-pulling" experiement in the case that the aspect ratio is AR = 1. In this graph, the pulling direction is vertical, while the initial alignment of the directors is horizontal. The size of the elastomer was 5 × 5mm 2 , with thickness 50 − 300 ± 5 µm. They found that at some critical elongation factor λ = 1.1, an opaque region started to emerge in the center of the elastomer. Using X-rays, they found that this region was made of stripes of about 15µm in width. In each stripe, the directors were aligned uniformly, while across the stripes, the directors were aligned in a zig-zag way, as shown in the second picture of Figure 1. When the elastomer was pulled further, the opaque region was broken into two regions symmetric about the center (the Figure 1: The stripe domain phenomenon. Reproduced from [20] third picture of Figure 1). And these two regions moved closer to the two ends as the elastomer was pulled further. Eventually at λ = 1.3 the two stripe domain regions reached the two ends (the last picture of Figure 1) and in this last stage, most of the directors at the center had rotated 90 degrees, and were now aligned in the pulling direction (vertical direction).
Another interesting phenomenon of liquid crystal elastomers is the so-called "semi-soft elasticity". Figure 2 shows the stress-strain graph of three different kinds of liquid crystal elastomers Figure 2: The semi-softness of liquid crystal elastomer. Reproduced from [19].
in the "clamped-pulling" experiment ( [14,7]). We can see that at the beginning, the nominal stress increases linearly with the strain, just like ordinary elastic materials. However, after the strain reaches a critical point, the stress-strain curve becomes relatively "flat", which we call the "plateau region". When the strain falls in this region, relatively small stress can induce relatively large deformation, which means the elastomer is "soft". We can also see from Figure 2 that when the strain continues to increase and reaches another critical value, the stress-strain curve becomes linear again.
The stored energy of liquid crystal elastomer proposed by Bladon, Warner and Terentjev (BTW) is [2] W (F, n) = µ |F | 2 − (1 − a)|F T n| 2 − 3a 1/3 (1) where F = I + ∇u(X) is the deformation gradient of the elastomer network, while n = n(X) is a unit vector denoting the orientation of the rod-like liquid crystal polymers, and X denotes an arbitrary point in the reference configuration. We assume that the elastomer is incompressible, which is equivalent to det(F ) = 1. Also µ is the elasticity constant, and a ∈ (0, 1) is a constant measuring the significance of the interaction between the displacement u and the orientation n. Notice that in the limit a = 1, the BTW energy degenerates to the stored energy of incompressible neo-Hookean materials. The closer the value a is to 1, the smaller the interaction between the displacement u and the orientation n. Similarly in 2D, the BTW energy can be defined as The only difference is that the constant term is changed from 3a 1/3 to 2a 1/2 . Proposition 1. The BTW energy W (F, n) in equation (1) and (2) are always non-negative.
When there is no shearing and the BTW energy (1) is zero for t = 0, 1 and strictly positive for 0 < t < 1. However, when there is shearing and the BTW energy (1) can be zero for 0 ≤ t ≤ 1.
When there is no shearing and the BTW energy (2) is zero for t = 0, 1 and is strictly positive for 0 < t < 1. However, when there is shearing and the BTW energy (2) can be zero for 0 ≤ t ≤ 1.
Corollary 3. If we take then W (F ) is a non-convex function of F .
Proof. Take 0 < t < 1, and λ(t) and F ± as defined Proposition 2. Then we have Therefore W (F ) is a non-convex function of F .
Proposition 2 tells us that introducing shearing might lower the BTW energy (1) or (2). On the other hand, because of the constraint of the clamps, global shearing is not possible, therefore, local shearing is developed in a zig-zag way, and that's why stripe domains occur.
In [8], DeSimone et al. did the finite element study of the clamped pulling of liquid crystal elastomer. They started from the 3D BTW energy (1), and got rid of the variable n by taking it to be always n = n(F ), the eigenvector of the largest eigenvalue a −2/3 of the matrix F F T . So they got It turned out that the energy function (11) is no longer a convex function of F (Corollary 3). Therefore, they replaced it by its quasi-convex envelope They pointed out, that the use of W qc (F ) in the numerical computations allows one to resolve only the macroscopic length scale, with the (possibly infinitesimal) microscopic scale already accounted for in W qc (F ). It turned out that their approach was very successful, in that they successfully captured the emerging and migrating of the stripe domain regions, in exactly the same way as those observed experimentally in [20]. However, their approach is limited for several reasons: • Although they were able to tell which part of the region is a "stripe-domain region", they couldn't tell how many stripes lie in that region. Actually, their model permits infinitely many stripes in that region, which is unphysical.
• They were only able to observe the "soft elasticity", not the "semi-soft elasticity". That is, the plateau region occurs immediately upon the pulling (Figure 3), instead of after some critical strain as observed experimentally. Figure 3: The stress-strain curve from Desimone's numerical simulation. Soft elasticity instead of semi-soft elasticity was observed. Reproduced from [8].
In this work, we aim to resolve the above issues. We combine the 2D BTW model (2) and the Oseen-Frank model ( [11]), and get the following energy functional for liquid crystal elastomers: where f is the body force, while g is the traction force on the Γ 2 part of the boundary ∂Ω. We assume u ∈ H 1 (Ω), and satisfies det(I + ∇u) = 1 a.e. in Ω; while n ∈ H 1 (Ω), and satisfies |n| = 1 a.e. in Ω. Let L be the characteristic length scale of the domain Ω, we get the following non-dimensionalized energỹ whereΩ is congruent with Ω with characteristic length scale 1, andΓ 2 ⊂ ∂Ω is the image of Γ 2 . Alsob = b cL 2 ,f = f cL , andg = g c . The rest of the paper is organized as the following. In section 2, we investigated the properties of the continuous problem, such as the existence of minimizer, the equilibrium equation and its linearization, the stress-free state, the well-posedness of the linearized equation etc. In section 3, we investigated problems related to the discretization, such as the existence of minimizer, the equilibrium equation and its linearization, the well-posedness of the linearized equation, and the existence and uniqueness of the Lagrange multipliers etc. In section 4, we presented the simulation results for the "clamped-pulling" experiment, including the results of the inf-sup tests, the convergence rates, and the stress-strain curve etc. Finally, in section 5, we discussed future directions.

Proof. Let
Then we have Then for any F 1 , F 2 ∈ M 3×3 , and 0 ≤ α ≤ 1, we have where we have used Proposition 1. Therefore L is a convex function of F .
Remark It's easy to see that Lemma 4 and Lemma 5 holds for 3D, too.

Theorem 6.
Let Ω be a nonempty, bounded, open subset of R d .
• If d = 2, suppose we have u k u in W 1,s with s > 4 3 , then we have det(I + ∇u k ) → det(I + ∇u) in D (Ω).
suppose we have u k u in W 1,s with s > 3 2 , then we have adj(I + ∇u k ) ij → adj(I + ∇u) ij in D (Ω).
suppose we have u k u in W 1,s , and adj(I + ∇u k ) adj(I + ∇u) in L q (Ω; M 3 ) with s > 1, q > 1 and 1 s + 1 q < 4 3 , then we have det(I + ∇u k ) → det(I + ∇u) in D (Ω). For the proof, please see John Ball [1]. Proof. Let m be the infimum of Π on A(u 0 , n 0 ), and let (u k , n k ) ∈ A(u 0 , n 0 ) be a minimizing sequence of Π. Obviously m < +∞. Thus Π(u k , n k ) is bounded above. By Lemma 4, we have where ε > 0 is small and C i > 0, i = 1, 2, 3 are constants and we have applied the generalized Poincaré inequality ( [4], p281) and the Trace Theorem ( [9], p258) in the last step. Therefore we have F k = I + ∇u k and ∇n k are bounded in L 2 . Now we have that ∇(u k − u 0 ) is bounded in L 2 , by the Poincaré inequality, we have u k is bounded in H 1 . On the other hand, since ∇n k is bounded in L 2 and n k is in H 1 (Ω, S 1 ), we can get that n k is bounded in H 1 . Now since H 1 is a reflexive Banach space and u k and n k are bounded in H 1 , we can find a subsequence of u k and a subsequence of n k such that they are weakly convergent in H 1 . We still denote them as (u k , n k ), and assume u k u, n k n. Now that u k u in H 1 , by Theorem 6, we have det(I + ∇u k ) → det(I + ∇u) in D (Ω). Since det(I + ∇u k ) = 1 a.e., we have det(I + ∇u) = 1 a.e. in Ω 1 . On the other hand, weak convergence in H 1 implies strong convergence in L 2 2 , thus we can find a subsequence of n k that converges point-wise almost everywhere. Therefore we have |n| = 1 almost everywhere. Finally since u k − u 0 ∈ H 1 0|Γ0 , and H 1 0|Γ0 is a closed linear subspace of H 1 , by the Mazur's Theorem, it's weakly closed. Therefore u − u 0 is also in H 1 0|Γ0 , which means u = u 0 on Γ 0 . Similarly we have n = n 0 on Γ 1 . Therefore we have (u, n) ∈ A(u 0 , n 0 ). By Lemma 5, the following function is a convex function of F and P . Therefore by Theorem 1 of section 8.2 of [9], Π is weakly lower semi-continuous. Thus we have Since m is the infimum of Π on A, we conclude that That is, (u, n) is the minimizer of Π on A.
Remark Theorem 6 is crucial in this existence proof. Since in the 3D case, we need the additional condition that adj(F k ) is weakly convergent to get the convergence of det(F k ), the above proof cannot be directly extended to 3D.

Equilibrium equation and its linearized system
Constrained minimization problems can usually be reduced to unconstrained minimization problems by introducing Lagrange multipliers. We consider the following non-dimensionalized energy functional where p is the Lagrange multiplier for the incompressibility constraint det(I +∇u) = 1, which can be interpreted as "pressure", while λ is the Lagrange multiplier for the unity constraint |n| = 1. Assume (u, n, p, λ) minimizes the non-dimensionalized energy (23). Then for any test function v, the function E(ε) = E(u + εv, p, n, λ) has a minimum at ε = 0. Thus we have from which we can get the following equation By similarly taking the variations n → n+εm, p → p+εq, or λ → λ+εµ, we get the following set of Euler-Lagrange equations and H 1 g|Γ0 is the corresponding vector version. We linearize the Euler-Lagrange equation around a solution (u, n, p, λ), and get the following linearized equations: We'll abuse the notation and use · −1 to also denote the norm in the H −1 Γ1 space. The bilinear forms in the linearized system are defined in the following equations

The stress-free state
From the Euler-Lagrange equations, we can get the following strong form partial differential equations: in Ω (37) and the following natural boundary conditions: where the Piola-Kirchhoff stress tensor is Notice that if we choose the displacement to be u ≡ 0 and the director to be n ≡ (0, 1) T , then we have If there is zero body force, equation (34) implies that p must be a constant. Then there is no way to make this σ zero if a = 1. However, it can be verified that the stress-free state can be achieved with and The above observation implies that the reference configuration (strain-free state) for the BTW model is different from the stress-free configuration. It can be verified that if we take the reference state to be the stress-free state instead, then the BTW model becomes the so-called "general BTW" model is the so-called step tensor, and L 0 = L(n 0 ) is the step tensor at the stress-free state.

Well-posedness of the linearized system
By adding (25) and (26) together, and adding (27)-(28) together, we can reduce the linearized system into a standard saddle point system The well-posedness of saddle point systems are well established. Here we quote the so-called Ladyzenskaya-Babuska-Brezzi theorem from [18].
with V and P given Hilbert spaces, L V and L P belonging to V and P respectively. Moreover, a and b are continuous bilinear forms defined on V × V and P × V, respectively. Define the operators Then the operator B is onto if and only if the spaces V and P satisfy the following inf-sup condition: inf Moreover, the mixed problem is well-posed if and only if B is onto and A is invertible.
Remark The operator A is invertible if and only if the following inf-sup condition is satisfied For some mixed system, we can prove the so called "ellipticity" condition which is a stronger condition than the inf-sup condition (48). Thus the ellipticity condition (49) together with (47), will give us a sufficient condition for the well-posedness of the saddle point system (45)-(46).
is satisfied if and only if the corresponding inf-sup conditions for b 1 (q, v) and b 2 (µ, m) are satisfied.
Proof. First if the inf-sup condition for b(q,ṽ) is satisfied, then by Theorem 8, we know that the operator is onto. Therefore, the operator and the operator are both onto. Thus by Theorem 8, the inf-sup conditions for b 1 (q, v) and b 2 (µ, m) are satisfied. Conversely, if the inf-sup conditions b 1 (q, v) and b 2 (µ, m) are satisfied, then the operators B 1 and B 2 are both onto. Thus the operator B is also onto. Therefore we have the inf-sup condition (32) is exactly the one in the incompressible elasticity [18], while b 2 (µ, m) in (33) is exactly the one in the harmonic map problem [12].
For the inf-sup condition for b 1 (q, v), it's well-known that it's satisfied at the strain-free state, where F = I, and the inf-sup condition is reduced to the one in the Stokes problem: Since the stress-free state has constant F matrix, by change of variables, it's easy to verify that the inf-sup condition for b 1 (q, v) is satisfied at the stress-free state, as well. In the general case that u = 0 and F is not a constant, analytical verification of the inf-sup condition for b 1 (q, v) can be very difficult.
As for the inf-sup condition for b 2 (µ, m), a slight modification of the proof in [12] gives us the following theorem: Finally, for the ellipticity condition for the bilinear form a(w,ṽ), it's generally very complicated to verify due to the complexity of the expressions of a 1 (·, ·), a 2 (·, ·) and a 3 (·, ·). However, it can be verified that at the stress-free state, the ellipticity condition for a(w,ṽ) is not satisfied. This doesn't necessarily mean that the linearized system (25)-(28) is not well-posed, though. After all, as remarked before, the ellipticity condition is a sufficient condition, not a necessary condition.

Discretization
In the elastomer problem, we have the following variables to solve: • The displacement vector field u and the pressure p, which is also the Lagrange multiplier for the incompressibility constraint det(F ) − 1 = 0, • The director vector field n and the Lagrange multiplier λ for the unity constraint |n| 2 −1 = 0.
The first pair (u, p) is similar to those in the incompressible elasticity, and we'll use the Taylor-Hood element P 2 − P 1 for (u, p). That is continuous piecewise quadratic finite element for u, and continuous piecewise linear finite element for p. The Taylor-Hoold element has been proved ( [18]) to be stable at least at the strain-free state u = 0. The second pair (n, λ) is similar to those in the harmonic map problem, and we will use piecewise linear finite element for both n and λ, as Winther et al. did in [12]. Notice that it's crucial to impose Dirichlet boundary conditions for n and λ at the same boundary to make sure that |n| = 1 is satisfied at all the mesh nodes. Let V h denote the space of continuous piecewise linear functions and V h,

Existence of the discrete minimization problem
Define the admissible set Notice that n h ∈ N h if and only if the function π h (|n h | 2 − 1) ∈ V h,0|Γ1 is identically 0, which means |n h | = 1 at all the mesh nodes. Let ϕ j , j = 1, · · · , N be a basis of V h , and ψ j , j = 1, · · · , M be a basis of V h,0|Γ1 . And define Then g j is a continuous function on W h,0|Γ0 + u 0h × V h,0|Γ1 + n 0h . Therefore A(u 0h , n 0h ) can be written as the intersection of reciprocal images of 0 of the continuous functions g j . Thus it's closed in W h,0|Γ0 + u 0h × V h,0|Γ1 + n 0h . Define the energy functional where F = I + ∇u. Then the discrete formulation of the minimization problem is Find (u h , n h ) ∈ A(u 0h , n 0h ), minimizing Π on A(u 0h , n 0h ).
Remark Similar arguments work for 3D, as well.
Theorem 12. There exists solution to the discrete minimization problem (57).
Proof. Take any (u h , n h ) ∈ A(u 0h , n 0h ), we have by Lemma 11 where ε > 0 is small and C i > 0, i = 1, 2, 3 are constants and we have applied the generalized Poincaré inequality ( [4], p281) and the Trace Theorem ( [9], p258) in the last step. Thus Π(u h , n h ) → ∞ as u h 1 or n h 1 goes to ∞. Therefore its minimum must be achieved at a bounded subset of A(u 0h , n 0h ). On the other hand, since A is the intersection of reciprocal images of 0 of the continuous functions g j , it's a closed set. Now we are minimizing a continuous function Π(u h , n h ) on a closed bounded finite-dimensional set, by the Weierstrass Theorem, we can find (u h , n h ) ∈ A(u 0h , n 0h ) minimizing Π on A(u 0h , n 0h ).
Remark This theorem is the discrete version of Theorem 7. Since we are dealing with finite dimensional spaces, the proof is much simpler. Plus, in the continuous case we were only able to prove Theorem 7 for 2D, but in the discrete case, it's easy to see that the proof above applies to 3D, too.

Equilibrium equation and its linearization
Similar to the continuous case, we start with the following energy functional with Lagrange multipliers Taking the first variation variation on (59) gives us the Euler-Lagrange equations (equilibrium equations) for the constrained minimization problem (57), and taking another variation on the Euler-Lagrange equations gives us the linearized problems.
Thus the equilibrium equations are, find for any test functions (v, m, q, µ) The corresponding linearized problem is, for a given ( is true for any (v, m, q, µ) Here a 1 , a 2 and b 1 are the same as in the continuous case, while and

Well-posedness of the linearized system
As in the continuous case, the well-posedness of the linearized system can be reduced to the verification of the inf-sup conditions for the bilinear forms b 1 (·, ·), b 2 (·, ·) and a(·, ·).
Since we used the Taylor-Hood element P 2 − P 1 for the (u, p) combination, the inf-sup condition for b 1 (·, ·) is at least satisfied at the strain-free state (Proposition 6.1 of [3]). Since at the stress-free state the deformation gradient F is a constant, by change of variables, it's easy to verify that the inf-sup condition for b 1 (·, ·) is satisfied at the stress-free state, as well.
For the inf-sup condition for b 2 (·, ·), we can use the results from [12]. A slight modification of the proof in [12] gives us the following result: Theorem 13. Assume n ∈ H 1 n0|Γ1 (Ω) W 1,∞ (Ω), and n h ∈ V h,n0|Γ1 satisfies |n h | ≥ C > 0 and n h − π h n 1 ≤ γ/| log(h)| 1/2 . Then we can find a positive constant β 2 , independent of h, such that In general, analytical verification of the inf-sup conditions may not be easy. However, in the discrete case, we can relate the inf-sup values to the singular values of certain matrices.
For the general inf-sup value β h in we have the following result ( [3]).
and q, v are the degrees of freedom of q h and v h respectively.
To verify the inf-sup condition or ellipticity condition for a(·, ·) on Ker(B h ), we are interested in calculating the inf-sup valueβ h in and the ellipticity constantα h in as well.
where u, v are the degrees of freedom of u h and v h respectively. Here B is assumed to be full-rank and the matrix Q is defined by the QR decomposition of (BT −1/2 ) T x TÃ y √ x T x y T y (77) SinceB is full-rank, the matrix R ∈ M m×m in the QR decompositioñ is non-singular. Let where w x ∈ R m and z x ∈ R n−m . Then it's easy to verify that Thus there will be no constraint on z x . Therefore where v ∈ R n is the degree of freedom for v h .
Proof. Let f : H −1 Γ1 → V h,0|Γ1 be the map taking v h tov h , and letφ i = f (ϕ i ). It's easy to see that By definition ofφ i , we have Sinceφ i ∈ V h,0|Γ1 , we can writeφ Substituting it into (82) gives

Existence and Uniqueness of the Lagrange multipliers
Le Tallec proved in [15] that if the inf-sup condition is satisfied at u h , then there exists a unique p h such that (u h , p h ) is a solution of the discrete equilibrium equations of the incompressible elasticity. In this section, we will prove similar results for our elastomer problem. The proof uses the following result of Clarke [6] on constrained minimization problems, where ∂g(x) is the generalized gradient introduced in [5].
Theorem 17. Denote J a finite set of integers. We suppose given: E a Banach space, g 0 , g j (j ∈ J) locally Lipschitz functions from E to R, and C a closed subset of E. We consider the following problem: subject to x ∈ C, g j (x) = 0, ∀j ∈ J. (83) Ifx is a local solution of (83) then there exist real numbers r 0 , s j not all zero, and a point ξ in the dual space E of E such that: where N c (x) is the normal cone at C inx, ∂g j is the generalized gradient of g j (x). Proof. Denote where the functions g j are defined as in (55). It's easy to see that Notice that where Dg 1 0 and Dg 2 0 are in W h,0|Γ0 + u 0h × V h,0|Γ1 + n 0h * . We have and Also, We also note that g j (v h , m h ) is continuously differentiable on W h,0|Γ0 + u 0h × V h,0|Γ1 + n 0h , and that Therefore applying Theorem 17, we have: There exists real numbers r 0 , s j not all zero, such that Using (88) and (89), equation (92) becomes Suppose now r 0 = 0. By linearity, and using equations (90), (91), we can then transform (93) to get Since at least one s j is nonzero, at least one of the equations (94) or (95) is in contradiction with the inf-sup conditions. Thus r 0 cannot be zero. We can then divide (93) by r 0 to get That is where we have denoted Equations (97) and (98) are precisely (60)-(61). Since we also have (u h , n h ) ∈ K h × N h , we conclude that (u h , n h , p h , λ h ) is a solution of (60)-(63). Finally, if there were two different p h , their difference would violate the inf-sup condition for b 1 ; and if there were two different λ h , their difference would violate the inf-sup condition for b 2 . So we have the uniqueness of both p h and λ h .

Implementation using FEniCS
All the computations in this project were done using FEniCS ( [10]), which is an open source software for automated solution of differential equations.
We want to solve the equilibrium equations (60)-(63). Assume (ū,n,p,λ) are the degrees of freedom for (u h , n h , p h , λ h ). Replacing the test functions (v, m, q, µ) by the corresponding basis functions, the equilibrium equations (60)-(63) can be looked as a nonlinear equation for (ū,n,p,λ): F (ū,n,p,λ) = 0 (101) Since it's nonlinear, we can use Newton's method to solve. The iterations are, for integers k ≥ 0, It can be verified that the matrix is exactly the matrix corresponding to the left side of the linearized system (64)-(67) (i.e., by replacing (w, l, o, γ) and (v, m, q, µ) by the corresponding basis functions, and replacing (u h , n h , p h , λ h ) by the solution at step k). Each step of (102) corresponds to a linear variational problem, and can be solved by FEniCS . Thus the whole problem can be solved. There is a minor issue, though, about the procedure above. In both the equilibrium equations (60)-(63) and the linearized equations (64)-(67), there are terms with the interpolation operator π h , which is not supported in the FEniCS form language. Thus we can't directly input those bilinear and linear forms in the form file for FEniCS. We solved this issue by first letting FEniCS assemble the matrix A and the right-side vector b without the π h terms, and then assemble the π h terms ourselves and update the matrix A and the vector b accordingly. It turns out that the π h terms are not very difficult to assemble.
to the stress-free state, and it's assumed that at this state the directors are aligned uniformly in the direction n = (0, 1) T The values of the displacement, pressure etc. at the stress-free state were mentioned in section 2.3. For completeness, we rewrite them here. The displacement at the stress-free state is given by The pressure at the stress-free state is given by And λ at the stress-free state is given by Suppose we want the aspect ratio at the stress-free state to be AR n , we need to set Notice that the whole problem is symmetric about both the vertical middle line X = 0.5AR and the horizontal middle line Y = 0.5. Hence we only need to do the computation at the upperright quarter [0.5AR, AR] × [0.5, 1]. By symmetry, we have the following Dirichlet boundary conditions: n = (0, 1) T on X = 0.5AR and Y = 0.5 (112) Also, we model the "clamped pulling" by specifying the following Dirichlet boundary conditions at X = AR: where t ∈ [0, 1], and 1 + M is the largest elongation factor (ratio of the length of the elastomer after stretching and the length of the elastomer before the stretching). We slowly increases t from 0 to 1, to increase the numerical stability of the Newton's method. To make sure that |n| = 1 at all the mesh nodes, we need to impose Dirichlet boundary conditions for λ at the same boundary as n. Therefore we impose λ = (1 − a)/ √ a on X = 0.5AR, X = AR and Y = 0.5.
We took a = 0.6, b = 0.0015, AR n = 1, M = 0.4, with "time" step t = 0.01. We have used uniform triangular meshes in the simulations. For example, Figure 5 shows the mesh with mesh size h = 2 −5 (notice that the mesh is 16 × 16, but recall that this computational domain is only 1/4 of the true domain).   Table 3 lists the inf-sup values at t = 0, where s(A |Ker(B)) is the inf-sup value of A |Ker(B), while e(A |Ker(B)) is the ellipticity constant of A |Ker(B). We can see that the inf-sup values of b 1 and b 2 don't change much with the mesh. This implies that the inf-sup conditions for b 1 and b 2 probably are always satisfied regardless of the mesh. Next, we can see that the inf-sup values for s(A |Ker(B)) are positive for all the 4 meshes, which means the inf-sup conditions are satisfied thus the discrete saddle point systems are well-posed for all the 4 meshes. However, it seems that the inf-sup values s(A |Ker(B)) tend to 0 as h goes to 0. Finally, we can see that the ellipticity constants e(A |Ker(B)) are negative for the given 4 meshes. This means the ellipticity conditions are not satisfied at the initial stress-free state. This is consistent with our analytical observation. Table 4 lists the inf-sup values at the 0.45 0.90 1.40 log 2 ( n 2h − n h 0 / n h − n h/2 0 ) 1.26 1.67 1.88 log 2 ( n 2h − n h 1 / n h − n h/2 1 ) 0.27 0.69 0.93      Figure 6 shows the graph of nominal stress vs. strain, where we have taken the stress-free state as the reference configuration. We can clearly see a plateau when the strain is in (0.10, 0.22). Comparing with Figure 2, we can say that we've recovered the semi-soft elasticity phenomenon. The plateau regions in Figure 2 are more flat than ours, because those graphs are for long and thin elastomers, while our graph is for the elastomer with the aspect ratio AR n = 1. Figure 7,8 and 9 show the directors configuration at the start, bottom, and end of the plateau region (elongation factor is about 1.10, 1.17 and 1.22 respectively). We can see that the directors just start to rotate at the start of the plateau region (Figure 7), while many of the directors have finish the rotations at the end of the plateau region ( Figure 9). Also, we can see from Figure  7,8 and 9 that during the plateau interval (0.10, 0.22), the elastomer domain is mostly blue (low BTW energy); while after the plateau region (Figure 10), the red part (hight BTW energy) starts to dominate the elastomer domain. This means by rotating the directors, the elastomer maintains relatively low BTW energy during the elongation; while after most of the directors have already finished the rotation, the BTW energy starts to increase. This suggests that the plateau region in the stress-strain graph is probably caused by the rotation of the directors. This agrees with the theory of liquid crystal elastomers [19]. Figure 10 shows the final configuration of the directors (elongation factor s = 1.4). We can see that most of the directors have finished the rotation in the final state.

Conclusion and Discussion
We have designed a numerical algorithm to simulate the 2D BTW+Oseen-Frank model of liquid crystal elastomer. We aimed to recover numerically the experimentally observed phenomena such as the stripe-domain phenomenon, the semi-soft elasticity etc. The existence, well-posedness, uniqueness and convergence etc. of the numerical algorithm have been investigated. We have successfully captured the semi-soft elasticity of liquid crystal elastomer, and found it to be related with the rotation of the directors. However, we didn't see the stripe-domain phenomenon. We believe this is mainly because the b value (the coefficient of the Oseen-Frank energy) 0.0015 is relatively large comparing to true value of b in practice, which we estimate to be probably less than 10 −6 . Since the b value is relatively large, the Oseen-Frank energy dominates the BTW energy, while the latter of which is the main driving force for the occurrence of stripe domains. On the other hand, to observe the stripe domain phenomenon, we should use meshes fine enough to resolve the stripes. For example, in the experiment of Zubarev and Finkelmann ([20]), the width of stripe domain divided by the width of the elastomer is around 15µm/5mm = 0.003, which is much smaller than the mesh size 2 −5 of our finest mesh. However, very small b values require small "time" steps to be stable, and this together with very fine meshes would be computationally very expensive. We have tried 512×512 uniform mesh, and the program ran out of memory. In the future, when the machines have much greater computing power, we can try very small b values and very fine mesh, and see whether the stripe domain occurs. Instead of using uniformly fine mesh, we can also use adaptive mesh to reduce the computational cost.
Another direction worth trying is to use Ericksen or Landau-de Gennes' model instead of the Oseen-Frank energy in the numerical simulation. Oseen-Frank energy only allows point defects, while Ericksen or Landau-de Gennes' model allows line and surface defects, as well ( [16]). The stripe domains might have line or surface defects in the transition area between the stripes, thus using Ericksen or Landau-de Gennes' model might have a better chance of capturing the stripe domain phenomenon.