Continuum Mechanics

Renan Cabrera
cabrer7@uwindsor.ca
Physics Department
University of Windsor
Ontario Canada

This notebook presents an application in Continuum Mechanics, based in the classic Landau's book Elasticity Theory.  A much more comprehensive application was done by Jean-Francois in his package Continuum Mechanics.

Initialization

In[1]:=

Needs["TensorCalculus4`Tensorial`"]

There is one animation at the beginning of The strain tensor section that requires the DrawGraphics package. If you don't have that package you can skip the animation.

In[2]:=

Needs["DrawGraphics`DrawingMaster`"]

Needs["DrawGraphics`DrawingArrows`"]

In[4]:=

DeclareIndexFlavor[{red, Red}, {blue, Blue}]

In[5]:=

DeclareBaseIndices[{1, 2, 3}] ;

In[6]:=

DefineTensorShortcuts[{{x, Curlu, u, zero}, 1}, {{δ, g, ℊ, γ, u, ℰ, σ, TubeStrain}, 2}, {{ε, Γ}, 3}] ;

In[7]:=

DeclareZeroTensor[zero]

In[8]:=

labs = {x, δ, g, Γ} ;

In[9]:=

SetTensorValues[δud[red @ i, red @ j], IdentityMatrix[3]]

SetTensorValues[εuuu[red @ i, red @ j, red @ k], PermutationPseudotensor[3]]

SetTensorValues[εddd[red @ i, red @ j, red @ k], PermutationPseudotensor[3]]

Forcing the identity in order to improve the efficiency of the computation

In[12]:=

CovariantD[Tensor[g, _, _], _] = 0 ;

Metric g_ (i j)^(i j) in Euclidean Cylindrical Coordinates

Defining an Euclidean metric in cylindrical coordinates

In[13]:=

(metric = DiagonalMatrix[{1, ρ^2, 1}])//MatrixForm

Out[13]//MatrixForm=

( {{1, 0, 0}, {0, ρ^2, 0}, {0, 0, 1}} )

In[14]:=

(tmetric = ToFlavor[red][CoordinatesToTensors[{ρ, θ, z}][metric]])//MatrixForm

