Corpo affine elastico vincolato

v. 2.08  (23/6/2003)   ©  A. Tatone  [Università dell'Aquila]

Inizializzazione

Off[General::"spell"];
Off[General::"spell1"];
Off[Solve::"svars"];

Prodotto tensoriale

prt[u_, v_] := Transpose[{v}] . {u}

Prodotto scalare

prs[u_, v_] := Flatten[u] . Flatten[v]

Traccia

tr[m_] := 
      (Subscript[m, [[1, 1]]]) + 
        (Subscript[m, [[2, 2]]]) + 
        (Subscript[m, [[3, 3]]])

Gradiente dello spostamento

mH = Array[ug[#1, #2]  &, {3, 3}];

MatrixForm[mH]

Rotazione infinitesima

mθ=(({{0, -θ[3], θ[2]}, {θ[3], 0, -θ[1]}, {-θ[2], +θ[1], 0}}));

Spostamento

u[{ζ1_, ζ2_, ζ3_}] := 
      u0 + mH . ({ζ1, ζ2, ζ3} - x0)

u0 := {u01, u02, u03}

Deformazione

ϕ[{ζ1_, ζ2_, ζ3_}] := {ζ1, ζ2, ζ3} + u[{ζ1, ζ2, ζ3}]

Gradiente dell'atto di moto

mG = Array[g[#1, #2]  &, {3, 3}];

Atto di moto affine

w[{ζ1_, ζ2_, ζ3_}] := 
      w0 + mG . ({ζ1, ζ2, ζ3} - x0)

w0 := {w01, w02, w03}

Vettori base

e1 = {1, 0, 0};
e2 = {0, 1, 0};
e3 = {0, 0, 1};

Matrice della identità

mI = {e1, e2, e3}

Relazioni tra i moduli elastici

lamé = (Subscript[(Solve[{Y == ((3  λ + 2  μ)/(λ + μ)) μ, ν == (λ/(2 (λ + μ)))}, {λ, μ}]), [[1]]])

Volume del parallelepipedo

vol := L1 L2 L3;

Origine delle coordinate e polo

Coordinate del centro del parallelepipedo (da fissare) (equivale alla scelta dell'origine)

xC = {0, 0, 0}

Polo (da scegliere) (è meglio usare una espressione del tipo   xC + ...  poichè risulta indipendente dalla scelta dell'origine)

x0 := xC - (L1/2) e1 - (L2/2) e2

Lunghezze spigoli

Eventuali valori o relazioni

Parametrizzazione delle facce

faccia1m = 
      xC - (L1/2) e1  +  ζ3 e3 + ζ2 e2

faccia1p = 
      xC + (L1/2) e1  +  ζ3 e3 + ζ2 e2

faccia3m = 
      xC - (L3/2) e3  +  ζ1 e1 + ζ2 e2

faccia3p = 
      xC + (L3/2) e3  +  ζ1 e1 + ζ2 e2

faccia2m = 
      xC - (L2/2) e2  +  ζ1 e1 + ζ3 e3

faccia2p = 
      xC + (L2/2) e2  +  ζ1 e1 + ζ3 e3

Parametrizzazione degli spigoli

spigolo1p2m = 
      xC - (L2/2) e2 + (L1/2) e1  +  ζ3 e3

spigolo1p2p = 
      xC + (L2/2) e2 + (L1/2) e1  +  ζ3 e3

spigolo1m2p = 
      xC - (L1/2) e1 + (L2/2) e2 + ζ3 e3

spigolo1m2m = 
      xC - (L1/2) e1 - (L2/2) e2 + ζ3 e3

Vincoli

Inizializzazione

vincoli = {};

Spigolo a destra

vincolo1 = 
      Complement[
          CoefficientList[(u[spigolo1p2m] . (e1) ) // 
                Simplify, {ζ1, ζ2, ζ3}] // Flatten, {0}] // 
        Union

vincoli = Join[vincoli, vincolo1]

Faccia inferiore

vincolo2 = 
      Complement[
          CoefficientList[(u[faccia2m] . (e2) ) // 
                Simplify, {ζ1, ζ2, ζ3}] // Flatten, {0}] // 
        Union

vincoli = Join[vincoli, vincolo2]

Lista dei vincoli su spostamenti e atti di moto

vincoli

Length[vincoli]

u0

Flatten[Join[u0, mH]]

Solve[Map[(# == 0) &, vincoli], Flatten[Join[mH, u0]]]

uvinc = %[[1]]

Length[uvinc]

wvinc = uvinc /. ug -> g /. {u01 -> w01, u02 -> w02, 
          u03 -> w03}

mH /. uvinc // MatrixForm

u0 /. uvinc

mG //. wvinc // MatrixForm

w0 //. wvinc

Risultante e momento delle forze attive

mMa = (∫_(-(L3/2))^((L3/2)))(prt[(spigolo1m2p - 
                      x0), q e1]) \[DifferentialD]ζ3 + (∫_(-(L3/2))^((L3/2)))(∫_(-(L1/2))^((L1/2)))(prt[(faccia2p - 
                        x0), -p e2]) \[DifferentialD]ζ1 \[DifferentialD]ζ3;

(1/vol) mMa // Simplify // MatrixForm

fa = (∫_(-(L3/2))^((L3/2)))(q e1) \[DifferentialD]ζ3 + (∫_(-(L3/2))^((L3/2)))(∫_(-(L1/2))^((L1/2)))(-p e2) \[DifferentialD]ζ1 \[DifferentialD]ζ3

Tensione

wpar = Join[Flatten[mG], w0]

mT=(({{σ[1, 1], σ[1, 2], σ[1, 3]}, {σ[1, 2], σ[2, 2], σ[2, 3]}, {σ[1, 3], σ[2, 3], σ[3, 3]}}));

MatrixForm[mT]

prs[fa, w0] //. wvinc

prs[fa, w0]

((prs[fa, w0] + prs[mMa, mG]) (1/vol) - prs[mT, mG]) //. 
      wvinc

Complement[
        Coefficient[((prs[fa, w0] + prs[mMa, mG]) (1/vol) - 
                prs[mT, mG]) //. wvinc, wpar], {0}] // Simplify

tens = (Subscript[(Solve[Map[(# == 0)  &, %], Union[Flatten[mT]]]), [[1]]])

mT /. tens // MatrixForm

Dilatazione infinitesima e rotazione infinitesima

Dalla tensione

mE = (1/(2  μ)) (mT - (λ/(3  λ + 
                          2  μ)) tr[mT] mI) //. tens // 
          Simplify;

MatrixForm[mE // Factor]

Dalle condizioni di vincolo

MatrixForm[mH]

MatrixForm[mH /. uvinc]

u0 /. uvinc

MatrixForm[mθ]

MatrixForm[mE + mθ]

MatrixForm[(mE + mθ - mH) /. uvinc]

var = 
      Intersection[Variables[(mE + mθ - mH) /. uvinc], 
        Join[Flatten[mH], Flatten[mE], Flatten[mθ], Flatten[mT], 
          u0]]

sol = (Subscript[(Solve[Map[(# == 0) &, Flatten[(mE + mθ - mH)] /. uvinc], var]), [[1]]])

MatrixForm[mH /. uvinc /. sol // Simplify]

MatrixForm[mE /. sol // Simplify]

MatrixForm[mθ /. sol // Simplify]

u0 /. uvinc /. sol /. lamé // Simplify

MatrixForm[mH /. uvinc /. sol /. lamé // 
        Simplify]

MatrixForm[mE /. sol /. lamé // Simplify]

MatrixForm[mT /. tens /. sol /. lamé // 
        Simplify]

Visualizzazione

u[{ζ1, ζ2, ζ3}] /. uvinc /. sol /. 
        lamé // Simplify

u[x0] /. uvinc /. sol /. lamé // Simplify

x0

uval = {u03 -> 0}

Sezione ortogonale a e3

Block[{q = 1, p = 1, L1 = 1, L2 = 1, L3 = 1, ν = 0.3, Y = 20, 
          lato1 = xC - (L2/2) e2 + ξ L1 e1, 
          lato2 = xC + (L1/2) e1 + ξ L2 e2, 
          lato3 = xC + (L2/2) e2 + ξ L1 e1, 
          lato4 = xC - (L1/2) e1 + ξ L2 e2}, 
        Block[{x0Q = (Take[#1, 2] &) /@ {lato1, lato2, lato3, lato4}, 
            xQ = Simplify[(Take[#1, 2] &) /@ {ϕ[lato1], ϕ[
                              lato2], ϕ[lato3], ϕ[
                              lato4]} /. uvinc /. sol /. uval /. lamé]}, 
          ParametricPlot[
            Evaluate[Join[xQ, x0Q]], {ξ, -(1/2), (1/2)}, 
            AspectRatio -> Automatic, 
            Axes -> False]]];

Sezione ortogonale a e2

Block[{q = 1, p = 1, L1 = 1, L2 = 1, L3 = 1, ν = 0.3, Y = 20, 
          u03 = 0, lato1 = xC + (L2/2) e2 + (L3/2) e3 + L1 ξ e1, 
          lato2 = xC + (L2/2) e2 - (L1/2) e1 + L3 ξ e3, 
          lato3 = xC + (L2/2) e2 - (L3/2) e3 + L1 ξ e1, 
          lato4 = xC + (L2/2) e2 + (L1/2) e1 + L3 ξ e3}, 
        Block[{x0Q = (Extract[#1, {{1}, {3}}] &) /@ {lato1, lato2, lato3, 
                  lato4}, 
            xQ = Simplify[(Extract[#1, {{1}, {3}}] &) /@ {ϕ[
                              lato1], ϕ[lato2], ϕ[lato3], ϕ[
                              lato4]} /. uvinc /. sol /. uval /. lamé]}, 
          ParametricPlot[
            Evaluate[Join[xQ, x0Q]], {ξ, -(1/2), (1/2)}, 
            AspectRatio -> Automatic, 
            Axes -> False]]];

Momento delle forze reattive

Da moltiplicare per il volume

mT  - (mMa/vol) /. tens /. sol // MatrixForm

fa