6. The Fundamental Equations
of Continuum Mechanics

Initialization

Needs["TensorCalculus4`Tensorial`"]

Needs["TContinuumMechanics2`TContinuumMechanics`"]

Needs["DrawGraphics`DrawingMaster`"] <br />

numequ = 1 ;

oldflavors = IndexFlavors ;

ClearIndexFlavor/@oldflavors ;

DeclareIndexFlavor[{black, Black}, {red, Red}, {green, ForestGreen}, {star, SuperStar}, {blue, Blue}, {hat, OverHat}, {tilde, OverTilde}, {bar, OverBar}]

DeclareBaseIndices[{1, 2, 3}] ;

SymbolSpaceDimension[] = NDim ;

SetScalarSingleCovariantD[False]

The following values are taken as an example of a (red) basis frame :

SetTensorValues[gdd[i, j], IdentityMatrix[3]]

RedBasisVectors = ({{1, 1, 0}, {0, 1, 2}, {1, 0, 1}}) ; RedMetric = ( {{2, 1, 1}, {1, 5, 2}, {1, 2, 2}} ) ;

SetTensorValueRules[d[red @ i], RedBasisVectors]

SetTensorValueRules[u[red @ i], Inverse[RedBasisVectors//Transpose]]

SetTensorValueRules[gdd[red @ i, red @ j], RedMetric, red]

TensorSymmetry[g, 2] = Symmetric[1, 2] ;

d := dru[red @ i] d[red @ i]

d := dsu[red @ j] d[red @ j]

d := dtu[red @ k] d[red @ k]

and from equation (2.5) :

Tensor[ℰ§, {Void, Void}, {i_, j_}] := (PartialD[ud[red @ i], red @ j] + PartialD[ud[red @ j], red @ i])/2

E (or E§) is a completely symmetric tensor :

TensorSymmetry[ℰ, 2] = Symmetric[1, 2] ;

TensorSymmetry[ℰ§, 2] = Symmetric[1, 2] ;

The ChristoffelSymbol Γ is also symmetric with respect to its secons and third indices :

TensorSymmetry[Γ, 3] = Symmetric[2, 3] ;

 = uu[red @ i] d[red @ i] ;

Kinematic Relations

The kinematic equation have been presented in the section Strain Tensor of Chapter 2, in a non tensorial form, involving partial derivatives. Using the results of Chapter 5 on the covariant derivatives, and the fact that we have tensors, it can be easily shown that we can identify in cartesian coordinates the partial derivatives with the covariant derivatives and then see that the kinematic equations take the form

      ℰ_ (ij)^(ij) == 1/2 (u_i^i_ (; j) + u_j^j_ (; i) + u_k^k_ (; j) u_k^k_ (; i))       (6.1)

and the linear version of the kinematic equation may be written,

(*6.2*)ResultFrame[res62 = ℰdd[red @ i, red @ j] == (ℰ§dd[red @ i_, red @ j_] = 1/2 (CovariantD[ud[red @ i], red @ j] + CovariantD[ud[red @ j], red @ i])) ]

      ℰ_ (ij)^(ij) == 1/2 (u_i^i_ (; j) + u_j^j_ (; i))       (6.2)

This equation is now valid in any coordinate system. The above equation is symmetrical in its indices, but what is the meaning of the antisymmetric expression 1/2 (u_ (i  ;  j)^(i  ;  j)-u_ (j  ;  i)^(j  ;  i)). Consider the tensor

ωu[k_] := -1/2CovariantD[ud[i], j] euuu[i, j, k]

and we have the identity :

(ωu[k]//ToFlavor[red, black]) d[red @ k] == 1/2TCurl[, g, red/@{j, i, k}, e][]

-1/2 u_i^i_ (; j) e_ (ijk)^(ijk) _k^k == 1/2 u_i^i_ (; j) e_ (jik)^(jik) _k^k

that is to say,

(*6.3*)ResultFrame[res63 = Tensor[ω, {red @ i}, {Void}] d[red @ i] == 1/2curl[] ]

      _i^i ω_i^i == 1/2 curl[u_i^i _i^i]       (6.3)

In cartesian coordinates,

Tensor[ω, {3}, {Void}] == (ωu[3]//SumExpansion[{i, j}]//LeviCivitaOrder[e]//LeviCivitaSimplify[e, ε][g, black])/.PermutationSymbolRule[ε]//Factor

ω_3^3 == -(u_i^i_ (; j) ε_ (3ij)^(3ij))/(2 (g)^1/2)

u_i^i_ (; j) ε_ (3ij)^(3ij) + 2 ω_3^3 == 0

In this expression, ω_3^3is the average rotation of a deformable volume element in the plane {1,2}.

Compatibility conditions

There are three displacements u_ i^i, but six different strain components ℰ_ (i  j)^(i  j). In fact if we calculate the second derivatives of ℰ_ (i  j)^(i  j), and there combination

res0 = CovariantD[ℰdd[red @ i, red @ j], {red @ k, red @ l}] euuu[red @ i, red @ k, red @ m] euuu[red @ j, red @ l, red @ n]

res = res0/.{CovariantD[ℰdd[red @ i, red @ j], {l1_, l2_}] :>CovariantD[(1/2 (CovariantD[ud[red @ i], red @ j] + CovariantD[ud[red @ j], red @ i])), {l1, l2}]}

ℰ_ (ij)^(ij) _ (; kl) e_ (ikm)^(ikm) e_ (jln)^(jln)

1/2 (u_i^i_ (; jkl) + u_j^j_ (; ikl)) e_ (ikm)^(ikm) e_ (jln)^(jln)

we find other conditions to which the ℰ_ (i  j)^(i  j) are subjected. In flat space, u_i^i_ (; jkl) is symmetrical with respect to its covariant derivatives, so that we verify that it gives zero,

We can use CovariantDOrdering if the space has a zero Gaussian curvature !

{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}

The compatibility conditions are then :

(*6.4*)ResultFrame[res64 = CovariantD[ℰdd[red @ i, red @ j], {red @ k, red @ l}] euuu[red @ i, red @ k, red @ m] euuu[red @ j, red @ l, red @ n] == 0]

      ℰ_ (ij)^(ij) _ (; kl) e_ (ikm)^(ikm) e_ (jln)^(jln) == 0      (6.4)

The above matrix {m,n} is symmetrical and the six equations corresponding to {m,n}={1,1},{2,2},{3,3},{1,2},{2,3},{3,1} give zero.

CovariantD[ℰdd[red @ i, red @ j], {red @ k, red @ l}] euuu[red @ i, red @ k, red @ m] euuu[red @ j, red @ l, red @ n] ;

restable =  {res[[1, 1]] == 0, res[[2, 2]] == 0, res[[3, 3]] == 0, res[[1, 2]] == 0, res[[2, 3]] == 0, res[[3, 1]] == 0}//TableForm

We can use CovariantDOrdering if the space has a zero Gaussian curvature !

ℰ_ (22)^(22) _ (; 33) - 2 ℰ_ (23)^(23) _ (; 23) + ℰ_ (33)^(33) _ (; 22) == 0
ℰ_ (11)^(11) _ (; 33) - 2 ℰ_ (13)^(13) _ (; 13) + ℰ_ (33)^(33) _ (; 11) == 0
ℰ_ (11)^(11) _ (; 22) - 2 ℰ_ (12)^(12) _ (; 12) + ℰ_ (22)^(22) _ (; 11) == 0
-ℰ_ (12)^(12) _ (; 33) + ℰ_ (13)^(13) _ (; 23) + ℰ_ (23)^(23) _ (; 13) - ℰ_ (33)^(33) _ (; 12) == 0
-ℰ_ (11)^(11) _ (; 23) + ℰ_ (12)^(12) _ (; 13) + ℰ_ (13)^(13) _ (; 12) - ℰ_ (23)^(23) _ (; 11) == 0
ℰ_ (12)^(12) _ (; 23) - ℰ_ (13)^(13) _ (; 22) - ℰ_ (22)^(22) _ (; 13) + ℰ_ (23)^(23) _ (; 12) == 0

In cartesian coordinates the compatibility equations are (the covariant derivatives being then identical to partial derivatives):

∂^2ℰ_y^y/∂z∂z + ∂^2ℰ_z^z/∂y∂y == ∂^2γ_ (yz)^(yz)/∂y∂z
∂^2ℰ_x^x/∂z∂z + ∂^2ℰ_z^z/∂x∂x == ∂^2γ_ (xz)^(xz)/∂x∂z
∂^2ℰ_x^x/∂y∂y + ∂^2ℰ_y^y/∂x∂x == ∂^2γ_ (xy)^(xy)/∂x∂y
2 ∂^2ℰ_z^z/∂x∂y + ∂^2γ_ (xy)^(xy)/∂z∂z == ∂^2γ_ (xz)^(xz)/∂y∂z + ∂^2γ_ (yz)^(yz)/∂x∂z
2 ∂^2ℰ_x^x/∂y∂z + ∂^2γ_ (yz)^(yz)/∂x∂x == ∂^2γ_ (xy)^(xy)/∂x∂z + ∂^2γ_ (xz)^(xz)/∂x∂y
2 ∂^2ℰ_y^y/∂x∂z + ∂^2γ_ (xz)^(xz)/∂y∂y == ∂^2γ_ (xy)^(xy)/∂y∂z + ∂^2γ_ (yz)^(yz)/∂x∂y

In the covariant form we would have written :

ℰ_y^y_ (; zz) + ℰ_z^z_ (; yy) == γ_ (yz)^(yz) _ (; yz)
ℰ_x^x_ (; zz) + ℰ_z^z_ (; xx) == γ_ (xz)^(xz) _ (; xz)
ℰ_x^x_ (; yy) + ℰ_y^y_ (; xx) == γ_ (xy)^(xy) _ (; xy)
2 ℰ_z^z_ (; xy) + γ_ (xy)^(xy) _ (; zz) == γ_ (xz)^(xz) _ (; yz) + γ_ (yz)^(yz) _ (; xz)
2 ℰ_x^x_ (; yz) + γ_ (yz)^(yz) _ (; xx) == γ_ (xy)^(xy) _ (; xz) + γ_ (xz)^(xz) _ (; xy)
2 ℰ_y^y_ (; xz) + γ_ (xz)^(xz) _ (; yy) == γ_ (xy)^(xy) _ (; yz) + γ_ (yz)^(yz) _ (; xy)

Condition of Equilibrium and Equation of Motion

Figure 6.2

[Graphics:HTMLFiles/index_82.gif]

Figure 6.2 Stress acting on a volume element.

The edges of the infinitesimal block of material (figure), are the vectors,

{d, d, d}

{dr_i^i _i^i, ds_j^j _j^j, dt_k^k _k^k}

The face on the right hand side is dA_ l^l equal to,

(d×d//LinearBreakout[Cross][d[_], u[_]]//CrossProductExpansion[, e, h]) . d[red @ l] ;

dd[red @ l_] = %//EvaluateDotProducts[, g]

ds_j^j dt_k^k e_ (jkl)^(jkl)

The difference between the forces applied on opposite faces are then given by,

CovariantD[dAd[red @ l] σuu[red @ l, red @ m] d[red @ m], red @ i] dru[red @ i]/.CovariantD[dAd[red @ l], red @ i] →0//CovariantDSimplify[, g, e]

force1 = Coefficient[%/.dAd[red @ l] →dd[red @ l], dru[red @ i] dsu[red @ j] dtu[red @ k]]

σ_ (lm)^(lm) _ (; i) dA_l^l dr_i^i _m^m

σ_ (lm)^(lm) _ (; i) e_ (jkl)^(jkl) _m^m

CovariantDSimplify[, g, e][dA_l^l dr_i^i (σ_ (lm)^(lm) _ (; i) _m^m + _m^m_ (; i) σ_ (lm)^(lm))]

σ_ (lm)^(lm) _ (; i) dA_l^l dr_i^i _m^m

Similar expression are obtained for the two others opposite faces, and globally we have

force2 = force1/.{i→j, j→k, k→i}

force3 = force2/.{i→j, j→k, k→i}

σ_ (lm)^(lm) _ (; j) e_ (kil)^(kil) _m^m

σ_ (lm)^(lm) _ (; k) e_ (ijl)^(ijl) _m^m

res1 = (force1 + force2 + force3) dru[red @ i] dsu[red @ j] dtu[red @ k]//Factor

dr_i^i ds_j^j dt_k^k (σ_ (lm)^(lm) _ (; k) e_ (ijl)^(ijl) + σ_ (lm)^(lm) _ (; i) e_ (jkl)^(jkl) + σ_ (lm)^(lm) _ (; j) e_ (kil)^(kil)) _m^m

This expression is identical to (see also the section treating the Gauss theorem in three dimensions, for a similar result)

res2 = dru[red @ i] dsu[red @ j] dtu[red @ k] CovariantD[σuu[red @ l, red @ m], red @ l] eddd[red @ i, red @ j, red @ k] d[red @ m]

σ_ (lm)^(lm) _ (; l) dr_i^i ds_j^j dt_k^k e_ (ijk)^(ijk) _m^m

and,

(res1//SumExpansion[red @ i, red @ j, red @ k, red @ l]//LeviCivitaOrder[e]//Simplify) == (res2//SumExpansion[red @ i, red @ j, red @ k, red @ l]//LeviCivitaOrder[e]//Simplify)

True

In addition, we see that dr_ i^i ds_ j^j dt_ k^k e_ (i  j  k)^(i  j  k) is the volume  dV of the block, so that the resultant of all the stresses is

TotalStress = (res2/.dru[red @ i] dsu[red @ j] dtu[red @ k] eddd[red @ i, red @ j, red @ k] →d)/.{m→i, l→j}

d σ_ (ji)^(ji) _ (; j) _i^i

Superposed to the stress we can have an external force

ExternalForce = d Xu[red @ i] d[red @ i]

d X_i^i _i^i

The system is in equilibrium if

res = TotalStress + ExternalForce == 0   //FullSimplify

d (σ_ (ji)^(ji) _ (; j) + X_i^i) _i^i == 0

Condition of Equilibrium :

(*6.5*)ResultFrame[res65 = (TotalForce = res[[1, 2]]) == 0]

      σ_ (ji)^(ji) _ (; j) + X_i^i == 0      (6.5)

The equilibrium of the moments is itself satisfied by the symmetry of the stress tensor (as shown in the above section on the stress tensor).

Dynamic equation

In the limit of small motion the Lagrangian or particle formulation and the Eulerian or field formulation coincide. We can associate a certain time dependent velocity to each particle, but since the particle always stay close to the position it has at rest, the same velocity is also associated with that point in space.

We define the usual notations,

u_i^i_ (, t) _i^i

Tensor[OverDot[u, 2], {red @ i}, {Void}] d[red @ i] == PartialD[Overscript[, .], t]/.PartialD[d[_], t] →0

_i^i Overscript[u, ..] _i^i == u_i^i_ (, tt) _i^i

If  ρ  is the mass density, we have the dynamic equation,

(*6.6*)ResultFrame[res66 = TotalForce[[1]] == (totalforce = ρ Tensor[OverDot[u, 2], {red @ i}, {Void}] - TotalForce[[2]]) ]

      σ_ (ji)^(ji) _ (; j) == -X_i^i + ρ Overscript[u, ..] _i^i      (6.6)

Fundamental Equations of the theory of Elasticity

We have the basic equations describing the equilibrium and the small motions of an elastic solid :

(1) the kinematic relation        ℰ_ (i  j)^(i  j)==1/2 (u_ (i  ;  j)^(i  ;  j)+u_ (j  ;  i)^(j  ;  i))

(2) the equilibrium condition        
X_ i^i+ σ_ (j  i  ;  j)^(j  i  ;  j)==0        or the dynamic equation    σ_ (j  i  ;  j)^(j  i  ;  j)==ρ Overscript[u, ..] _ i^i-X_ i^i
        
(3) the Hooke's laws        
σ_ (i  j)^(i  j) = ℰ_ (l  m)^(l  m) _ (i  j  l  m)^(i  j  l  m)

    and    
σ_ (i  p)^(i  p)==(Ε (ℰ_ (i  p)^(i  p) + (ν g_ (i  p)^(i  p) ℰ_ (m  m)^(m  m))/(1 - 2 ν)))/(ν + 1)    or    σ_ (i  p)^(i  p)==2 μ ℰ_ (i  p)^(i  p)g_ (i  p)^(i  p) ℰ_ (m  m)^(m  m)
    
They form 6+3+6 = 15 component equations, which contain 15 unknowns (6 ℰ_ (i  j)^(i  j), 6 σ_ (i  j)^(i  j), and 3 u_ i^i) .

We recall first here the Hooke's laws (equations (4.3) and (4.11)):

σ§uu[i_, j_] := Tensor[ℰ, {Void, Void}, {l, m}] Tensor[, {i, j, l, m}, {Void, Void, Void, Void}]

If the material is anisotropic, we have from (2) and (3) in the cell grey box,

res = (CovariantD[(σ§uu[i, j]//ToFlavor[red, black])/.ℰdd[red @ i_, red @ j_] →ℰ§dd[red @ i, red @ j], red @ j]//FullSimplify) == totalforce

1/2 ((u_l^l_ (; m) + u_m^m_ (; l)) _ (ijlm)^(ijlm) _ (; j) + (u_l^l_ (; mj) + u_m^m_ (; lj)) _ (ijlm)^(ijlm)) == -X_i^i + ρ Overscript[u, ..] _i^i

If the material is anisotropic but homogeneous, the covariant derivative of the elastic modulus vanishes,

(*6.7*)ResultFrame[res67 = (res/.CovariantD[uuuu[i_, j_, k_, l_], m_] →0)//Simplify]

      X_i^i + 1/2 (u_l^l_ (; mj) + u_m^m_ (; lj)) _ (ijlm)^(ijlm) == ρ Overscript[u, ..] _i^i      (6.7)

If the material is isotropic, we have (we first calculate σ_ (i  j)^(i  j) as a function of the displacements),

(Ε (u_i^i^(; j) + u_j^j^(; i)))/(2 (1 + ν)) - (Ε ν u_p^p_ (; p) g_ (ij)^(ij))/(-1 + ν + 2 ν^2)

res1 = σuu[red @ i, red @ j] == res

σ_ (ij)^(ij) == (Ε (u_i^i^(; j) + u_j^j^(; i)))/(2 (1 + ν)) - (Ε ν u_p^p_ (; p) g_ (ij)^(ij))/(-1 + ν + 2 ν^2)

Introducing the dynamic equation σ_ (j  i  ;  j)^(j  i  ;  j)==ρ Overscript[u, ..] _ i^i-X_ i^i == totalforce :

res = (CovariantD[res1[[2]], red @ j]//CovariantDSimplify[, g, e]//MetricSimplifyD[g]//SymmetricStandardOrder[u, {2, 3}])/.p→j//Simplify ; (*6.8*)

ResultFrame[res68 = res == totalforce]

This fundamental equation can also be expressed using Lame moduli (equation (4.7)),

 LameRule := {Ε→ (μ (3 λ + 2 μ))/(λ + μ), ν→λ/(2 (λ + μ))}

(*6.9*)ResultFrame[res69 = (res68/.LameRule//Simplify) ]

      μ (u_i^i^(; j)) _ (; j) + (λ + μ) (u_j^j^(; i)) _ (; j) + X_i^i == ρ Overscript[u, ..] _i^i      (6.9)

Taking the divergence of the above equation,

res1 = (res69[[1]] - res69[[2]]) d[red @ i]//Expand//TDiv[, g, red @ k]//CovariantDSimplify[, g, e]//KroneckerAbsorb[g]//Simplify

μ (u_i^i^(; j)) _ (; ji) + (λ + μ) (u_j^j^(; i)) _ (; ji) + X_k^k_ (; k) - ρ Overscript[u, ..] _k^k_ (; k)

divkin = res1[[1]] + (res1[[2]]/.{i→j, j→i}) + res1[[Range[3, 4]]]/.{k→i}//CovariantDOrdering//Simplify

We can use CovariantDOrdering if the space has a zero Gaussian curvature !

(λ + 2 μ) (u_i^i^(; j)) _ (; ij) + X_i^i_ (; i) - ρ Overscript[u, ..] _i^i_ (; i)

We can then introduce the Laplace operator notation  Δ  and write using the preceding identification,

so that the preceding equation is simply in absence of external force,

((divkin//LaplaceOp[Δ])/.CovariantD[Xu[red @ i], red @ i] →0) == 0(*6.1*)

ResultFrame[res610 = %[[1, 2]]/(λ + 2 μ) == -%[[1, 1]]/(λ + 2 μ) == -c^2 %[[1, 1]]/ρ]

-ρ Overscript[u, ..] _i^i_ (; i) + (λ + 2 μ) Δ[u] _i^i_ (; i) == 0

where  c has the dimension of a velocity.
A second important equation is obtained when taking the curl of equation (6.9) :

(res69[[1]] - res69[[2]]) d[red @ i]//Expand//TCurl[, g, red/@{k, l, h}, e]//CovariantDSimplify[, g, e]

res1 = %[[1]] (%[[2]]//MetricSimplifyD[g]) %[[3]]//Simplify

(μ (u_l^l^(; j)) _ (; jk) + (λ + μ) u_j^j_ (; ljk) + X_l^l_ (; k) - ρ Overscript[u, ..] _l^l_ (; k)) e_ (klh)^(klh) _h^h

The term u_j^j_ (; ljk) in the parenthesis, symmetrical, together with e_ (klh)^(klh) gives zero, after summation on k and l,

res = (res1/.res1[[1, 2, 2]] →0) == 0

(μ (u_l^l^(; j)) _ (; jk) + X_l^l_ (; k) - ρ Overscript[u, ..] _l^l_ (; k)) e_ (klh)^(klh) _h^h == 0

(*6.11*)ResultFrame[res611 = res/._h^h →1//LaplaceOp[Δ]]

       (μ (u_l^l^(; j)) _ (; jk) + X_l^l_ (; k) - ρ Overscript[u, ..] _l^l_ (; k)) e_ (klh)^(klh) == 0      (6.11)

An up-down swap of the indices gives an equivalent expression

(*6.12*)ResultFrame[res612 = (res611//ExpandAll//UpDownSwapD[red/@{k, l, h}])//Simplify]

       (μ (u_l^l^(; j)) _ (; j)^(; k) + X_l^l^(; k) - ρ Overscript[u, ..] _l^l^(; k)) e_ (klh)^(klh) == 0      (6.12)

When there is no force in equ. (6.11)

res0 = res612/.X_l^l^(; k) →0 ; (*6.13*)

ResultFrame[res613 = (-res0[[1]]/μ/.ρ→μ Overscript[c, ~]^2//Simplify) == 0]

       -((u_l^l^(; j)) _ (; j)^(; k) - Overscript[u, ..] _l^l^(; k) Overscript[c, ~]^2) e_ (klh)^(klh) == 0      (6.13)

with  

Overscript[c, ~] == ρ/μ^(1/2) ;

1) Any solution of equ. (6.9) without external force, must satisfy equ. (6.10) and (6.11).
    Equation (6.13) has a trivial solution u_ (l  ;  k)^(l  ;  k)ε_ (k  l  m)^(k  l  m)= 0, or  curl [u] = 0.
A vector field which satisfies this condition derives from a potential φ, the displacement potential.

Tφ = Tensor[φ] ;

 := TGrad[, red @ i][Tφ]

u[red @ j_] :=  . u[red @ j]//EvaluateDotProducts[, g]//MetricSimplifyD[g]

ContravariantD[uu[red @ j], red @ k] == ContravariantD[u[red @ j], red @ k] == ContravariantD[u[red @ k], red @ j] == ContravariantD[uu[red @ k], red @ j]

u_j^j^(; k) == φ^(; jk) == φ^(; kj) == u_k^k^(; j)

If we introduce this result in (6.9) when X_ i^i= 0, and take into account the relation u_j^j^(; k)==u_k^k^(; j).

res69/.Xu[red @ i] →0 ;

((%[[1, 1]] + (%[[1, 2]]/.u_a_^a_^(; b_) :→u_b^b^(; a)) - %[[2]])/(λ + 2 μ)/.ρ-> (λ + 2 μ) c^2//Simplify) == 0 ; (*6.14*)

ResultFrame[res614 = (%//LaplaceOp[Δ])]

       -c^2 Overscript[u, ..] _i^i + Δ[u] _i^i == 0      (6.14)

2) On the other hand, a solution of equ. (6.10) has a trivial solution u_ (i  ;  i)^(i  ;  i)= 0, or  div [u] = 0. Carried into (6.9) gives

res69/.Xu[red @ i] →0/. u_j^j^(; i) _ (; j) → 0/.ρ-> μ Overscript[c, ~]^2//Simplify ; <br />(*6.15*)

ResultFrame[res615 = (%[[1]]/μ//LaplaceOp[Δ]) - %[[2]]/μ == 0]

       -Overscript[c, ~]^2 Overscript[u, ..] _i^i + Δ[u] _i^i == 0      (6.15)

The general solutions of both equations (6.14) and (6.15) is in cartesian coordinates (u_i^i==u_i^i),

 ud[i] == Ad[i] fd[i][x - c t]

u_i^i == A_i^i f_i^i[-c t + x]

with three arbitrary functions f_i^i. This represents an elastic wave propagating with the velocity c (or Overscript[c, ~] in the case of eq.(6.14)). The solutions with curl [u] = 0 are dilatational waves, while those with   div [u] = 0 are shear waves.

Flow of Viscous Fluids

We define first the operation TotalTimeDerivation[labs, i][T, t] which is the sum of the partial derivative of a tensor T with respect to time t, and of the spatial coordinates in the parameter defined by the labels labs (see the equation for p and a below).

TotalTimeDerivation[labs_, i_][T_, t_] := PartialD[labs][T, t] + (TotalD[T, t]//ExpandTotalD[labs, i])

We adopt here the Eulerian formulation : coordinate system tied to the fluid, so that a particle of the fluid has a fixed position x_ i^i. In the lab frame, this position is y_ i^i and depends generally on time. Its velocity v_ i^i is given by,

vu[red @ i] == TotalD[yu[red @ i], t]

v_i^i == y_i^i/t

and any field quantity p is

TotalTimeDerivation[{y, _, _, _}, red @ i][Tensor[p], t] ;

%//PartialDToDif

p_ (, t) + p_ (, i) y_i^i/t

Note : here we cannot use PartialSum, because time t do not belong to the base indices (it is not relativistic!).

The acceleration is,

Clear[] ;

TotalTimeDerivation[{y, _, _, _}, red @ j][Tensor[], t] ;

 == %/.TotalD[yu[i_], t] ->vu[i]//PartialDToDif

 == _ (, t) + _ (, j) v_j^j

with the components,

(*6.16*)ResultFrame[res616 = au[red @ i] == PartialD[vu[red @ i], t] + vu[red @ j] CovariantD[vu[red @ i], red @ j]]

      a_i^i == v_i^i_ (, t) + v_i^i_ (; j) v_j^j      (6.16)

where we have used the fact that,

() _ (, j) == v_i^i_ (; j) _i^i == v_i^i_ (, j) _i^i + v_i^i _k^k Γ_ (kij)^(kij)

Result (6.15)  contains two terms, v_ (i  , t)^(i  , t) represents the change of velocity when the flow is non stationary, while v_ j^j v_ (i  ;  j)^(i  ;  j) represents the change of velocity when the flow moves during dt to another point where the velocity is different.

Dynamic equation of the fluid flow :

We use result (6.5) derived for solids, but now Overscript[u, ..] _ i^i is the acceleration a_ i^i of the particles of fluid

σ_ (ji)^(ji) _ (; j) == -X_i^i + ρ Overscript[u, ..] _i^i ;   (*  equ . 6.5  *)

%/.Overscript[u, ..] _i^i→a_i^i(*6.17*)

ResultFrame[res617 = %/.a_i^i→v_i^i_ (, t) + v_i^i_ (; j) v_j^j ]

σ_ (ji)^(ji) _ (; j) == ρ a_i^i - X_i^i

      σ_ (ji)^(ji) _ (; j) == ρ (v_i^i_ (, t) + v_i^i_ (; j) v_j^j) - X_i^i      (6.17)

On the other hand, time derivative of the kinematic relation (6.2) ℰ_ (i  j)^(i  j)==1/2 (u_ (i  ;  j)^(i  ;  j)+u_ (j  ;  i)^(j  ;  i)) gives,

TotalD[ℰdd[red @ i, red @ j], t] == 1/2 (CovariantD[vd[red @ i], red @ j] + CovariantD[vd[red @ j], red @ i])

ℰ_ (ij)^(ij)/t == 1/2 (v_i^i_ (; j) + v_j^j_ (; i))

It is interesting here to split this equation into a rate of dilatation e and a rate of distortion .
In chapter 4, we have defined the parameters e and s, and we introduce,

HyddilDEF = {σud[i_, i_] ->3s, ℰud[i_, i_] ->3e} ;

      3 e/t == ℰ_ (mm)^(mm)/t == v_m^m_ (; m)       (6.18)

We define,

Tensor[ e§, {i_, Void}, {Void, j_}] := ℰud[i, j] - e  δud[i, j]

      e_ (ij)^(ij)/t == 1/2 (v_i^i_ (; j) + v_j^j^(; i)) - 1/3 v_m^m_ (; m) δ_ (ij)^(ij)       (6.19)

If the fluid is incompressible, e/t=0 , or v_m^m_ (; m)= 0, so that in this case equ. (6.18) becomes,

(*6.2*)ResultFrame[res620 = distortion[[1]] == (distortion[[2]]/.{CovariantD[vu[red @ m], red @ m] →0}) == TotalD[ℰud[red @ i, red @ j], t]]

      e_ (ij)^(ij)/t == 1/2 (v_i^i_ (; j) + v_j^j^(; i)) == ℰ_ (ij)^(ij)/t      (6.20)

and there is no difference between e_ (i  j)^(i  j) and ℰ_ (i  j)^(i  j).

In a fluid at rest, or also in an inviscid fluid in motion, the pressure is isotropic, while in a viscous fluid in motion, the volume elements are undergoing a deformation and the compressible stresses have different magnitude in different directions. If we still want to define a pressure in the fluid, it can be the average of normal stress, counted positive when compressive. From the rule HyddilDEF which defines the hydrostatic stress s and the cubic dilation e, we have :

Tensor[p] == -s == -1/3σud[red @ m, red @ m]

p == -s == -1/3 σ_ (mm)^(mm)

Introducing the above relation in s_ (i  j)^(i  j) given in chapter 4

sud[i_, j_] = σud[i, j] - s  δud[i, j]/.s→ -Tensor[p] ;

we have,

Tensor[s, {red @ i, Void}, {Void, red @ j}] == sud[red @ i, red @ j]

s_ (ij)^(ij) == p δ_ (ij)^(ij) + σ_ (ij)^(ij)

Finally, when we introduce the constitutive equation (4.16) s_ (i  j)^(i  j)==2 μ Overscript[e, .] _ (i  j)^(i  j),

res = sud[red @ i, red @ j] == 2μ distortion[[1]] == (2μ distortion[[2]]/.CovariantD[vu[i_], i_] →0)

p δ_ (ij)^(ij) + σ_ (ij)^(ij) == 2 μ e_ (ij)^(ij)/t == μ (v_i^i_ (; j) + v_j^j^(; i))

      σ_ (ij)^(ij) →μ (v_i^i^(; j) + v_j^j^(; i)) - p g_ (ij)^(ij)       (6.21)

The Navier-Stokes Equation

The elimination of σ_ (i  j)^(i  j) between equation (6.21) and equation (6.17) leads to the important Navier-Stokes equation for the kinetic evolution of a viscous fluid in presence of external forces :

First we note that we consider an incompressible fluid, e/t=0 , or v_m^m_ (; m)= 0,  so that :

CovariantD[ContravariantD[vu[red @ j], red @ i], red @ j] →0

(v_j^j^(; i)) _ (; j) →0

res617<br />(*6.22*)

ResultFrame[(res622 = (%%[[2]] - %[[2]]//Expand)/.PartialD[vu[red @ i], t] →Tensor[Overscript[v, .], {red @ i}, {Void}]) == 0]

σ_ (ij)^(ij) _ (; j) == μ (v_i^i^(; j)) _ (; j) - p^(; i)

σ_ (ji)^(ji) _ (; j) == ρ (v_i^i_ (, t) + v_i^i_ (; j) v_j^j) - X_i^i

      μ (v_i^i^(; j)) _ (; j) - p^(; i) - ρ v_i^i_ (; j) v_j^j + X_i^i - ρ Overscript[v, .] _i^i == 0      (6.22)

This equation, which contains four unknown, v_ i^i and p, together with the continuity equation v_m^m_ (; m)= 0, determines the flow of the fluid.

Streamlines

The streamlines of the velocity field v_ i^i are the lines everywhere tangent to the vector v. The streamlines passing through all the points of a closed curve C, form a tube called a stream tube.
Let us choose a coordinate system y_ i^i such that,

vu[red @ 2] := 0

vu[red @ 3] := 0

  := vu[red @ 1] d[red @ 1]

d := dru[red @ 2] d[red @ 2]

d := dsu[red @ 3] d[red @ 3]

where the considered curve C is the quadrilateral defined by {dr, ds}. Its area is,

d = d×d//LinearBreakout[Cross][d[_], u[_]]//CrossProductExpansion[, e, 1]

dr_2^2 ds_3^3 e_ (231)^(231) _1^1

and the flux of fluid through C is,

res = d . //EvaluateDotProducts[, g]

dr_2^2 ds_3^3 e_ (231)^(231) v_1^1

If we consider another cross section y_ 1^1 = const, of the stream tube, this cross section is again dA. Along the steamlines, dr_ 2^2 and ds_ 3^3 do not vary, and

(CovariantD[res, red @ 1]/.{CovariantD[dru[red @ 2], red @ 1] →0, CovariantD[dsu[red @ 3], red @ 1] →0}//LeviCivitaOrder[e]) == 0

(dr_2^2 ds_3^3 e_ (123)^(123) v_1^1) _ (; 1) == 0

Due to the continuity equation v_ (m  ;  m)^(m  ;  m) = 0, this is zero. The same amount of fluid passes through all the cross sections of a stream tube. It is interesting to notice that this result valid in the obvious case of an incompressible fluid, needs only the condition v_ (m  ;  m)^(m  ;  m) = div v = 0, to be valid. We will see later a case of application which is not obvious.
    If the fluid flow is stationary, Overscript[v, .] _ i^i= 0, the streamlines are the paths followed by fluid particles during their motion in the flow field. If the flow is not stationary, the shape of the streamlines changes with time. The pathlines x_ j^j = const, are then different fromthe streamlines.

When the fluid is at rest, all the velocity terms cancel in the Navier-Stokes equation. This determines the hydrostatic pressure,

(res622/.{Tensor[v, __] →0, Tensor[Overscript[v, .], __] →0}) == 0

 or     CovariantD[Tensor[p], red @ i] == Xd[red @ i]

-p^(; i) + X_i^i == 0

or  p_ (; i) == X_i^i

Since p is a scalar, X_ i^imust be a gradient field, and its curl is zero, otherwise, the fluid cannot be in equilibrium.

 Xd[red @ i] == CovariantD[Tensor[Ω], red @ i] curl X == 0  or    TCurl[, g, red/@{i, j, k}, e][Xd[red @ i] u[red @ i]][[Range[1, 2]]] == 0

(curl X == 0  or   X_j^j_ (; i) e_ (ijk)^(ijk) == 0) (X_i^i == Ω_ (; i))

When the viscosity is negligeable (like in water or air), the Navier-Stokes equation reads,

res622/.μ→0 ;

-%[[2]] - %[[4]] == %[[3]] + %[[1]]

ρ v_i^i_ (; j) v_j^j + ρ Overscript[v, .] _i^i == -p^(; i) + X_i^i

or from (6.16)   a_ i^i==v_ j^j v_ (i  ;  j)^(i  ;  j)+v_ (i  , t)^(i  , t) == ( v_ i^i)/(t), and lowering the indices,

(*6.23*)ResultFrame[res623 = ρ TotalD[vd[red @ i], t] == Xd[red @ i] - CovariantD[Tensor[p], red @ i]]

      ρ v_i^i/t == -p_ (; i) + X_i^i      (6.23)

Circulation integral of the velocity on a closed curve, vorticity

"∮ .d" == "∮"vd[red @ i]    Tensor[dy, {red[i]}, {Void}]

       res = --------- %[[2]] == TotalD[%[[2]], t]       t

∮ .d == ∮ dy_i^i v_i^i

 --------- ∮ dy_i^i v_i^i == ∮ (v_i^i dy_i^i/t + dy_i^i v_i^i/t) t

∮ v_i^i dy_i^i/t == ∮ d[y_j^j/t] v_i^i == ∮ d[v_i^i] v_i^i == 1/2 {v_i^i v_i^i} _C == 0

Only the second term in the integration remains, so that we have using equ. (6.23)

Finally, we see that for an inviscid fluid moving in a conservative force field,

(*6.24*)ResultFrame[res624 = "∮"vd[red @ i] Tensor[dy, {red @ i}, {Void}] == const]

      ∮ dy_i^i v_i^i == const      (6.24)

for any closed curve moving with the particles of fluid. In particular, this is true for an infinitesimal curve bounding an area dA. From Stokes theorem, (equ.5.8), we deduce that,

=. ; d=. ;

curl[] . d == const

or introducing the vorticity              (*6.25*)

ResultFrame[res625 = ω == curl[]] (*6.26*)

ResultFrame[res626 = ω . d == const ]

curl[] . d == const

or introducing the vorticity   

      ω == curl[]       (6.25)

      ω . d == const      (6.26)

Similarly to the velocity field and streamlines of the flow,  we can define the vorticity field and the vorticity lines or vortex filaments.
The product ω.dA is the flux of vorticity through the cross section dA of a vortex tube passing through all the points of the perimeter C of dA. This flux is constant.
    Unsteady flows frequently start at rest; steady flows frequently approach from infinity with a uniform velocity toward an obstacle. In both case curl(v) = 0, and from equ.(6.26) valid for any area, curl(v) remains throughout zero : the flow field is a gradient field

(*6.27*)ResultFrame[res627 = vd[red @ i] == CovariantD[Tensor[Φ], red @ i]]

      v_i^i == Φ_ (; i)       (6.27)

This kind of velocity field is called a potential flow. Using the continuity equation v_ (m  ;  m)^(m  ;  m) = 0 leads to

TLaplacian[, g, red @ i][Tensor[Φ]] == 0

(Φ^(; i)) _ (; i) == 0

Introducing (6.27) into (6.22) gives

res = (res622/.μ→0/.vu[i_] :→ContravariantD[Tensor[Φ], i]/.Tensor[Overscript[v, .], {red @ i}, {Void}] :→ContravariantD[Tensor[Overscript[Φ, .]], red @ i]) == 0

-p^(; i) - ρ (Φ^(; i)) _ (; j) Φ^(; j) - ρ (Overscript[Φ, .])^(; i) + X_i^i == 0

or equivalently :

res1 = res//LowerIndexD[i, i]

-p_ (; i) - ρ (Overscript[Φ, .]) _ (; i) - ρ Φ_ (; ij) Φ^(; j) + X_i^i == 0

so that Φ_ (; ji) being symmetrical  Φ_ (; ji)= Φ_ (; ij),

res0 = Φ^(; j) Φ_ (; j)//HoldOp[CovariantD]

(res2 = CovariantD[res0, red @ i]) == (res3 = CovariantD[res0, red @ i]//ReleaseHold)

res2 == (res = res3[[1]] + (res3[[2]]//UpDownSwapD[red @ j])//SymmetricStandardOrder[Φ])

Φ^(; j) Φ_ (; j)

(Φ^(; j) Φ_ (; j)) _ (; i) == (Φ^(; j)) _ (; i) Φ_ (; j) + Φ_ (; ji) Φ^(; j)

(Φ^(; j) Φ_ (; j)) _ (; i) == 2 (Φ^(; j)) _ (; i) Φ_ (; j)

result = p_ (; i) == (Solve[res1/.Φ^(; j) Φ_ (; ij) →1/2 (Φ^(; j) Φ_ (; j)) _ (; i), p_ (; i)][[1, 1, 2]]//Expand)

p_ (; i) == -1/2 ρ (Φ^(; j) Φ_ (; j)) _ (; i) - ρ (Overscript[Φ, .]) _ (; i) + X_i^i

which gives an equation for the pressure p

(*6.28*)ResultFrame[res628 = result[[1]] == result[[2, 3]] - ρ CovariantD[Tensor[1/2res0 + Tensor[Overscript[Φ, .]]], red @ i]]

      p_ (; i) == -ρ (1/2 Φ^(; j) Φ_ (; j) + Overscript[Φ, .]) _ (; i) + X_i^i      (6.28)

The boundary condition for (6.28) concerns the velocity component normal to the boundary S of the domain of the fluid flow. If n_ i^i is the unit normal to S, the normal component of v is

vd[red @ i] nu[red @ i] == CovariantD[Tensor[Ω], red @ i] nu[red @ i]

n_i^i v_i^i == Ω_ (; i) n_i^i

This is subject to an integrability condition by the divergence theorem, equ.(5.7)  so that n_ i^i v_ i^i must be chosen such that,

(*6.29*)ResultFrame[res629 = ∫   vu[red @ i] dAd[red @ i] == 0]                                   \S

       ∫   dA_i^i v_i^i == 0      (6.29)                                            S

Seepage Flow

This section is concerned by the seepage flow through a porous medium. It is a randomly multiconnected medium whose channels are randomly obstructed. The quantity that measures how "holed" the medium is due to the presence of these channels is called the porosity of the medium. Another important quantity, analogous to the conductivity of a network of resistors, is the Darcy permeability: it decreases continuously as the porosity decreases, until a critical porosity reached when it vanishes. Close to this critical porosity, the regime is fractal (see Physics and Fractal Structures, Jean-François Gouyet, Springer-Verlag, Berlin, New York, and Masson, 1996), and the seepage flow is described by the model of invasion percolation. All the present studies are supposed to take place far from the critical region, so that the mean-field approaches remain valid.
The flow in the channels are mostly dominated by viscous friction, and its velocity is proportional to the pressure gradient.
Let X_ i^i be the force acting on a unit of the fluid volume, and v the mean velocity (filter velocity) of the fluid flow with respect to an infinitesimal bulk cross section. To make the fluid move, X_ i^i-p_ (;  i)^(;  i) must be different from zero. If the pore geometry has no preferential direction, v follow the direction of the force field X_ i^i-p_ (;  i)^(;  i).

Darcy's law

(*6.3*)ResultFrame[res630 = vd[red @ i] == k (Xd[red @ i] - CovariantD[Tensor[p], red @ i]) == k CovariantD[Tensor[Tensor[Ω] - Tensor[p]], red @ i] ]

      v_i^i == k (-p_ (; i) + X_i^i) == k (-p + Ω) _ (; i)       (6.30)

Some porous media are anisotropic (linear channels, layers,...) and the Darcy's law takes the more general expression

      v_i^i == k_ (ij)^(ij) (-p_ (; j) + X_j^j) == (-p + Ω) _ (; j) k_ (ij)^(ij)       (6.31)

and the velocity field satisfies the continuity equation,

 := vu[red @ i] d[red @ i] (*6.32*)

ResultFrame[res632 = TDiv[, g, red @ i][] == 0  ]

      v_i^i_ (; i) == 0      (6.32)

Combining this condition with equ.(6.31) gives,

res = TDiv[, g, red @ i][kuu[red @ h, red @ j] (Xd[red @ j] - CovariantD[Tensor[p], red @ j]) d[red @ h]] == 0//FullSimplify

(p_ (; ji) - X_j^j_ (; i)) k_ (ij)^(ij) + k_ (ij)^(ij) _ (; i) (p_ (; j) - X_j^j) == 0

But in an homogeneous medium, the Darcy coefficient do not vary in space and k_ (i  j  ;  i)^(i  j  ;  i)=0.

(*6.33*)ResultFrame[res633 = Expand/@(res/.res[[1, 2]] →0)//SymmetricStandardOrder[p, {2}]]

      p_ (; ij) k_ (ij)^(ij) - X_j^j_ (; i) k_ (ij)^(ij) == 0      (6.33)

and finally in the isotropic case,

rule = kuu[red @ i_, red @ j_] →k guu[red @ i, red @ j]

k_ (i_j_)^(i_j_) →k g_ (ij)^(ij)

res = res633/.rule//MetricSimplifyD[g]

k (p^(; j)) _ (; j) - k X_i^i_ (; i) == 0

This gives the Poisson equation,

(*6.34*)ResultFrame[res634 = res634 = res[[1, 2]] == -res[[1, 1]]/.j→i]

       -k X_i^i_ (; i) == -k (p^(; i)) _ (; i)       (6.34)

or in absence of local force, the Laplace equation,

(*6.35*)ResultFrame[res635 = (-res634[[2]]/k//LaplaceOp[Δ]) == 0]

      Δ[p] == 0      (6.35)

Flügge discusses here the fact that k_ (i  j)^(i  j) is  always found symmetric on practical examples, but emphasizes that there is a lack of convincing general proof.
We will admit the symmetry of k.

(*6.36*)ResultFrame[res636 = kuu[red @ i, red @ j] == kuu[red @ j, red @ i]]

      k_ (ij)^(ij) == k_ (ji)^(ji)       (6.36)

Action of the fluid moving in a porous medium, on the solid matrix.

[Graphics:HTMLFiles/index_396.gif]

Figure 6.4  Volume element of a porous medium

The volume in the above figure is supposed infinitesimal, but its dimension must be much larger than the size of the pores. This allows, far from the critical porosity, to make an effective field approximation. The surface of pores in the section {ds, dt} is  φ || ds×dt ||. From equation (4.1), the force due to the the stress acting on the cross section {ds, dt} is,

σuu[red @ l, red @ m] dAd[red @ l] d[red @ m]

dA_l^l _m^m σ_ (lm)^(lm)

the force due to the fluid pressure and acting on the cross section φ dA_ i^i _ i^i is,

p φ dAd[red @ i] u[red @ i]

p φ dA_i^i _i^i

so that the total force acting on the face {ds, dt} is,

d == (σ§ji = σuu[red @ j, red @ i] - guu[red @ j, red @ i] Tensor[Tensor[p  ] Tensor[ φ ]]) dAd[red @ j] d[red @ i]//Factor

d == dA_j^j _i^i (-(p φ) g_ (ji)^(ji) + σ_ (ji)^(ji))

Replacing the TotalForce in equation (6.5) by the above force components we have  (we define σ§ji  = σ_ (ji)^(ji)-g_ (ji)^(ji) p_^ φ_^),

(TotalForce/.σuu[red @ j, red @ i] →σ§ji) == 0 ;

res = %//CovariantDSimplify[, g, e]//MetricSimplifyD[g]

σ_ (ji)^(ji) _ (; j) - (p φ)^(; i) + X_i^i == 0

(*6.37*)ResultFrame[res637 = res[[1, 1]] == -res[[1, 2]] - res[[1, 3]]]

      σ_ (ji)^(ji) _ (; j) == (p φ)^(; i) - X_i^i      (6.37)

Note this equation shows that a constant pressure produces a stress if the porosity φ is variable.
To study a stress problem in presence of seepage flow, we have to consider the three equations, (4.11), (6.2), and (6.37).


Created by Mathematica  (November 27, 2007) Valid XHTML 1.1!