Corpo affine elastico vincolato

v. 2.13 (9/7/2004) © A.Tatone [Università dell'Aquila]

Id

nbid=NotebookInformation[SelectedNotebook[]]

nbid[[1, 2, 2]]

$System

$Version

Date[]

nbdir = 
      nbid[[1, 2, 1]] // 
        ToFileName

Directory[]

SetDirectory[nbdir]

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]]])

Tensione

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

Origine delle coordinate e polo

Coordinate del centro del parallelepipedo  (fissare queste coordinate 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 - (L2/2) e2 - (L1/2) e1

Volume del solido

vol := L1 L2 L3 /2;

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

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

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

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

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

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

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

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

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

spigolo3p1p = 
      xC + (L3/2) e3 + (L1/2) e1  +  ζ2 e2

spigolo3p1m = 
      xC + (L3/2) e3 - (L1/2) e1  +  ζ2 e2

spigolo3m1p = 
      xC - (L3/2) e3 + (L1/2) e1  +  ζ2 e2

spigolo3m1m = 
      xC - (L3/2) e3 - (L1/2) e1  +  ζ2 e2

Sinonimi

spigolo2p1p := spigolo1p2p; spigolo2m1p := spigolo1p2m; 
    spigolo2p1m := spigolo1m2p; 
    spigolo2m1m := spigolo1m2m;
spigolo3p2p := spigolo2p3p; spigolo3m2p := spigolo2p3m; 
    spigolo3p2m := spigolo2m3p; 
    spigolo3m2m := spigolo2m3m;
spigolo1p3p := spigolo3p1p; spigolo1m3p := spigolo3p1m; 
    spigolo1p3m := spigolo3m1p; 
    spigolo1m3m := spigolo3m1m;

Contorno della forma del corpo

Tutto il parallelepipedo

corpo = {spigolo1p2p, spigolo1p2m, spigolo1m2p, 
          spigolo1m2m, 
spigolo2p3p, spigolo2p3m, 
          spigolo2m3p, spigolo2m3m, 
spigolo3p1p, 
          spigolo3p1m, spigolo3m1p, 
          spigolo3m1m} /. {ζ1 -> ξ L1, ζ2 -> ξ L2, ζ3 -> ξ L3}

Altri spigoli

spigolo3pd := xC + (L3/2) e3 + (-e1 + e2) ξ

spigolo3md := xC - (L3/2) e3 + (-e1 + e2) ξ

Parte del parallelepipedo

corpo = {spigolo3pd, spigolo3md, spigolo1p2m, spigolo1m2p, 
          spigolo1m2m, 
spigolo2m3p, 
          spigolo2m3m, 
spigolo3p1m, 
          spigolo3m1m} /. {ζ1 -> ξ L1, ζ2 -> ξ L2, ζ3 -> ξ L3}

Vincoli

Inizializzazione

vincoli = {};

Spigolo a destra in basso

n = -Sin[(π/4)] e1 + Cos[(π/4)] e2;

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

vincoli = Join[vincoli, vincolo1]

Spigolo a sinistra in basso

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

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

vincoli = Join[vincoli, vincolo2, vincolo3]

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

b := p e1

mMa = (∫_(-(L2/2))^((L2/2)))(∫_(-(L3/2))^((L3/2)))(prt[(faccia1m - x0), 
                  b]) \[DifferentialD]ζ2 \[DifferentialD]ζ3;

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

fa = (∫_(-(L1/2))^((L1/2)))(∫_(-(L3/2))^((L3/2)))b \[DifferentialD]ζ1 \[DifferentialD]ζ3

Tensione

wpar = Join[Flatten[mG], w0]

MatrixForm[mT]

prs[fa, w0] //. wvinc

prs[fa, w0]

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

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

Deformazione infinitesima

Dilatazione infinitesima (dalla funzione di risposta)

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

MatrixForm[mE // Factor]

Condizioni di vincolo

MatrixForm[mH]

MatrixForm[mH /. uvinc]

u0 /. uvinc

MatrixForm[mθ]

MatrixForm[mE + mθ]

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

Descrittori dei campi di spostamento affini vincolati e componenti dalla tensione non determinate dalle equazioni di bilancio

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]]])

Soluzione

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]

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

Visualizzazione

Eliminazione spostamento rigido

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

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

mH /. uvinc /. sol // Simplify

Complement[{u[x0], mH} /. uvinc /. sol // Flatten // 
          Simplify // Variables, {L1, L2, L3, ℒ, p, 
        q, λ, μ}]

urig = Map[# -> 0  &, %]

Funzioni di visualizzazione

projection[n_] := 
      Drop[If[n == 3 || 
              n == 
                2, #, {#[[1]], -#[[3]], #[[2]]}], {n}] &

proj[n_] := 
        Block[{q = 1, p = 1, L1 = 1, L2 = 1, 
            L3 = 1, ℒ = 1, ν = 0.3, Y = 20}, 
          Block[{x0Q = projection[n] /@ corpo, 
              xQ = Simplify[projection[n] /@ ϕ /@  
                              corpo /. uvinc /. sol /. urig /. lamé]}, 
            plt0 = ParametricPlot[
                Evaluate[x0Q], {ξ, -(1/2), (1/2)}, 
                PlotStyle -> {Hue[0.2]}, AspectRatio -> Automatic, 
                Axes -> False, DisplayFunction -> Identity]; 
            plt1 = ParametricPlot[Evaluate[xQ], {ξ, -(1/2), (1/2)}, 
                PlotStyle -> {Hue[0.9]}, AspectRatio -> Automatic, 
                Axes -> False, DisplayFunction -> Identity]]; 
          Show[plt0, plt1, 
            DisplayFunction -> $DisplayFunction]];

wframe := 
        Block[{q = 1, p = 1, L1 = 1, L2 = 1, 
            L3 = 1, ℒ = 1, ν = 0.3, Y = 20}, 
          Block[{x0Q = corpo, 
              xQ = Simplify[ϕ /@  
                            corpo /. uvinc /. sol /. urig /. lamé]}, 
            plt0 = ParametricPlot3D[
                Evaluate[x0Q], {ξ, -(1/2), (1/2)}, 
                AspectRatio -> Automatic, Axes -> False, 
                Boxed -> False, DisplayFunction -> Identity]; 
            plt1 = ParametricPlot3D[
                Evaluate[xQ], {ξ, -(1/2), (1/2)}, 
                AspectRatio -> Automatic, Axes -> False, 
                Boxed -> False, DisplayFunction -> Identity]]; 
          Show[plt0, plt1, 
            DisplayFunction -> $DisplayFunction]];

wframe;

Sezione ortogonale a e3

proj[3];

Sezione ortogonale a e2

proj[2];

Sezione ortogonale a e1

proj[1];

Momento delle forze reattive

Da moltiplicare per il volume

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

fa