MapThread[SetTensorValues[#1, #2] &, {{gdd[red[μ], red[ν]], guu[red[μ], red[ν]]}, {tmetric, Inverse[tmetric]}}] ;

Out[14]//MatrixForm=

( {{1, 0, 0}, {0, (x_1^1)^2, 0}, {0, 0, 1}} )

Calculating the Christoffel symbols from the metric

In[16]:=

{downΓ, upΓ} = CalculateChristoffels[labs, red, Simplify[#, Trig→False] &] ;

In[17]:=

SetTensorValues[Γddd[red @ a, red @ b, red @ c], downΓ] ;

SetTensorValues[Γudd[red @ a, red @ b, red @ c], upΓ] ;

The strain tensor ℰ_ (i  j)^(i  j)=1/2 (u_ i^i _ (; j)+u_ j^j _ (; i))

The strain tensor is related to the change in the distance between two generic points when they are subject to a displacement. The following amimation shows this distance represented by a the blue vector. (Select and evaluate the thin closed cell for the animation.)

[Graphics:HTMLFiles/index_26.gif]

The strain tensor could be calculated in terms of the covariant derivatives of the displacement vector u_ i^i

In[27]:=

ExpandedStrainTensorRule = udd[red[i_], red[j_]] →ToFlavor[red][1/2 (CovariantD[ud[i], j] + CovariantD[ud[j], i])]

Out[27]=

u_ (i_j_)^(i_j_) →1/2 (u_i^i_ (; j) + u_j^j_ (; i))

An explicit expansion of the strain tensor in cylindrical coordinates is

In[28]:=

ExpandCovariantD[labs, red[a]][udd[red[i], red[j]]/.ExpandedStrainTensorRule]//Expand

SumExpansion[red[a]][ArrayExpansion[red[i], red[j]][%]]

SetTensorValueRules[udd[red[i], red[j]], %]

Out[28]=

-1/2 u_a^a Γ_ (aij)^(aij) - 1/2 u_a^a Γ_ (aji)^(aji) + 1/2 ∂u_i^i/∂x_j^j + 1/2 ∂u_j^j/∂x_i^i

Out[29]=

where we stored the values in a rule associated to the stress tensor

In[31]:=

TensorValueRules[u]

Out[31]=

A comparison with the result found in Landau's book reveals that they are different. They are different because Landau works with the so called "physical components" and here we work with either contravariant of covariant components.

In a next section we are going to need an explicit expression for the divergence of the displacement  u_ (b b)^(b b)

In[32]:=

udd[red @ i, red @ j] * guu[red @ j, red @ i]

%/.ExpandedStrainTensorRule//ExpandCovariantD[labs, red @ a]

DisplacementDivergenceRule = udu[red @ b_Symbol, red @ b_Symbol] → (Expand[%]//SumExpansion[red @ a, red @ i, red @ j])

Out[32]=

g_ (ji)^(ji) u_ (ij)^(ij)

Out[33]=

1/2 g_ (ji)^(ji) (-u_a^a Γ_ (aji)^(aji) + ∂u_i^i/∂x_j^j) + 1/2 g_ (ji)^(ji) (-u_a^a Γ_ (aij)^(aij) + ∂u_j^j/∂x_i^i)

Out[34]=

u_ (b_Symbolb_Symbol)^(b_Symbolb_Symbol) →u_1^1/x_1^1 + ∂u_1^1/∂x_1^1 + ∂u_2^2/∂x_2^2/(x_1^1)^2 + ∂u_3^3/∂x_3^3

Compatibility equations  ℰ_ (i j)^(i j) _ (; k l) ε_ (i k m)^(i k m) ε_ (j l n)^(j l n)=0

The definition of the stress tensor does not allow to calculate the displacements given the values of the stress tensor by itself. More equations are needed to constrain the six values of the stress tensor.

The compatibility equations in cylindrical coordinates are such that each term of the following array must be zero. Only three of them are independent. This operation may take more than a few seconds.

In[35]:=

stepa = CovariantD[ℰdd[red @ i, red @ j], {red @ k, red @ l}] * εuuu[red @ i, red @ k, red @ m] * εuuu[red @ j, red @ l, red @ n]

Out[35]=

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

Out[36]=

Thermodynamic potentials

Three thermodynamical potentials can be defined in terms of the  stress tensor σ_ (i  j)^(i  j)

Internal Energy

dU = T dS + σ_ (i  j)^(i  j)  u_ (i  j)^(i  j)

Helmholtz Free Energy

dF = σ_ (i  j)^(i  j)  u_ (i  j)^(i  j) - S dT

Gibbs Energy

dG = -S dT - u_ (i  j)^(i  j)  σ_ (i  j)^(i  j)

Free energy expansion for isotropic materials in terms of the Lamé moduli μ and λ

In[37]:=

Attributes[μ] = {Constant} ; Attributes[λ] = {Constant} ; Attributes[K] = {Constant} ;

In[38]:=

Out[38]=

F→μ g_ (ik)^(ik) g_ (jm)^(jm) u_ (ij)^(ij) u_ (km)^(km) + 1/2 λ g_ (ij)^(ij) g_ (km)^(km) u_ (ij)^(ij) u_ (km)^(km)

The (generalized) trace describes a pure Hydrostatic compression

In[39]:=

1/3 guu[red[i], red[j]] gdd[red[p], red[q]] uuu[red[p], red[q]]

Out[39]=

1/3 g_ (pq)^(pq) g_ (ij)^(ij) u_ (pq)^(pq)

The traceless part of the strain tensor describes a pure shear stress

In[40]:=

uuu[red[i], red[j]] - 1/3 guu[red[i], red[j]] gdd[red[p], red[q]] uuu[red[p], red[q]]

Out[40]=

u_ (ij)^(ij) - 1/3 g_ (pq)^(pq) g_ (ij)^(ij) u_ (pq)^(pq)

Hooke's Law

Hooke's Law is obtained from the free energy expansion for isotropic materials (in terms of μ and the Hydrostatic modulus K) applied to the equilibrium equation

In[41]:=

Out[41]=

The stress tensor can be found taking the partial derivative of the Free energy respect to the strain tensor σ_ (r s)^(r s)=∂F/∂u_ (r s)^(r s).

The metric is independent of the displacement

In[42]:=

HoldPattern[PartialD[_][Tensor[g, _, _], Tensor[u, _, _]]] = 0 ;

Metric simplifying the entire expression and using TensorSimplify.

In[43]:=

Expand[PartialD[{u, δ, g, Γ}][Tensor[F], uuu[red[r], red[s]]]/.FreeEnergy[μ, K]]//KroneckerAbsorb[δ]//MetricSimplify[g]

DisplayForm[FrameBox[IsotropicStressTensorRule = σdd[red[r_], red[s_]] →Collect[SimplifyTensorSum[TensorSimplify[%]/. {gdu[b_, b_] →3, gud[b_, b_] →3}], {μ}]]]

Out[43]=

Out[44]//DisplayForm=

σ_ (r_s_)^(r_s_) →K g_ (rs)^(rs) u_ (bb)^(bb) + μ (2 u_ (rs)^(rs) - 2/3 g_ (rs)^(rs) u_ (bb)^(bb))

Equilibrium  Equation

The equilibrium equation that comes from Newtons's Law is    σ_ (i j)^(i j) _ (; j)f^i==0.

Neglecting the volume forces we have only the first term.  It is convenient to expand in terms of the covariant indices of the stress tensor, so the equilibrium equation becomes

In[45]:=

equilibrium = guu[j, k] CovariantD[σdd[i, j], k] == zerod[i]//ToFlavor[red]

Out[45]=

σ_ (ij)^(ij) _ (; k) g_ (jk)^(jk) == zero_i^i

The zero tensor balances the free index on both sides of the equation, behaves as a zero, and expands to a zero array. Expanding and expressing the stress tensor in terms of the strain tensor as first step

In[46]:=

equilibrium

%//ExpandCovariantD[labs, red @ a]

%/.IsotropicStressTensorRule ;

step1 = %//SumExpansion[red @ k, red @ j, red @ a]//ArrayExpansion[red @ i]

Out[46]=

σ_ (ij)^(ij) _ (; k) g_ (jk)^(jk) == zero_i^i

Out[47]=

g_ (jk)^(jk) (-Γ_ (aki)^(aki) σ_ (aj)^(aj) - Γ_ (akj)^(akj) σ_ (ia)^(ia) + ∂σ_ (ij)^(ij)/∂x_k^k) == zero_i^i

Out[49]=

and finally we obtain the homogeneous equation for the displacement in cylindrical coordinates

In[50]:=

step2 = step1/.DisplacementDivergenceRule/.TensorValueRules[u]//ExpandAll

Out[50]=

Example1: Cylindrical tube under an internal constant  hydrostatic pressure

This is an example solved in Landau's book Elasticity Theory. Chapter 7 problem 4

Using the expression of the deformation in cylindrical coordinates. The deformation depends only on the radial coordinate so we define the rule

In[51]:=

ConditionExample1 = {PartialD[_][_, xu[i : red[2] | red[3]]] →0, PartialD[_][_, HoldPattern[{___, Tensor[x, {red[2] | red[3]}, {Void}], ___}]] →0}

Out[51]=

{PartialD[_][_, x_ (i : 2 | 3)^(i : 2 | 3)] →0, PartialD[_][_, HoldPattern[{___, x_ (2 | 3)^(2 | 3), ___}]] →0}

which simplifies the equilibrium equations to

In[52]:=

step3 = Simplify[step2/.ConditionExample1]

Out[52]=

Paying attention to the first equilibrium equation we realize that it is an ODE (Ordinary Differential Equation) that can be easily solved

In[53]:=

First[step3]

%/.{ud[red @ 1] →U[xu[red @ 1]]}/.xu[red @ 1] →ρ

DSolve[%, U[ρ], ρ] _[[1, 1]]

Out[53]=

((3 K + 4 μ) (-u_1^1 + x_1^1 (x_1^1 ∂^2u_1^1/∂x_1^1∂x_1^1 + ∂u_1^1/∂x_1^1)))/x_1^1 == 0

Out[54]=

((3 K + 4 μ) (-U[ρ] + ρ (U^′[ρ] + ρ U^′′[ρ])))/ρ == 0

Out[55]=

U[ρ] → ((1 + ρ^2) C[1])/(2 ρ) + ( (-1 + ρ^2) C[2])/(2 ρ)

Out[56]=

{u_1^1→C[3]/x_1^1 + C[4] x_1^1, u_2^2→0, u_3^3→0}

Alternatively, we can use the new feature of Tensorial 4.0, a copy paste , to solve the differential equation of the first equilibrium equation

In[57]:=

-u_1^1 + x_1^1 (x_1^1 ∂^2u_1^1/∂x_1^1∂x_1^1 + ∂u_1^1/∂x_1^1) == 0

%/.{ud[red @ 1] →U[xu[red @ 1]]}/.xu[red @ 1] →ρ

DSolve[%, U[ρ], ρ] _[[1, 1]]

Out[57]=

-u_1^1 + x_1^1 (x_1^1 ∂^2u_1^1/∂x_1^1∂x_1^1 + ∂u_1^1/∂x_1^1) == 0

Out[58]=

-U[ρ] + ρ (U^′[ρ] + ρ U^′′[ρ]) == 0

Out[59]=

U[ρ] → ((1 + ρ^2) C[1])/(2 ρ) + ( (-1 + ρ^2) C[2])/(2 ρ)

Out[60]=

{u_1^1→C[3]/x_1^1 + C[4] x_1^1, u_2^2→0, u_3^3→0}

The remaining two equilibrium equations are zero as it shown in the next subsection.

Curl of the Curl of the Displacement

In[61]:=

HoldPattern[PartialD[_][Tensor[ε, _, _], _]] = 0

Out[61]=

0

The curl of the displacement tensor is

In[62]:=

(gdd[red[l], red[k]] εuuu[red[i], red[j], red[k]] CovariantD[ud[red[i]], red[j]])/Tensor[g^(1/2)]

ExpandCovariantD[labs, red[a]][%] <br />

(ArrayExpansion[red[l]][SumExpansion[red[i], red[j], red[a], red[k]][%]]/.ConditionExample1)/.Tensor[g^(1/2)] :→ xu[red[1]]

SetTensorValues[Curlud[red[p]], %]

Out[62]=

(u_i^i_ (; j) g_ (lk)^(lk) ε_ (ijk)^(ijk))/g^(1/2)

Out[63]=

(g_ (lk)^(lk) ε_ (ijk)^(ijk) (-u_a^a Γ_ (aji)^(aji) + ∂u_i^i/∂x_j^j))/g^(1/2)

Out[64]=

{0, x_1^1 ∂u_3^3/∂x_1^1, -∂u_2^2/∂x_1^1/x_1^1}

Defining a special display format for Curlu, and taking the curls once more on top of our previous result:

In[66]:=

TensorLabelFormat[Curlu, "(∇×u)"]

Taking the curls once more on the top of our previous result

In[67]:=

(gdd[red[l], red[k]] εuuu[red[i], red[j], red[k]] CovariantD[Curlud[red[i]], red[j]])/Tensor[g^(1/2)]

ExpandCovariantD[labs, red[a]][%]//Expand

Expand[ArrayExpansion[red[l]][SumExpansion[red[i], red[j], red[a], red[k]][-%]]/.TensorValueRules[Curlu]/.Tensor[g^(1/2)] →xu[red[1]]]

%/.ConditionExample1

Out[67]=

((∇×u)_i^i_ (; j) g_ (lk)^(lk) ε_ (ijk)^(ijk))/g^(1/2)

Out[68]=

-((∇×u)_a^a g_ (lk)^(lk) Γ_ (aji)^(aji) ε_ (ijk)^(ijk))/g^(1/2) + (g_ (lk)^(lk) ε_ (ijk)^(ijk) ∂(∇×u)_i^i/∂x_j^j)/g^(1/2)

Out[69]=

Out[70]=

{0, ∂^2u_2^2/∂x_1^1∂x_1^1 - ∂u_2^2/∂x_1^1/x_1^1, ∂^2u_3^3/∂x_1^1∂x_1^1 + ∂u_3^3/∂x_1^1/x_1^1}

which are proportional to the remaining two equilibrium equations.  The curl is interpreted as the measure of an infinitesimal rotation. A static cylinder under hydrostatic pressure does not generate rotations, so the last two  equilibrium equations are identical to zero.

The conditions of a relative pressure p inside the tube defines the values of the constants

In[71]:=

ConstantsC4C3 = {C[4] → (p r^2 (1 + σ) (1 - 2 σ))/((R^2 - r^2) "E"), C[3] → (p r^2 R^2 (1 + σ))/((R^2 - r^2) "E")}

Out[71]=

{C[4] → (p r^2 (1 - 2 σ) (1 + σ))/(E (-r^2 + R^2)), C[3] → (p r^2 R^2 (1 + σ))/(E (-r^2 + R^2))}

where E and σ are the Young and Poisson moduli that are easier to measure in the experiment

Calculation of the Strain tensor

In[72]:=

ArrayExpansion[red[i], red[j]][udd[red[i], red[j]]]/.TensorValueRules[u]

Out[72]=

In our case the Strain tensor is

In[73]:=

ArrayExpansion[red[i], red[j]][udd[red[i], red[j]]]/.TensorValueRules[u]/.SolDisplacement/.ConstantsC4C3

SetTensorValueRules[TubeStraindd[red[i], red[j]], %]

Out[73]=

The divergence of the displacement tensor (trace of the strain tensor) is

In[75]:=

udu[red[k], red[k]]

TraceTubeStrainRule = udu[red[b_], red[b_Symbol]] →%/.DisplacementDivergenceRule/.SolDisplacement/.ConstantsC4C3

Out[75]=

u_ (kk)^(kk)

Out[76]=

The following rule defines the Lame moduli in terms of the Young Poisson moduli.

In[77]:=

ToPoissonYoung = Simplify[Solve[{"E" == (9 K μ)/(3 K + μ), 2 σ == (3 K - 2 μ)/(3 K + μ)}, {K, μ}][[1]]]

Out[77]=

{K→E/(3 - 6 σ), μ→E/(2 + 2 σ)}

Calculating the Stress Tensor

From the previous parts we get everything necessary to calculate the Stress Tensor

In[78]:=

σdd[red[i], red[j]]

%/.IsotropicStressTensorRule/.TraceTubeStrainRule

TubeStressTensor = Simplify[ArrayExpansion[red[i], red[j]][%]/.u→TubeStrain/.TensorValueRules[TubeStrain]/.ToPoissonYoung]

Out[78]=

σ_ (ij)^(ij)

Out[79]=

Out[80]=

{{(p r^2 (R^2 - (x_1^1)^2))/((r^2 - R^2) (x_1^1)^2), 0, 0}, {0, -(p r^2 (R^2 + (x_1^1)^2))/(r^2 - R^2), 0}, {0, 0, -(2 p r^2 σ)/(r^2 - R^2)}}

which agrees with Landau's solution after  σ_ (2 2)^(2 2) is divided by the square of x_ 1^1. This is because Landau uses the physical components instead of the covariant components as we use here.


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