Mathematical analysis and homogenization of the torsion problem

In this paper we address the limitations of the classical formulation of the torsion problem and give a self-contained survey of the function spaces and formulations that are more suitable for analysing the torsional behavior of composite materials. We also prove some homogenization results for the torsion problem in both the periodic and stochastic setting. Our theoretical results are illustrated by numerical examples.


Introduction
The first investigation on the torsion problem goes back to the works of Coulomb and Navier (which included some erroneous conclusions).The most classical results, including applications to a number of cases of technical importance, is mainly due to Saint-Venant.The generalization of the Saint-Venant theory to torsion of bars consisting of different materials joined along their side surfaces has also been known for a long time (see [13], [14] and [15]).This classical theory, which still is popular especially in the elasticity community, requires that the surfaces separating the different materials are sufficiently smooth.Weak and variational formulations of the torsion problem defined on Sobolev spaces enable us to relax this requirement, but are also suitable for other reasons, in particular for obtaining bounds on the torsional rigidity and analyzing the limit torsional behavior of composite materials when the characteristic length of the microstructures (e.g. the diameter of the fibres) approaches 0. This analysis is based on the homogenization method, which was initiated by De Giorgi and Spagnolo in the late 60's and further developed during the last 30 years by many well-known mathematicians (for more information concerning the homogenization theory see e.g. the books [4], [6] and [17]).Estimating this limit behavior is of great practical importance since any direct numerical treatment leads to difficulties when the characteristic size of the inhomogeneity is small compared with the global geometry (see Figure 1).In this paper we address the limitations of the classical formulation of the torsion problem and discuss formulations that are more suitable for analyzing the mechanical behavior of composites materials (see section 2).We also study some function spaces (see Section 3) which are necessary in connection with the analysis of the corresponding dual problem.The rest of the paper considers homogenization of the torsion problem in both the periodic and the stochastic setting (see Section 5 and Section 6).In order to illustrate our theoretical results we even present some concrete numerical computations.
For the purpose of making the paper readable to a broader audience, including scientists within the mathematical elasticity community, as well as pure mathematicians, we have made an effort making the presentation as self-contained as possible starting from the standard equations of elastic equilibrium.

