The content of these notes is based on the lectures by Prof. Graeme W. Milton (University of Utah) given in a course on metamaterials in Spring 2007.
Recall the material laminated in the
x
1
{\displaystyle x_{1}}
direction as shown in Figure 1.
[ 1]
Figure 1. A laminate with direction of lamination
x
1
{\displaystyle x_{1}}
.
We showed that we could write
[
−
E
1
D
t
]
=
L
(
x
1
)
⋅
[
−
D
1
E
t
]
{\displaystyle {\begin{bmatrix}-E_{1}\\\mathbf {D} _{t}\end{bmatrix}}={\boldsymbol {L}}(x_{1})\cdot {\begin{bmatrix}-D_{1}\\\mathbf {E} _{t}\end{bmatrix}}}
or,
[
−
⟨
E
1
⟩
⟨
D
t
⟩
]
=
⟨
L
(
x
1
)
⟩
⋅
[
−
D
1
E
t
]
{\displaystyle {\begin{bmatrix}-\langle E_{1}\rangle \\\langle \mathbf {D} _{t}\rangle \end{bmatrix}}=\langle {\boldsymbol {L}}(x_{1})\rangle \cdot {\begin{bmatrix}-D_{1}\\\mathbf {E} _{t}\end{bmatrix}}}
where
L
(
x
1
)
:=
1
ϵ
11
[
−
1
ϵ
1
t
ϵ
t
1
ϵ
11
ϵ
t
t
−
ϵ
t
1
⋅
ϵ
1
t
]
.
{\displaystyle {\boldsymbol {L}}(x_{1}):={\cfrac {1}{\epsilon _{11}}}~{\begin{bmatrix}-1&{\boldsymbol {\epsilon }}_{1t}\\{\boldsymbol {\epsilon }}_{t1}&~~\epsilon _{11}~{\boldsymbol {\epsilon }}_{tt}-{\boldsymbol {\epsilon }}_{t1}\cdot {\boldsymbol {\epsilon }}_{1t}\end{bmatrix}}~.}
We also showed that
[
−
⟨
E
1
⟩
⟨
D
t
⟩
]
=
L
eff
⋅
[
−
D
1
E
t
]
{\displaystyle {\begin{bmatrix}-\langle E_{1}\rangle \\\langle \mathbf {D} _{t}\rangle \end{bmatrix}}={\boldsymbol {L}}_{\text{eff}}\cdot {\begin{bmatrix}-D_{1}\\\mathbf {E} _{t}\end{bmatrix}}}
where
L
eff
:=
1
ϵ
11
eff
[
−
1
ϵ
1
t
eff
ϵ
t
1
eff
ϵ
11
ϵ
t
t
eff
−
ϵ
t
1
eff
⋅
ϵ
1
t
]
.
{\displaystyle {\boldsymbol {L}}_{\text{eff}}:={\cfrac {1}{\epsilon _{11}^{\text{eff}}}}~{\begin{bmatrix}-1&{\boldsymbol {\epsilon }}_{1t}^{\text{eff}}\\{\boldsymbol {\epsilon }}_{t1}^{\text{eff}}&~~\epsilon _{11}~{\boldsymbol {\epsilon }}_{tt}^{\text{eff}}-{\boldsymbol {\epsilon }}_{t1}^{\text{eff}}\cdot {\boldsymbol {\epsilon }}_{1t}\end{bmatrix}}~.}
Therefore, we got the relation
L
eff
=
⟨
L
(
x
1
)
⟩
.
{\displaystyle {\boldsymbol {L}}_{\text{eff}}=\langle {\boldsymbol {L}}(x_{1})\rangle ~.}
This provides the relations
ϵ
11
eff
=
⟨
1
ϵ
11
⟩
−
1
ϵ
1
j
eff
=
⟨
1
ϵ
11
⟩
−
1
⟨
ϵ
1
j
ϵ
11
⟩
ϵ
i
j
eff
=
⟨
ϵ
i
j
−
ϵ
i
1
ϵ
1
j
ϵ
11
⟩
+
⟨
1
ϵ
11
⟩
−
1
⟨
ϵ
i
1
ϵ
11
⟩
⟨
ϵ
1
j
ϵ
11
⟩
,
i
,
j
≠
1
{\displaystyle {\begin{aligned}\epsilon _{11}^{\text{eff}}&=\langle {\cfrac {1}{\epsilon _{11}}}\rangle ^{-1}\\\epsilon _{1j}^{\text{eff}}&=\langle {\cfrac {1}{\epsilon _{11}}}\rangle ^{-1}\langle {\cfrac {\epsilon _{1j}}{\epsilon _{11}}}\rangle \\\epsilon _{ij}^{\text{eff}}&=\langle \epsilon _{ij}-{\cfrac {\epsilon _{i1}~\epsilon _{1j}}{\epsilon _{11}}}\rangle +\langle {\cfrac {1}{\epsilon _{11}}}\rangle ^{-1}\langle {\cfrac {\epsilon _{i1}}{\epsilon _{11}}}\rangle \langle {\cfrac {\epsilon _{1j}}{\epsilon _{11}}}\rangle ~~,i,j\neq 1\end{aligned}}}
where
⟨
f
(
x
1
)
⟩
=
1
V
∫
0
l
f
(
x
1
)
d
x
1
.
{\displaystyle \langle f(x_{1})\rangle ={\cfrac {1}{V}}~\int _{0}^{l}f(x_{1})~{\text{d}}x_{1}~.}
When the off diagonal elements vanish, we get
ϵ
11
eff
=
⟨
1
ϵ
11
⟩
−
1
(Harmonic average)
ϵ
j
j
eff
=
⟨
ϵ
j
j
⟩
(Arithmetic average)
ϵ
a
b
eff
=
0
if
a
≠
b
.
{\displaystyle {\begin{aligned}\epsilon _{11}^{\text{eff}}&=\langle {\cfrac {1}{\epsilon _{11}}}\rangle ^{-1}&&\qquad {\text{(Harmonic average)}}\\\epsilon _{jj}^{\text{eff}}&=\langle \epsilon _{jj\rangle }&&\qquad {\text{(Arithmetic average)}}\\\epsilon _{ab}^{\text{eff}}&=0&&\qquad {\text{if}}~a\neq b~.\end{aligned}}}
The harmonic average corresponds to a situation in which each layer may be
thought of as a capacitor in series while the arithmetic average corresponds
to a situation where the capacitors are in parallel.
Effective Elastic Properties of a Layered Medium
edit
In this section we use the approach of arranging constant fields to find the effective elastic properties of a layered medium in which each layer is anisotropic. The approach used is that of Backus (Backus62 ).
Consider the layered medium shown in Figure 2. In this case the displacement field is continuous across the interfaces between the layers.
Figure 2. An elastic layered medium with direction of lamination
x
1
{\displaystyle x_{1}}
.
The strain field (
ε
{\displaystyle {\boldsymbol {\varepsilon }}}
) is given by
ε
(
x
)
=
1
2
[
∇
u
+
(
∇
u
)
T
]
.
{\displaystyle {\boldsymbol {\varepsilon }}(\mathbf {x} )={\frac {1}{2}}\left[{\boldsymbol {\nabla }}\mathbf {u} +({\boldsymbol {\nabla }}\mathbf {u} )^{T}\right]~.}
In a Cartesian basis (
e
1
,
e
2
,
e
3
{\displaystyle \mathbf {e} _{1},\mathbf {e} _{2},\mathbf {e} _{3}}
) with coordinates
(
x
1
,
x
2
,
x
3
{\displaystyle x_{1},x_{2},x_{3}}
) we can write
ε
i
j
=
1
2
(
∂
u
i
∂
x
j
+
∂
u
j
∂
x
i
)
.
{\displaystyle \varepsilon _{ij}={\frac {1}{2}}\left({\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}}\right)~.}
Therefore the strain components
ε
22
,
ε
33
,
ε
23
{\displaystyle \varepsilon _{22},\varepsilon _{33},\varepsilon _{23}}
are
also continuous across the interfaces. Moreover, these components of
strain are also constant in each layer. Since a piecewise constant field
that is also continuous must be constant, the strain components
ε
22
,
ε
33
,
ε
23
{\displaystyle \varepsilon _{22},\varepsilon _{33},\varepsilon _{23}}
must be constant throughout the laminate.
The tractions (normal components of the stress) at each interface are given by
t
=
σ
⋅
n
⟹
t
i
=
σ
i
j
n
j
{\displaystyle \mathbf {t} ={\boldsymbol {\sigma }}\cdot \mathbf {n} \qquad \implies \qquad t_{i}=\sigma _{ij}~n_{j}}
where
σ
{\displaystyle {\boldsymbol {\sigma }}}
is the stress and
n
{\displaystyle \mathbf {n} }
is the normal at the interface.
Now the tractions must be continuous at the interfaces. Since the normal
components of the stress are piecewise constant in each layer, this implies
that the normal components of the stress must also be constant throughout
the laminate.
We have chosen
n
{\displaystyle \mathbf {n} }
such that
n
=
(
1
,
0
,
0
)
{\displaystyle \mathbf {n} =(1,0,0)}
. Therefore the stress
components
σ
11
,
σ
12
,
σ
13
{\displaystyle \sigma _{11},\sigma _{12},\sigma _{13}}
must be constant.
Recall that the constitutive relation for an anisotropic elastic material
is given by
σ
=
C
:
ε
.
{\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {\mathsf {C}}}:{\boldsymbol {\varepsilon }}~.}
Following the approach that we used for the permittivity, we now write
the constitutive relation in the form
[
σ
n
σ
t
]
=
[
C
n
n
C
n
t
C
t
n
C
t
t
]
[
ε
n
ε
t
]
{\displaystyle {\begin{bmatrix}{\boldsymbol {\sigma }}_{n}\\{\boldsymbol {\sigma }}_{t}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {C}}_{nn}&{\boldsymbol {C}}_{nt}\\{\boldsymbol {C}}_{tn}&{\boldsymbol {C}}_{tt}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {\varepsilon }}_{n}\\{\boldsymbol {\varepsilon }}_{t}\end{bmatrix}}}
where
σ
n
=
[
σ
11
2
σ
12
2
σ
13
]
;
σ
t
=
[
σ
22
σ
33
2
σ
23
]
;
ε
n
=
[
ε
11
2
ε
12
2
ε
13
]
;
ε
t
=
[
ε
22
ε
33
2
ε
23
]
{\displaystyle {\boldsymbol {\sigma }}_{n}={\begin{bmatrix}\sigma _{11}\\{\sqrt {2}}~\sigma _{12}\\{\sqrt {2}}~\sigma _{13}\end{bmatrix}}~;~~{\boldsymbol {\sigma }}_{t}={\begin{bmatrix}\sigma _{22}\\\sigma _{33}\\{\sqrt {2}}~\sigma _{23}\end{bmatrix}}~;~~{\boldsymbol {\varepsilon }}_{n}={\begin{bmatrix}\varepsilon _{11}\\{\sqrt {2}}~\varepsilon _{12}\\{\sqrt {2}}~\varepsilon _{13}\end{bmatrix}}~;~~{\boldsymbol {\varepsilon }}_{t}={\begin{bmatrix}\varepsilon _{22}\\\varepsilon _{33}\\{\sqrt {2}}~\varepsilon _{23}\end{bmatrix}}}
and
C
n
n
=
[
C
1111
2
C
1112
2
C
1113
2
C
1211
2
C
1212
2
C
1213
2
C
1311
2
C
1312
2
C
1313
]
;
C
n
t
=
[
C
1122
C
1133
2
C
1123
2
C
1222
2
C
1233
2
C
1223
2
C
1322
2
C
1333
2
C
1323
]
;
{\displaystyle {\boldsymbol {C}}_{nn}={\begin{bmatrix}C_{1111}&{\sqrt {2}}~C_{1112}&{\sqrt {2}}~C_{1113}\\{\sqrt {2}}~C_{1211}&2~C_{1212}&2~C_{1213}\\{\sqrt {2}}~C_{1311}&2~C_{1312}&2~C_{1313}\end{bmatrix}}~;~~{\boldsymbol {C}}_{nt}={\begin{bmatrix}C_{1122}&C_{1133}&{\sqrt {2}}~C_{1123}\\{\sqrt {2}}~C_{1222}&{\sqrt {2}}~C_{1233}&2~C_{1223}\\{\sqrt {2}}~C_{1322}&{\sqrt {2}}~C_{1333}&2~C_{1323}\end{bmatrix}}~;~~}
C
t
n
=
[
C
2211
2
C
2212
2
C
2213
C
3311
2
C
3312
2
C
3313
2
C
2311
2
C
2312
2
C
2313
]
;
C
t
t
=
[
C
2222
C
2233
2
C
2223
C
3322
C
3333
2
C
3323
2
C
2322
2
C
2333
2
C
2323
]
.
{\displaystyle {\boldsymbol {C}}_{tn}={\begin{bmatrix}C_{2211}&{\sqrt {2}}~C_{2212}&{\sqrt {2}}~C_{2213}\\C_{3311}&{\sqrt {2}}~C_{3312}&{\sqrt {2}}~C_{3313}\\{\sqrt {2}}~C_{2311}&2~C_{2312}&2~C_{2313}\end{bmatrix}}~;~~{\boldsymbol {C}}_{tt}={\begin{bmatrix}C_{2222}&C_{2233}&{\sqrt {2}}~C_{2223}\\C_{3322}&C_{3333}&{\sqrt {2}}~C_{3323}\\{\sqrt {2}}~C_{2322}&{\sqrt {2}}~C_{2333}&2~C_{2323}\end{bmatrix}}~.}
From the major symmetry of
C
{\displaystyle {\boldsymbol {\mathsf {C}}}}
, we see that
C
n
t
=
C
t
n
T
{\displaystyle {\boldsymbol {C}}_{nt}={\boldsymbol {C}}_{tn}^{T}}
. Also,
C
n
n
{\displaystyle {\boldsymbol {C}}_{nn}}
and
C
t
t
{\displaystyle {\boldsymbol {C}}_{tt}}
are symmetric.
Writing the first row out, we get
(1)
σ
n
=
C
n
n
⋅
ε
n
+
C
n
t
⋅
ε
t
⟹
ε
n
=
C
n
n
−
1
⋅
σ
n
−
[
C
n
n
−
1
⋅
C
n
t
]
⋅
ε
t
.
{\displaystyle {\text{(1)}}\qquad {\boldsymbol {\sigma }}_{n}={\boldsymbol {C}}_{nn}\cdot {\boldsymbol {\varepsilon }}_{n}+{\boldsymbol {C}}_{nt}\cdot {\boldsymbol {\varepsilon }}_{t}\qquad \implies \qquad {\boldsymbol {\varepsilon }}_{n}={\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {\sigma }}_{n}-\left[{\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {C}}_{nt}\right]\cdot {\boldsymbol {\varepsilon }}_{t}~.}
From the second row we get
σ
t
=
C
n
t
T
⋅
ε
n
+
C
t
t
⋅
ε
t
.
{\displaystyle {\boldsymbol {\sigma }}_{t}={\boldsymbol {C}}_{nt}^{T}\cdot {\boldsymbol {\varepsilon }}_{n}+{\boldsymbol {C}}_{tt}\cdot {\boldsymbol {\varepsilon }}_{t}~.}
Substituting the expression for
ε
n
{\displaystyle {\boldsymbol {\varepsilon }}_{n}}
from the first row, we get
(2)
σ
t
=
C
n
t
T
⋅
[
C
n
n
−
1
⋅
σ
n
−
C
n
n
−
1
⋅
C
n
t
⋅
ε
t
]
+
C
t
t
⋅
ε
t
=
[
C
n
t
T
⋅
C
n
n
−
1
]
⋅
σ
n
+
[
C
t
t
−
C
n
t
T
⋅
C
n
n
−
1
⋅
C
n
t
]
⋅
ε
t
.
{\displaystyle {\text{(2)}}\qquad {\boldsymbol {\sigma }}_{t}={\boldsymbol {C}}_{nt}^{T}\cdot \left[{\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {\sigma }}_{n}-{\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {C}}_{nt}\cdot {\boldsymbol {\varepsilon }}_{t}\right]+{\boldsymbol {C}}_{tt}\cdot {\boldsymbol {\varepsilon }}_{t}=\left[{\boldsymbol {C}}_{nt}^{T}\cdot {\boldsymbol {C}}_{nn}^{-1}\right]\cdot {\boldsymbol {\sigma }}_{n}+\left[{\boldsymbol {C}}_{tt}-{\boldsymbol {C}}_{nt}^{T}\cdot {\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {C}}_{nt}\right]\cdot {\boldsymbol {\varepsilon }}_{t}~.}
Collecting (1) and (2) we get
[
−
ε
n
σ
t
]
=
[
−
C
n
n
−
1
C
n
n
−
1
⋅
C
n
t
C
n
t
T
⋅
C
n
n
−
1
C
t
t
−
C
n
t
T
⋅
C
n
n
−
1
⋅
C
n
t
]
[
−
σ
n
ε
t
]
.
{\displaystyle {\begin{bmatrix}-{\boldsymbol {\varepsilon }}_{n}\\\\{\boldsymbol {\sigma }}_{t}\end{bmatrix}}={\begin{bmatrix}-{\boldsymbol {C}}_{nn}^{-1}&{\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {C}}_{nt}\\\\{\boldsymbol {C}}_{nt}^{T}\cdot {\boldsymbol {C}}_{nn}^{-1}&~~~{\boldsymbol {C}}_{tt}-{\boldsymbol {C}}_{nt}^{T}\cdot {\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {C}}_{nt}\end{bmatrix}}{\begin{bmatrix}-{\boldsymbol {\sigma }}_{n}\\\\{\boldsymbol {\varepsilon }}_{t}\end{bmatrix}}~.}
Taking a volume average gives
(3)
[
−
⟨
ε
n
⟩
⟨
σ
t
⟩
]
=
[
−
⟨
C
n
n
−
1
⟩
⟨
C
n
n
−
1
⋅
C
n
t
⟩
⟨
C
n
t
T
⋅
C
n
n
−
1
⟩
⟨
C
t
t
−
C
n
t
T
⋅
C
n
n
−
1
⋅
C
n
t
⟩
]
[
−
σ
n
ε
t
]
.
{\displaystyle {\text{(3)}}\qquad {\begin{bmatrix}-\langle {\boldsymbol {\varepsilon }}_{n}\rangle \\\\\langle {\boldsymbol {\sigma }}_{t}\rangle \end{bmatrix}}={\begin{bmatrix}-\langle {\boldsymbol {C}}_{nn}^{-1}\rangle &\langle {\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {C}}_{nt}\rangle \\\\\langle {\boldsymbol {C}}_{nt}^{T}\cdot {\boldsymbol {C}}_{nn}^{-1}\rangle &~~~\langle {\boldsymbol {C}}_{tt}-{\boldsymbol {C}}_{nt}^{T}\cdot {\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {C}}_{nt}\rangle \end{bmatrix}}{\begin{bmatrix}-{\boldsymbol {\sigma }}_{n}\\\\{\boldsymbol {\varepsilon }}_{t}\end{bmatrix}}~.}
If the effective stiffness of the material is given by
⟨
σ
⟩
=
C
eff
:
⟨
ε
⟩
{\displaystyle \langle {\boldsymbol {\sigma }}\rangle ={\boldsymbol {\mathsf {C}}}_{\text{eff}}:\langle {\boldsymbol {\varepsilon }}\rangle }
we can also show that
(4)
[
−
⟨
ε
n
⟩
⟨
σ
t
⟩
]
=
[
−
(
C
n
n
eff
)
−
1
(
C
n
n
eff
)
−
1
⋅
C
n
t
eff
(
C
n
t
eff
)
T
⋅
(
C
n
n
eff
)
−
1
C
t
t
eff
−
(
C
n
t
eff
)
T
⋅
(
C
n
n
eff
)
−
1
⋅
C
n
t
eff
]
[
−
σ
n
ε
t
]
.
{\displaystyle {\text{(4)}}\qquad {\begin{bmatrix}-\langle {\boldsymbol {\varepsilon }}_{n}\rangle \\\\\langle {\boldsymbol {\sigma }}_{t}\rangle \end{bmatrix}}={\begin{bmatrix}-({\boldsymbol {C}}_{nn}^{\text{eff}})^{-1}&({\boldsymbol {C}}_{nn}^{\text{eff}})^{-1}\cdot {\boldsymbol {C}}_{nt}^{\text{eff}}\\\\({\boldsymbol {C}}_{nt}^{\text{eff}})^{T}\cdot ({\boldsymbol {C}}_{nn}^{\text{eff}})^{-1}&~~~{\boldsymbol {C}}_{tt}^{\text{eff}}-({\boldsymbol {C}}_{nt}^{\text{eff}})^{T}\cdot ({\boldsymbol {C}}_{nn}^{\text{eff}})^{-1}\cdot {\boldsymbol {C}}_{nt}^{\text{eff}}\end{bmatrix}}{\begin{bmatrix}-{\boldsymbol {\sigma }}_{n}\\\\{\boldsymbol {\varepsilon }}_{t}\end{bmatrix}}~.}
Comparing (3) and (4) we can show that
C
n
n
eff
=
⟨
C
n
n
−
1
⟩
−
1
C
n
t
eff
=
⟨
C
n
n
−
1
⟩
−
1
⋅
⟨
C
n
n
−
1
⋅
C
n
t
⟩
C
t
t
eff
=
⟨
C
t
t
−
C
t
n
⋅
C
n
n
−
1
⋅
C
n
t
⟩
+
⟨
C
t
n
⋅
C
n
n
−
1
⟩
⋅
⟨
C
n
n
−
1
⟩
−
1
⋅
⟨
C
n
n
−
1
⋅
C
n
t
⟩
.
{\displaystyle {\begin{aligned}{\boldsymbol {C}}_{nn}^{\text{eff}}&=\left\langle {\boldsymbol {C}}_{nn}^{-1}\right\rangle ^{-1}\\{\boldsymbol {C}}_{nt}^{\text{eff}}&=\left\langle {\boldsymbol {C}}_{nn}^{-1}\right\rangle ^{-1}\cdot \left\langle {\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {C}}_{nt}\right\rangle \\{\boldsymbol {C}}_{tt}^{\text{eff}}&=\left\langle {\boldsymbol {C}}_{tt}-{\boldsymbol {C}}_{tn}\cdot {\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {C}}_{nt}\right\rangle +\left\langle {\boldsymbol {C}}_{tn}\cdot {\boldsymbol {C}}_{nn}^{-1}\right\rangle \cdot \left\langle {\boldsymbol {C}}_{nn}^{-1}\right\rangle ^{-1}\cdot \left\langle {\boldsymbol {C}}_{nn}^{-1}\cdot {\boldsymbol {C}}_{nt}\right\rangle ~.\end{aligned}}}
If the material in each layer is isotropic, then the constitutive relation
is
σ
=
C
:
ϵ
=
λ
(
x
1
)
tr
(
ε
)
1
+
2
μ
(
x
1
)
ε
{\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {\mathsf {C}}}:{\boldsymbol {\epsilon }}=\lambda (x_{1})~{\text{tr}}({\boldsymbol {\varepsilon }})~{\boldsymbol {\mathit {1}}}+2~\mu (x_{1})~{\boldsymbol {\varepsilon }}}
where
λ
{\displaystyle \lambda }
is the Lame modulus and
μ
{\displaystyle \mu }
is the shear modulus. In
that case the effective properties of the laminate are
C
1111
eff
=
⟨
1
λ
+
2
μ
⟩
−
1
;
C
1122
eff
=
C
1133
eff
=
⟨
λ
λ
+
2
μ
⟩
⟨
1
λ
+
2
μ
⟩
−
1
C
1212
eff
=
C
1313
eff
=
⟨
1
μ
⟩
−
1
;
C
2323
eff
=
⟨
μ
⟩
;
C
1112
eff
=
C
1113
eff
=
C
1123
eff
=
0
C
2222
eff
=
C
3333
eff
=
⟨
4
μ
λ
(
λ
+
μ
)
(
λ
+
2
μ
)
2
⟩
+
⟨
λ
λ
+
2
μ
⟩
⟨
1
λ
+
2
μ
⟩
−
1
C
2233
eff
=
⟨
2
μ
λ
λ
+
2
μ
⟩
+
⟨
λ
λ
+
2
μ
⟩
2
⟨
1
λ
+
2
μ
⟩
−
1
{\displaystyle {\begin{aligned}C_{1111}^{\text{eff}}&=\langle {\cfrac {1}{\lambda +2~\mu }}\rangle ^{-1}~;~~C_{1122}^{\text{eff}}=C_{1133}^{\text{eff}}=\langle {\cfrac {\lambda }{\lambda +2~\mu }}\rangle ~\langle {\cfrac {1}{\lambda +2~\mu }}\rangle ^{-1}\\C_{1212}^{\text{eff}}&=C_{1313}^{\text{eff}}=\langle {\cfrac {1}{\mu }}\rangle ^{-1}~;~~C_{2323}^{\text{eff}}=\langle \mu \rangle ~;~~C_{1112}^{\text{eff}}=C_{1113}^{\text{eff}}=C_{1123}^{\text{eff}}=0\\C_{2222}^{\text{eff}}&=C_{3333}^{\text{eff}}=\langle {\cfrac {4~\mu ~\lambda ~(\lambda +\mu )}{(\lambda +2~\mu )^{2}}}\rangle +\langle {\cfrac {\lambda }{\lambda +2~\mu }}\rangle ~\langle {\cfrac {1}{\lambda +2~\mu }}\rangle ^{-1}\\C_{2233}^{\text{eff}}&=\langle {\cfrac {2~\mu ~\lambda }{\lambda +2~\mu }}\rangle +\langle {\cfrac {\lambda }{\lambda +2~\mu }}\rangle ^{2}~\langle {\cfrac {1}{\lambda +2~\mu }}\rangle ^{-1}\end{aligned}}}
Laminates with Arbitrary Direction of Lamination n
edit
So far we have dealt with laminates with a single direction of lamination
that was oriented with the
x
1
{\displaystyle x_{1}}
axis. In this section we generalize our
approach to deal with laminates with a normal
n
{\displaystyle \mathbf {n} }
which is not parallel
to the
x
1
{\displaystyle x_{1}}
axis.
Recall that the normal component of
D
{\displaystyle \mathbf {D} }
, i.e.,
D
⋅
n
{\displaystyle \mathbf {D} \cdot \mathbf {n} }
, is
constant and the tangential components of
E
{\displaystyle \mathbf {E} }
are constant throughout
the entire laminate (if there is only one direction of lamination).
Let us introduce the second-order tensor basis
Γ
1
(
n
)
=
n
⊗
n
;
Γ
2
(
n
)
=
1
−
n
⊗
n
=
1
−
Γ
1
(
n
)
.
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )=\mathbf {n} \otimes \mathbf {n} ~;~~{\boldsymbol {\mathit {\Gamma }}}_{2}(\mathbf {n} )={\boldsymbol {\mathit {1}}}-\mathbf {n} \otimes \mathbf {n} ={\boldsymbol {\mathit {1}}}-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )~.}
These are useful because
(5)
Γ
1
(
n
)
⋅
D
(
x
)
=
Γ
1
(
n
)
⋅
⟨
D
⟩
normal component
Γ
2
(
n
)
⋅
E
(
x
)
=
Γ
2
(
n
)
⋅
⟨
E
⟩
tangential components.
{\displaystyle {\text{(5)}}\qquad {\begin{aligned}{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \mathbf {D} (\mathbf {x} )&={\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \langle \mathbf {D} \rangle &&\quad {\text{normal component}}\\{\boldsymbol {\mathit {\Gamma }}}_{2}(\mathbf {n} )\cdot \mathbf {E} (\mathbf {x} )&={\boldsymbol {\mathit {\Gamma }}}_{2}(\mathbf {n} )\cdot \langle \mathbf {E} \rangle &&\quad {\text{tangential components.}}\end{aligned}}}
Therefore,
(6)
Γ
1
(
n
)
⋅
E
(
x
)
=
E
(
x
)
−
⟨
E
⟩
+
Γ
1
(
n
)
⋅
⟨
E
⟩
.
{\displaystyle {\text{(6)}}\qquad {\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \mathbf {E} (\mathbf {x} )=\mathbf {E} (\mathbf {x} )-\langle \mathbf {E} \rangle +{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \langle \mathbf {E} \rangle ~.}
Let us now introduce a polarization field
(7)
P
(
x
)
=
[
ϵ
(
x
)
−
ϵ
0
1
]
⋅
E
(
x
)
=
D
(
x
)
−
ϵ
0
E
(
x
)
{\displaystyle {\text{(7)}}\qquad \mathbf {P} (\mathbf {x} )=[{\boldsymbol {\epsilon }}(\mathbf {x} )-\epsilon _{0}~{\boldsymbol {\mathit {1}}}]\cdot \mathbf {E} (\mathbf {x} )=\mathbf {D} (\mathbf {x} )-\epsilon _{0}~\mathbf {E} (\mathbf {x} )}
where
ϵ
0
{\displaystyle \epsilon _{0}}
is an arbitrary constant.
The volume averaged polarization field is given by
(8)
⟨
P
⟩
=
ϵ
eff
⋅
⟨
E
⟩
−
ϵ
0
⟨
E
⟩
=
(
ϵ
eff
−
ϵ
0
1
)
⋅
⟨
E
⟩
.
{\displaystyle {\text{(8)}}\qquad \langle \mathbf {P} \rangle ={\boldsymbol {\epsilon }}_{\text{eff}}\cdot \langle \mathbf {E} \rangle -\epsilon _{0}~\langle \mathbf {E} \rangle =({\boldsymbol {\epsilon }}_{\text{eff}}-\epsilon _{0}~{\boldsymbol {\mathit {1}}})\cdot \langle \mathbf {E} \rangle ~.}
Define
S
(
x
)
:=
ϵ
0
[
ϵ
0
1
−
ϵ
(
x
)
]
−
1
S
eff
:=
ϵ
0
[
ϵ
0
1
−
ϵ
eff
]
−
1
.
{\displaystyle {\begin{aligned}{\boldsymbol {S}}(\mathbf {x} )&:=\epsilon _{0}~[\epsilon _{0}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}(\mathbf {x} )]^{-1}\\{\boldsymbol {S}}_{\text{eff}}&:=\epsilon _{0}~[\epsilon _{0}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}_{\text{eff}}]^{-1}~.\end{aligned}}}
Then,
(9)
S
(
x
)
⋅
P
(
x
)
=
−
(
ϵ
0
[
ϵ
0
1
−
ϵ
(
x
)
]
−
1
)
⋅
(
[
ϵ
0
1
−
ϵ
(
x
)
]
⋅
E
(
x
)
)
=
−
ϵ
0
E
(
x
)
S
eff
⋅
⟨
P
⟩
=
−
(
ϵ
0
[
ϵ
0
1
−
ϵ
eff
]
−
1
)
⋅
(
[
ϵ
0
1
−
ϵ
eff
]
⋅
⟨
E
⟩
)
=
−
ϵ
0
⟨
E
⟩
.
{\displaystyle {\text{(9)}}\qquad {\begin{aligned}{\boldsymbol {S}}(\mathbf {x} )\cdot \mathbf {P} (\mathbf {x} )&=-\left(\epsilon _{0}~[\epsilon _{0}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}(\mathbf {x} )]^{-1}\right)\cdot \left([\epsilon _{0}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}(\mathbf {x} )]\cdot \mathbf {E} (\mathbf {x} )\right)=-\epsilon _{0}~\mathbf {E} (\mathbf {x} )\\{\boldsymbol {S}}_{\text{eff}}\cdot \langle \mathbf {P} \rangle &=-\left(\epsilon _{0}~[\epsilon _{0}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}_{\text{eff}}]^{-1}\right)\cdot \left([\epsilon _{0}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}_{\text{eff}}]\cdot \langle \mathbf {E} \rangle \right)=-\epsilon _{0}~\langle \mathbf {E} \rangle ~.\end{aligned}}}
Applying the projection
Γ
1
(
n
)
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )}
to (7), we get
Γ
1
(
n
)
⋅
P
(
x
)
=
Γ
1
(
n
)
⋅
D
(
x
)
−
ϵ
0
Γ
1
(
n
)
⋅
E
(
x
)
.
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \mathbf {P} (\mathbf {x} )={\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \mathbf {D} (\mathbf {x} )-\epsilon _{0}~{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \mathbf {E} (\mathbf {x} )~.}
Using equations (5)
1
{\displaystyle _{1}}
and (6), we have
Γ
1
(
n
)
⋅
P
(
x
)
=
Γ
1
(
n
)
⋅
⟨
D
⟩
−
ϵ
0
E
(
x
)
+
ϵ
0
⟨
E
⟩
−
ϵ
0
Γ
1
(
n
)
⋅
⟨
E
⟩
.
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \mathbf {P} (\mathbf {x} )={\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \langle \mathbf {D} \rangle -\epsilon _{0}~\mathbf {E} (\mathbf {x} )+\epsilon _{0}~\langle \mathbf {E} \rangle -\epsilon _{0}~{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \langle \mathbf {E} \rangle ~.}
From the definitions (9) we can then write
Γ
1
(
n
)
⋅
P
(
x
)
=
Γ
1
(
n
)
⋅
⟨
D
⟩
+
S
(
x
)
⋅
P
(
x
)
−
S
eff
⋅
⟨
P
⟩
−
ϵ
0
Γ
1
(
n
)
⋅
⟨
E
⟩
.
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \mathbf {P} (\mathbf {x} )={\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \langle \mathbf {D} \rangle +{\boldsymbol {S}}(\mathbf {x} )\cdot \mathbf {P} (\mathbf {x} )-{\boldsymbol {S}}_{\text{eff}}\cdot \langle \mathbf {P} \rangle -\epsilon _{0}~{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \langle \mathbf {E} \rangle ~.}
Define
(10)
V
(
n
,
x
)
:=
Γ
1
(
n
)
⋅
⟨
D
⟩
−
S
eff
⋅
⟨
P
⟩
−
ϵ
0
Γ
1
(
n
)
⋅
⟨
E
⟩
.
{\displaystyle {\text{(10)}}\qquad \mathbf {V} (\mathbf {n} ,\mathbf {x} ):={\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \langle \mathbf {D} \rangle -{\boldsymbol {S}}_{\text{eff}}\cdot \langle \mathbf {P} \rangle -\epsilon _{0}~{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \langle \mathbf {E} \rangle ~.}
Then we have
Γ
1
(
n
)
⋅
P
(
x
)
=
S
(
x
)
⋅
P
(
x
)
+
V
(
n
,
x
)
{\displaystyle {\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \mathbf {P} (\mathbf {x} )={\boldsymbol {S}}(\mathbf {x} )\cdot \mathbf {P} (\mathbf {x} )+\mathbf {V} (\mathbf {n} ,\mathbf {x} )}
or
(11)
V
(
n
,
x
)
=
−
[
S
(
x
)
−
Γ
1
(
n
)
]
⋅
P
(
x
)
.
{\displaystyle {\text{(11)}}\qquad \mathbf {V} (\mathbf {n} ,\mathbf {x} )=-\left[{\boldsymbol {S}}(\mathbf {x} )-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]\cdot \mathbf {P} (\mathbf {x} )~.}
Also, form equations (10) and (8) we have
V
(
n
,
x
)
=
Γ
1
(
n
)
⋅
[
⟨
D
⟩
−
ϵ
0
⋅
⟨
E
⟩
]
−
S
eff
⋅
⟨
P
⟩
=
Γ
1
(
n
)
⋅
⟨
P
⟩
−
S
eff
⋅
⟨
P
⟩
{\displaystyle \mathbf {V} (\mathbf {n} ,\mathbf {x} )={\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \left[\langle \mathbf {D} \rangle -\epsilon _{0}~\cdot \langle \mathbf {E} \rangle \right]-{\boldsymbol {S}}_{\text{eff}}\cdot \langle \mathbf {P} \rangle ={\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\cdot \langle \mathbf {P} \rangle -{\boldsymbol {S}}_{\text{eff}}\cdot \langle \mathbf {P} \rangle }
or,
(12)
V
(
n
,
x
)
=
−
[
S
eff
−
Γ
1
(
n
)
]
⋅
⟨
P
⟩
.
{\displaystyle {\text{(12)}}\qquad \mathbf {V} (\mathbf {n} ,\mathbf {x} )=-\left[{\boldsymbol {S}}_{\text{eff}}-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]\cdot \langle \mathbf {P} \rangle ~.}
Inverting (11) and (12) we have
(13)
P
(
x
)
=
−
[
S
(
x
)
−
Γ
1
(
n
)
]
−
1
⋅
V
(
n
,
x
)
⟨
P
⟩
=
−
[
S
eff
−
Γ
1
(
n
)
]
−
1
⋅
V
(
n
,
x
)
.
{\displaystyle {\text{(13)}}\qquad {\begin{aligned}\mathbf {P} (\mathbf {x} )&=-\left[{\boldsymbol {S}}(\mathbf {x} )-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]^{-1}\cdot \mathbf {V} (\mathbf {n} ,\mathbf {x} )\\\langle \mathbf {P} \rangle &=-\left[{\boldsymbol {S}}_{\text{eff}}-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]^{-1}\cdot \mathbf {V} (\mathbf {n} ,\mathbf {x} )~.\end{aligned}}}
Also, taking the volume average of (13)
1
{\displaystyle _{1}}
, we have
(14)
⟨
P
⟩
=
−
⟨
[
S
(
x
)
−
Γ
1
(
n
)
]
−
1
⟩
⋅
V
(
n
,
x
)
.
{\displaystyle {\text{(14)}}\qquad \langle \mathbf {P} \rangle =-\langle \left[{\boldsymbol {S}}(\mathbf {x} )-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]^{-1}\rangle \cdot \mathbf {V} (\mathbf {n} ,\mathbf {x} )~.}
Therefore, comparing (13)
2
{\displaystyle _{2}}
and (14) and invoking
the arbitrary nature
ϵ
0
{\displaystyle \epsilon _{0}}
, we have
[
S
eff
−
Γ
1
(
n
)
]
−
1
=
⟨
[
S
(
x
)
−
Γ
1
(
n
)
]
−
1
⟩
.
{\displaystyle {\left[{\boldsymbol {S}}_{\text{eff}}-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]^{-1}=\langle \left[{\boldsymbol {S}}(\mathbf {x} )-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]^{-1}\rangle ~.}}
This relation provides us with a means of computing the effective
permittivity of a layered medium oriented at a random angle (given by
the normal
n
{\displaystyle \mathbf {n} }
) with respect to the coordinate basis.
Case 1: Simple or Rank-1 Laminate
edit
Consider the Rank-1 laminate shown in Figure 3. The layers have permittivities alternating between
ϵ
1
{\displaystyle {\boldsymbol {\epsilon }}_{1}}
and
ϵ
2
1
{\displaystyle \epsilon _{2}~{\boldsymbol {\mathit {1}}}}
. The volume fraction of phase
1
{\displaystyle 1}
is
f
1
{\displaystyle f_{1}}
while that of phase
2
{\displaystyle 2}
is
f
2
{\displaystyle f_{2}}
such that
f
1
+
f
2
=
1
{\displaystyle f_{1}+f_{2}=1}
.
Figure 3. A Rank-1 laminate.
Recall that
S
(
x
)
=
ϵ
0
[
ϵ
0
1
−
ϵ
(
x
)
]
−
1
.
{\displaystyle {\boldsymbol {S}}(\mathbf {x} )=\epsilon _{0}~[\epsilon _{0}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}(\mathbf {x} )]^{-1}~.}
Let us take the limit as
ϵ
0
→
ϵ
2
{\displaystyle \epsilon _{0}\rightarrow \epsilon _{2}}
. Since
ϵ
(
x
)
=
ϵ
2
1
{\displaystyle \epsilon (\mathbf {x} )=\epsilon _{2}~{\boldsymbol {\mathit {1}}}}
in phase
2
{\displaystyle 2}
, we have
S
(
x
)
→
∞
in phase
2
.
{\displaystyle {\boldsymbol {S}}(\mathbf {x} )\rightarrow \infty \qquad {\text{in phase}}~2~.}
Therefore,
[
S
(
x
)
−
Γ
1
(
n
)
]
−
1
→
0
in phase
2
.
{\displaystyle \left[{\boldsymbol {S}}(\mathbf {x} )-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]^{-1}\rightarrow 0\qquad {\text{in phase}}~2~.}
Hence, right hand side of
[
S
eff
−
Γ
1
(
n
)
]
−
1
=
⟨
[
S
(
x
)
−
Γ
1
(
n
)
]
−
1
⟩
{\displaystyle \left[{\boldsymbol {S}}_{\text{eff}}-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]^{-1}=\langle \left[{\boldsymbol {S}}(\mathbf {x} )-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]^{-1}\rangle }
reduces to an average only over phase
1
{\displaystyle 1}
. If we define
S
1
:=
ϵ
0
[
ϵ
0
1
−
ϵ
1
]
−
1
→
ϵ
2
[
ϵ
2
1
−
ϵ
1
]
−
1
{\displaystyle {\boldsymbol {S}}_{1}:=\epsilon _{0}~[\epsilon _{0}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}_{1}]^{-1}\rightarrow \epsilon _{2}~[\epsilon _{2}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}_{1}]^{-1}}
we get
(15)
[
S
eff
−
Γ
1
(
n
)
]
−
1
=
f
1
[
S
1
−
Γ
1
(
n
)
]
−
1
.
{\displaystyle {\text{(15)}}\qquad \left[{\boldsymbol {S}}_{\text{eff}}-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]^{-1}=f_{1}~\left[{\boldsymbol {S}}_{1}-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]^{-1}~.}
Taking the inverse of both sides of (15) gives
f
1
[
S
eff
−
Γ
1
(
n
)
]
=
S
1
−
Γ
1
(
n
)
{\displaystyle f_{1}~\left[{\boldsymbol {S}}_{\text{eff}}-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )\right]={\boldsymbol {S}}_{1}-{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )}
or,
f
1
S
eff
=
S
1
−
(
1
−
f
1
)
Γ
1
(
n
)
=
S
1
−
f
2
Γ
1
(
n
)
.
{\displaystyle f_{1}~{\boldsymbol {S}}_{\text{eff}}={\boldsymbol {S}}_{1}-(1-f_{1})~{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )={\boldsymbol {S}}_{1}-f_{2}~{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )~.}
Since
S
eff
=
ϵ
0
[
ϵ
0
1
−
ϵ
eff
]
−
1
→
ϵ
2
[
ϵ
2
1
−
ϵ
eff
]
−
1
{\displaystyle {\boldsymbol {S}}_{\text{eff}}=\epsilon _{0}~[\epsilon _{0}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}_{\text{eff}}]^{-1}\rightarrow \epsilon _{2}~[\epsilon _{2}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}_{\text{eff}}]^{-1}}
we then have
f
1
ϵ
2
[
ϵ
2
1
−
ϵ
eff
]
−
1
=
ϵ
2
[
ϵ
2
1
−
ϵ
1
]
−
1
−
f
2
Γ
1
(
n
)
.
{\displaystyle {f_{1}~\epsilon _{2}~[\epsilon _{2}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}_{\text{eff}}]^{-1}=\epsilon _{2}~[\epsilon _{2}~{\boldsymbol {\mathit {1}}}-{\boldsymbol {\epsilon }}_{1}]^{-1}-f_{2}~{\boldsymbol {\mathit {\Gamma }}}_{1}(\mathbf {n} )~.}}
This is the formula of Tartar, Murat, Lurie, and Cherkaev and can be shown to be equivalent to the Backus formula.
↑ The discussion in this lecture is based on Milton02 . Please consult that book for more details and references. The method of Backus (Backus62 ) (see also Postma55 and Tartar76 ) has been used.
[Backus62] G. E. Backus. Long-wave elastic anisotropy produced by horizontal layering. J. Geophys. Res. , 67:4427--4440, 1962.
[Milton02] G. W. Milton. Theory of Composites . Cambridge University Press, New York, 2002.
[Postma55] G. W. Postma. Wave propagation in a stratified medium. Geophysics , 20:780--806, 1955.
[Tartar76] L. Tartar. Estimation de coefficients homogeneises. In R. Glowinski and J. L. Lions, editors, Computer Methods in Applied Sciences and Engineering , pages 136--212. Springer-Verlag, Berlin, 1976.