Classical and weak formulations
We consider a locally isotropic composite bar with Lamé constants λ = λ(x, y), μ = (x, y) and cross-section Ω ⊂ R 2 .Let Ω be an open, bounded and connected set of class C 2 with boundary ∂Ω consisting of a finite number of connected components L 1 , L 2 ,...,L m+1 (disconnected from each other) where L m+1 surrounds the others.The torsion problem consists of determining the function ϕ(x, y) such that the displacement components (u, v, w) at each point (x, y, z) are given by ( 1) together with the boundary conditions on the side surface where n = (n x , n y ) is the outward normal to the boundary ∂Ω.In addition we have the symmetry property σ xy = σ yx , σ xz = σ zx and σ yz = σ zy .
Inserting (1) into ( 2) and ( 3) we obtain the stress components This problem can only be understood in classical sense in parts of Ω where the shear modulus μ = μ(x, y) is sufficiently smooth.However, even pioneer works in the mathematical theory of elasticity deal with cases where μ may be discontinuous.For example the book Muskhelisvili [15] considers the case of bars consisting of prismatic (cylindrical) connected components made of different materials and joined along their side surfaces (see also [16]).Assuming that these surfaces are sufficiently smooth, it makes sense to define the normal component (σ xz , σ yz ) • n = σ xz n x + σ yz n y of the stress vector (σ xz , σ yz ) on these surfaces.Then, using the condition that the forces acting on any surfaces, in particular the surfaces that separates different materials, are equal in magnitude and opposite in direction, we find that the normal component (σ Multiplying the equation ( 8) with an arbitrary function v, which is smooth on Ω, we obtain by the Green formula that But by ( 9) and (6) we see that So we conclude that This gives us the following weak formulation of the problem: for all v ∈ H 1 (Ω) (the twist-constant may certainly be cancelled).This problem has a unique solution up to an arbitrary additive constant (which is easily seen by the Lax-Milgram Lemma), even under much more general conditions for the shear modulus μ than described above.In fact it is enough to assume that μ : Ω → R is Lebesgue measurable and that there exist constants α, β such that for all (x, y) ∈ Ω.

On some function spaces
Under the general condition (11) the boundary condition (6) and the equilibrium equation (8) do not have any meaning in the classical sense.This makes it necessary to understand the concept of divergence and normal derivative at the boundary in a more abstract way.For any vector In accordance with the terminology of the distribution theory, we say that in which case we put div u = f (without causing any ambiguity since f necessarily must be unique).
There exists a linear continuous operator γ 0 : Moreover, there exists a linear continuous operator l Ω : For more information about traces of functions belonging to H m (Ω), see Lions and Magenes [8].Let H − 1 2 (∂Ω) denote the dual space of H 1 2 (∂Ω) and let E(Ω) denote the set of all u ∈ L 2 (Ω) such that div u ∈ L 2 (Ω).The following trace result (see Temam [18, p. 9]) gives a generalized meaning of the normal derivative at the boundary u • n| ∂Ω and a generalization of the classical Greens formula.

Lemma 1. There exists a continuous linear operator γ
In the two dimensional case it is usual to define curl u as follows: If u is not a differentiable function, these definitions must be understood in sense of distributions.More precisely, if Here, H −1 (Ω) denotes the dual of the space Similarly as for the divergence, we say that curl in which case we put curl u = f (without causing any ambiguity since f necessarily must be unique).
Let V and H be the spaces defined by where D(Ω) is the space of smooth functions with compact support in Ω.It turns out that H equals the closure of V in L 2 (Ω) (see [18]).Moreover, H ⊥ , the orthogonal complement of H in L 2 (Ω), is given by ( 14) (see [18, p. 15]).The space can be characterized as follows (15) curl where v i ∈ D(Ω) has the property that v i = 1 on L i and v i = 0 on L j for all j = i .This result is the same as that found in [18, Appendix 1] with the only difference that γ n u, v i = 0 is used instead of Li u • n ds = 0.In order to see that such a function v i exists, let {O k } N k=1 be a finite collections of open disks covering L i but small enough such that non of them intersects L j for all j = i.Then, due to the theorem of infinitely differentiable partitions of unity (see e.g.[1, Theorem 3.14]), we can construct functions ψ k ∈ D(R 2 ), with support in O k , such that ψ(x) := N k=1 ψ k (x) = 1 for all x ∈ L i .Thus, the function v i = ψ| Ω has the desired properties.
Of independent interest we note that if Ω is simply connected then (15) reduces to curl This follows directly from ( 12) by putting w = v 1 = 1.
Let us now prove the following characterization of the space H.
Theorem 2. The space H defined in (13) can alternatively be defined as follows: Proof.By (15) we observe that for some positive constant C, so ( 16) If u ∈ H then, due to the fact that V is dense in H (with respect to the L 2 -norm), we can find a sequence . By ( 16) ϕ k is a Cauchy-sequence in H 1 (Ω), which (by the completeness of H 1 (Ω)) converges to some function ϕ ∈ H 1 (Ω).Since we obtain that u = curl ϕ.The continuity of the trace operator γ 0 now gives that Thus, γ 0 (ϕ) = lim k γ 0 (ϕ k ), so, since all γ 0 (ϕ k ) are constant on each L i , the limit γ 0 (ϕ) must possess the same property.The proof is complete.

Dual formulation of the torsion problem
When the shear modulus satisfies the general condition (11) we may derive the weak formulation (10) directly by using the generalized Green formula from Lemma 1 by using the assumption that the stress vector for some ϕ ∈ H 1 (Ω).The fact that (10) Assuming that ϕ is twice differentiable in (17), we obtain that Thus, the classical formulation corresponding to (10) takes the form (20) where k 1 ,...,k m+1 are unknown constants (to be determined as part of the problem).Since the solution of (20) only is determined within an arbitrary additive constant, we might specify one of the constants k 1 ,...,k m+1 .For example, when Ω is simply connected (i.e.m + 1 = 1, as in Figure 1) the constant value of each v ∈ W along the boundary ∂Ω may be set equal to 0. Thus, in this case Note that the resultant torsion moment is given by M = τD, where The constant D, which is called the torsional rigidity, is independent of τ since the solution ϕ of (10) is independent of τ .It is possible to show that where E * is the "energy" (24) In actual calculations of D , using commercially available FE-software, we replace the above space W by a finite dimensional subspace W h ⊂ W and compute the value where ψ h is the FE-solution, i.e. the solution satisfying the dual problem (18) with W replaced by W h .In order to illustrate, let us consider the simple case when Ω = −2, 2 2 and let us compute E * h by using the software Ansys.We let W h be a function space consisting of quadratic polynomials (8 node elements called plane77 in Ansys with element edge length 0.07).Using this space we obtain the value E * h = 17.9939, i.e.D ≈ 35.988 .Note that this is close to the exact value (up to three decimals) found by using Fourier-series, which is 35.994(see e.g.[19, p. 313]).
Proof.Using (11) where and Q is an open bounded domain such that Q ⊃ Ω κ for all κ.This shows that {grad ϕ κ + (−y, x)} and hence {grad ϕ κ } is bounded in L 2 (Ω ∞ ).By subtracting the average value we can always assume that ϕ κ belongs to the space is an equivalent norm on V, by the Poincare inequality, it follows that {ϕ κ } is bounded in V. Thus, by the reflexivity of V there exists a subsequence, still denoted {ϕ κ } , such that Let us now define it follows that {η κ } is bounded in L 2 (Ω ∞ ).Thus, there exists a subsequence, still denoted by {η κ } , and some function

Stochastic homogenization
Let (M, F , λ) be a probability space and T z (z = (x, y) ∈ R 2 ) an ergodic dynamic system on M. A vector field f ∈ L 2 (M) 2 is said to be potential (respectively solenoidal) if its generic realization f (T z ω) is a potential (respectively solenoidal) vector field defined on R 2 .We denote by L 2 pot (M) (respectively L 2 sol (M)) the subspace of L 2 (M) 2 formed by potential (respectively solenoidal) vector fields.We can now define the following space of vector fields with vanishing mean value Let a = a(ω) be a measurable function such that 0 < α ≤ a(ω) ≤ β < ∞ for every ω ∈ M. Now, consider the case when μ = μ κ,ω is of the form μ κ,ω (x, y) = a(T (κx,κy) ω).
Similarly as for (25) we may consider the problem: Find ϕ κ,ω ∈ W Let b be the constant 2 × 2 matrix defined for every η ∈ R 2 by where v η is the solution of the auxiliary problem: Find We are now in the position to state the following stochastic version of (3).

Theorem 4. For almost every
Roughly speaking the proof of this theorem relies on mimicing the proof above in the periodic setting for fixed realizations (with suitable adjustments) and replacing convergences of the type (35) with those obtained from the use of Birkhoff ergodicity theorem (in the form which is usual in connection with stochastic homogenization, see e.g.[6]).
The homogenized matrix b given by (37) is usually impossible to calculate directly.Fortunately, using results given by Bourgeat and Piatnitski [2] we are able to find b approximately by using so-called periodic approximation.Indeed, let the function μ ω given by μ ω (x, y) = a(T (x,y) ω) be restricted to the cube Y δ = [0, δ] 2 and then extended by periodicity to R 2 .We denote this extension μ δ ω .The homogenized matrix b δ ω corresponding to this periodic structure is then given by Note also that since b δ ω is symmetric it is also possible to find this matrix by using that Concerning this fact, see e.g.[6].According to [2] (42 for almost every ω ∈ M.This remarkable fact makes it possible to relate stochastic homogenization to periodic homogenization.For each fixed δ > 0 let b δ denote the expectation of b δ ω , i.e. ( By using that α ≤ μ δ ω ≤ β we obtain from (41) that Hence, by Lebesque dominated convergence theorem, (42) and ( 43) we find that  As an example, let 0 < ρ < 1/2 and consider the probability space (M, F , λ) defined as follows.Each ω ∈ M is identified as a union of discs in R 2 with radius ρ distributed in such a way that there exists a 1-periodic lattice (with period equal to 1 in both directions) where each square contains exactly one disc.It is easy to see that M can be identified with the space Γ = [0, 1] 2 .Moreover, we define the dynamic system T z : M → M by letting T z ω be the translated set T z ω = ω + z .By results from ergodic theory it is possible to show that this dynamic system is ergodic (see e.g. the book [5], page 180).Now, let a = a(ω) be defined by Therefore it is reasonable to use this matrix as an approximation also for the homogenized matrix, i.e. putting b ≈ (1.8)I .This matrix can now be used to find an approximate value D ∞ for torsional rigidity D κ,ω by solving the global problem (39) and inserting this solution into (38).This is certainly a much simpler problem to solve than calculating D κ,ω directly, especially when the number of fibres is large.As an example we may consider the case when Ω κ = Ω ∞ = [0, 1] 2 .Using that the torsional rigidity for a unit square with unit shear modulus is 2.2496 (see [19, p. 313]), we obtain that D ∞ = 2.2496 • (1.8) = 4.05.Thus, by the above homogenization result we may assume that D κ,ω ≈ 4.05 when κ (or equivalently the number of fibers κ 2 ) is sufficiently large.In Figure 2, we have illustrated one such structure for κ = 17 which is randomly generated by using a simple computer program according to the above description of M. For similar examples of numerical computations of random structures we refer to [3].

Figure 1 .
Figure 1.Torsion of composite bar with simply connected cross section

2
and λ with the product measure on Γ, i.e. the measure λ = m × ∞ n=1 m, where m denotes the Lebesgue measure on [0, 1] The realization μ κ,ω (x, y) = a(T (κx,κy) ω), can then be written as(45) μ κ,ω = αχ ω + βχ R 2 \ω ,where χ A denotes the characteristic function (indicator function) of the set A. It is clear from (45) that each ω ∈ M corresponds to a fixed twocomponent composite material with shear modulus equal to α on the discs (fibers) and β on the remaining part (the matrix).The matrix b δ can be found approximately by first computing b δ ω (numerically using the finite element method) for a large number N of different realizations and then taking the average of these values.By letting ρ = 0.3, α = 1, β = 1000 and using this way of calculating b δ for N = 400 realizations, it appears that b 2 ≈ b 4 ≈ b 6 ≈ b 8 ≈ (1.8)I , where I is the identity matrix.The variations between these four matrixes are found to be in third decimal.
On a point of the common boundary between S i and S k with normal n i (pointing out of S i ), the normal component (σ xz , σ yz ) • n i on the set S i , denoted (σ xz , σ yz ) • n i i , must equal on the set S k , which is the negative value of the normal component (σ xz , σ yz ) • n k on the set S k where n k = −n i (the normal is pointing out of S k ), i.e.
(14)e grad ϕ ∈ H ⊥ by(14)and curl v ∈ H by Theorem 2), this integral vanishes.Hence, putting ψ = ψ/τ (for convenience) we find that the corresponding weak minimum formulation takes the form: Find ψ ∈ W ∈ W. Let Ω i denote the region inside the contour L i , and note that the outward normal to the boundary of Ω i is −n except for i = m + 1, where it is n.Then we obtain that has a solution shows that such a stress vector exists.By Theorem 2, (σ xz , σ yz ) = curl ψ for some function ψ ∈ W .Moreover, we observe that yz − x = grad ϕ ∈ H ⊥ .Thus, for all v ∈ V (⊂ W ) we have thatΩ 1 τμ grad ψ • grad v dxdy − Ω grad ϕ • curl v dxdy.