Modelare Matematica in Elasticitate PDF [PDF]

Zitiervorschau

ˇ ˆIN MODELARE MATEMATICA ELASTICITATE Andaluzia Matei

Cuprins 1 Preliminarii 1.1 Tensori cartezieni de ordinul doi . . . . . . . . . . . . . . . . 1.2 Spat¸ii de funct¸ii . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Ecuat¸ii variat¸ionale eliptice; aproximarea solut¸iilor . . . . . .

1 1 6 11

2 Geometria deformat¸iei. Principii generale 2.1 Configurat¸ii. Coordonate. Derivata materialˇa. . . . . . . . . 2.2 Vector deplasare. Tensorul deformat¸iilor infinitezimale. . . . 2.3 Principii generale . . . . . . . . . . . . . . . . . . . . . . . .

15 15 17 32

3 Modelare prin ecuat¸ii cu derivate part¸iale 3.1 Ecuat¸iile Cauchy. Condit¸ii la limitˇa . . . . . . . . . . . . . . 3.2 Legi de comportament . . . . . . . . . . . . . . . . . . . . . 3.3 Modele simplificate . . . . . . . . . . . . . . . . . . . . . . .

39 39 50 62

4 Probleme clasice ˆın elastostaticˇ a 4.1 Comprimarea uniform˘a . . . . . . . . . . . . . . . . . . . . . 4.2 Tract¸iunea simpl˘a . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Torsiunea unui arbore cilindric . . . . . . . . . . . . . . . . .

69 69 76 82

5 Solut¸ii slabe. Aproximarea solut¸iilor slabe 93 5.1 Probleme antiplane . . . . . . . . . . . . . . . . . . . . . . . 93 5.2 Probleme plane . . . . . . . . . . . . . . . . . . . . . . . . . 104 i

5.3 Probleme 3D . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6 Anexˇ a 117 6.1 Coduri MATLAB . . . . . . . . . . . . . . . . . . . . . . . . 117 6.2 Tablouri de date . . . . . . . . . . . . . . . . . . . . . . . . 131

Capitolul 1 Preliminarii

1.1

Tensori cartezieni de ordinul doi

Tensorul cartezian de ordinul doi poate fi definit ca o aplicat¸ie liniar˘a a unui spat¸iu vectorial euclidian ˆın el ˆınsu¸si. Aceastˇa definit¸ie a fost adoptatˇa de exemplu ˆın lucrˇarile [41, 42, 28].1 Fie L un spat¸iu vectorial euclidian ¸si T : L → L un tensor cartezian de ordinul doi. Astfel, pentru fiecare v ∈ L, T v = u ∈ L.

(1.1)

Dac˘a spat¸iul L este raportat la o baz˘a ortonormat˘a {e1 , e2 , ..., en }, 1

Pentru definitii echivalente ale tensorului cartezian de ordinul doi, cˇaci exist˘a mai

multe definit¸ii ale tensorului cartezian ˆın general ¸si ale tensorului cartezian de ordin doi ˆın particular, se pot consulta de exemplu lucrarile [6, 21, 22].

1

atunci orice vector T ej se scrie: T ej = Tij ei .

(1.2)

ˆIn scrierea relat¸iei (1.2) s-a utilizat convent¸ia de sumare a lui Einstein (convent¸ia de sumare dupˇa indicele care se repetˇa); vom utiliza aceastˇa convent¸ie pe parcursul ˆıntregii lucrˇari. Numerele Tij se vor numi componentele tensorului ˆın baza dat˘a. Evident ele sunt ˆın numar de n2 iar (Tij )1≤i,j≤n este o matrice patratic˘a de n × n elemente. Ea caracterizeaz˘a tensorul ˆın baza dat˘a. Utilizˆand (1.1) ¸si (1.2), putem scrie, ui ei = T (vj ej ) = Tij vj ei ¸si de aici, deducem, ui = Tij vj .

(1.3)

Tensorul nul, notat O, este aplicat¸ia care face s˘a corespund˘a oric˘arui vector, vectorul nul 0, Ov = 0 iar tensorul unitate, notat I, este aplicat¸ia ce las˘a neschimbat orice vector v, Iv = v. Suma a doi tensori T ¸si S este un tensor G de componente Gij = Tij + Sij

1 ≤ i, j ≤ n.

Prin compunerea a doi tensori T ¸si S obt¸inem un tensor T S de componente (Tik Skj ) ce corespund produsului matricial; ˆın general, T S ̸= ST . Mult¸imea tensorilor cartezieni de ordinul doi ˆınzestratˇa cu operat¸iile de adunare si compunere anunt¸ate mai sus este un inel unitar. Puterile unui tensor T se definesc dupˇa cum urmeazˇa: T 0 = I, T 1 = T , T 2 = T T , ... .

Un tensor T este inversabil dac˘a existˇa un tensor T −1 numit tensorul invers astfel ˆıncˆat T T −1 = T −1 T = I. Produsul dintre un tensor T ¸si un numˇar real λ este un tensor G de componente Gij = λTij 1 ≤ i, j ≤ n. Mult¸imea tensorilor cartezieni de ordinul doi ˆınzestratˇa cu operat¸iile de adunare si ˆınmult¸ire cu scalari anunt¸ate mai sus este un spat¸iu vectorial real. ˆIn plus, putem defini un produs scalar, notat ” · ”, T · S = tr(T S T ) = tr(ST T ) = Tij Sij , unde T ¸si S sunt tensori cartezieni de ordinul doi. Norma tensorului T notatˇa |T |, este norma indusˇa de produsul scalar, |T | =

√

T · T.

(1.4)

ˆInzestrat cu acest produs scalar, spat¸iul vectorial real al tensorilor cartezieni de ordinul doi este un spat¸iu Hilbert. Tensorul transpus se noteazˇa cu T T ¸si este definit de relat¸ia: u · Tv = TTu · v

∀u, v ∈ L.

Au loc urmˇatoarele proprietˇa¸ti (T 1 + T 2 )T = T T1 + T T2 , (T 1 T 2 )T = T T2 T T1 , (T n )T = (T T )n ;

(T −1 )T = (T T )−1 ;

(T T )T = T .

Un tensor S se nume¸ste simetric dac˘a S = S T . Notˇam cˇa spat¸iul vectorial real al tensorilor simetrici cartezieni de ordinul doi, este un spat¸iu Hilbert. Un tensor Ω se nume¸ste antisimetric dac˘a Ω = −ΩT .

Orice tensor T se descompune ˆın mod unic ˆın suma a doi tensori, unul simetric S ¸si altul antisimetric Ω. ˆIntr-adevˇar, T =S+Ω unde

1 S = (T + T T ); Ω= 2 Un tensor R se nume¸ste ortogonal dac˘a invariant produsul scalar, adicˇa Ru · Rv = u · v

1 (T − T T ). 2 aplicat¸ia care ˆıl define¸ste las˘a

∀u, v ∈ L.

Se poate demonstra cˇa R este ortogonal dac˘a ¸si numai dac˘a RRT = RT R = I sau, echivalent, dac˘a ¸si numai dac˘a R−1 = RT . Urma tensorului T se noteaz˘a tr T ¸si se define¸ste astfel: tr T = Tii Evident tr T = tr T T . ˆIn plus, se poate verifica prin induct¸ie cˇa tr(T T )n = trT n . Un tensor simetric S se nume¸ste pozitiv definit dac˘a u · Su > 0 ∀u ̸= 0. Un tensor simetric S se nume¸ste pozitiv semidefinit dac˘a u · Su ≥ 0 ∀u ∈ L.

Fie T un tensor cartezian de ordinul doi ¸si fie ecuat¸ia: T u = λu. Un vector u (nenul) ce satisface aceast˘a ecuat¸ie se nume¸ste vector propriu iar scalarul λ pentru care se verificˇa ecuat¸ia se nume¸ste valoare proprie. ˆIntr-o baz˘a dat˘a, (Tij − λδij )uj = 0

i ∈ {1, 2, ..., n},

unde δij este simbolul lui Kronecker,     1 dacˇa δij =    0 dacˇa

i = j; i ̸= j.

Ca acest sistem s˘a admit˘a solut¸ii nebanale este necesar ¸si suficient s˘a avem: det(Tij − λδij ) = 0, unde Pn (λ) = det(T − λI) se nume¸ste polinom caracteristic iar Pn (λ) = 0

(1.5)

se nume¸ste ecuat¸ie caracteristic˘a. Valorile proprii sunt reale dac˘a ¸si numai dac˘a tensorul T este simetric; valorile proprii sunt pozitive dac˘a ¸si numai dac˘a tensorul T este pozitiv definit. Dac˘a T este simetric ¸si dac˘a r˘ad˘acinile ecuat¸iei caracteristice sunt distincte, exist˘a n vectori proprii ortogonali doi cˆate doi (ce formeaz˘a o baz˘a); dac˘a r˘ad˘acinile sunt multiple se pot determina de asemenea n vectori proprii ortogonali doi cˆate doi. Baza determinat˘a de vectorii proprii se nume¸ste baz˘a proprie. Numim invariant al tensorului de ordin doi orice aplicat¸ie definit˘a pe mult¸imea acestor tensori cu valori scalare. Fie I1 = trT = Tii

I2 = trT 2 = Tij Tji In = trT n = ... Ace¸sti invariant¸i se numesc invariant¸i principali ai tensorului T . Cu ajutorul lor putem scrie Pn (λ) = det(T − λI) = λn − I1 λn−1 + I2 λn−2 − ... ± In . Un tensor de forma αI, unde α ∈ R, se nume¸ste tensor sferic. Un tensor simetric avˆand urma zero se nume¸ste tensor deviator. Orice tensor simetric de ordinul doi S se reprezint˘a ˆın mod unic ca suma dintre un tensor sferic ¸si un tensor deviator. ˆIntr-adev˘ar, S=

( trS trS ) I+ S− I . n n

Dat fiind un tensor de ordin doi T , numim deviatorul lui T tensorul T D definit dupˇa cum urmeazˇa T D := T −

trT I. n

(1.6)

ˆIn paginile prezentei lucrˇari, vom nota cu Sd , d ∈ {2, 3}, spat¸iul tensorilor simetrici cartezieni de ordinul doi, definit¸i pe Rd .

1.2

Spat¸ii de funct¸ii

Spat¸ii de funct¸ii scalare Fie Ω ⊂ R2 un domeniu deschis, mˇarginit, conex, cu frontierˇa netedˇa. ˆIncepem prin a reaminti pe scurt cˆateva rezultate asupra spat¸iilor Sobolev, H 1 (Ω) = {u ∈ L2 (Ω) | ∂i u ∈ L2 (Ω), i = 1, 2}, unde ∂i u pentru i ∈ 1, 2 desemneazˇa derivata ˆın sens slab; pentru definit¸ia derivatei ˆın sens slab a se consulta de exemplu [7].

H 1 (Ω) este un spat¸iu Hilbert ˆınzestrat cu produsul scalar (u, v)H 1 (Ω) = (u, v)L2 (Ω) + (∂i u, ∂i v)L2 (Ω) ¸si norma asociatˇa 1

∥u∥H 1 (Ω) = (u, u)H2 1 (Ω) . Amintim cˇa: • C ∞ (Ω) este dens ˆın H 1 (Ω); • (Teorema lui Rellich ) H 1 (Ω) ⊂ L2 (Ω) inject¸ie compactˇa; • (Teorema lui Sobolev ) Existˇa o aplicat¸ie liniarˇa ¸si continuˇa γ : H 1 (Ω) → L2 (Γ) astfel ˆıncˆat γ u = u |Γ ∀u ∈ C 1 (Ω). Aplicat¸ia γ se nume¸ste aplicat¸ia urmˇa; ea este definitˇa ca prelungirea prin densitate a aplicat¸iei u → u|Γ definitˇa pentru u ∈ C 1 (Ω.) Imaginea lui 1 H 1 (Ω) prin aceastˇa aplicat¸ie, notatˇa H 2 (Γ), este un subspat¸iu al lui L2 (Γ) care este Hilbert ˆın raport cu structura transportatˇa de γ. ˆIn plus, notˇam cˇa • γ : H 1 (Ω) → L2 (Γ) nu este aplicat¸ie surjectivˇa ¸si nici injectivˇa; 1

• H 2 (Γ) ⊂ L2 (Γ) inject¸ie continuˇa; 1

• γ : H 1 (Ω) → H 2 (Γ) este o aplicat¸ie liniarˇa, continuˇa ¸si surjectivˇa. • γ : H 1 (Ω) → L2 (Γ) este un operator compact. Pentru detalii, recomandˇam de exemplu [1]. Introducem subspat¸iul ˆınchis al lui H 1 (Ω) definit prin V = {v ∈ H 1 (Ω) | γv = 0 a.p.t. pe Γ1 }.

(1.7)

Ca subspat¸iu ˆınchis al unui spat¸iu Hilbert, spat¸iul V este un spat¸iu Hilbert ˆınzestrat cu produsul scalar (·, ·)H 1 (Ω) . Amintim urmˇatoarea inegalitate de tip Poincar´e.

Teorema 1.1. Fie Γ1 o port¸iune a frontierei lui Ω astfel ˆıncˆ at m(Γ1 ) > 0. Atunci existˇa o constantˇ a CP = CP (Ω, Γ1 ) > 0 astfel ˆıncˆ at ∥u∥L2 (Ω) ≤ CP ∥∇u∥L2 (Ω)2

∀u ∈ V.

(1.8)

Pentru demonstrat¸ia acestui enunt¸ recomandˇam [31]. ∫ Notˇam cˇa aplicat¸ia (u, v) 7→ Ω ∇u · ∇v dx este un produs scalar pe V care induce norma ∥∇u∥L2 (Ω)2 . Utilizˆand Teorema 1.1 se poate demonstra cˇa norma ∥∇u∥L2 (Ω)2 este echivalentˇa cu norma ∥u∥H 1 (Ω) . Prin urmare, V este un spat¸iu Hilbert real, ˆınzestrat cu produsul scalar ∫ ∇u · ∇v dx

(u, v)V =

∀u, v ∈ V

(1.9)

Ω

¸si norma asociatˇa ∥u∥V = ∥∇ u∥L2 (Ω)2

∀u ∈ V.

(1.10)

De asemenea, notˇam cˇa existˇa o constantˇa C0 = C0 (Ω, Γ1 ) > 0 astfel ˆıncˆat ∥v∥L2 (Γ) ≤ C0 ∥v∥V ,

∀v ∈ V.

(1.11)

Spat¸ii de funct¸ii vectoriale Fie Ω ⊂ Rd , d ∈ {2, 3} un domeniu deschis, mˇarginit, conex, cu frontierˇa netedˇa. Introducem urmˇatoarele spat¸ii de funct¸ii H = {u = (ui ) | ui ∈ L2 (Ω), i ∈ 1, d}, H = {σ = (σij ) | σij = σji ∈ L2 (Ω), i, j ∈ 1, d}, ∂ui ∂uj + ∈ L2 (Ω), i, j ∈ 1, d} H1 = {u ∈ H | ∂xj ∂xi H1 = {σ ∈ H | σij,j ∈ L2 (Ω), i, j ∈ 1, d}.

Spat¸iile H, H, H1 ¸si H1 sunt spat¸ii reale Hilbert ˆınzestrate cu produsele scalare urmˇatoare ∫ (u, v)H = ui vi dx, (u, v)H1 = (u, v)H + (ε(u), ε(v))H , ∫Ω (σ, τ )H = σij τij dx, (σ, τ )H1 = (σ, τ )H + (Div σ, Div τ )H Ω

unde ε : H1 → H ¸si Div : H1 → H divergent¸ˇa, definit¸i prin ε(u) = (εij (u)),

sunt operatorii deformat¸ie ¸si

1 εij (u) = (ui,j + uj,i ), 2

Div σ = (σij,j ).

Normele pe spat¸iile H, H, H1 ¸si H1 sunt notate prin ∥ · ∥H , ∥ · ∥H , ∥ · ∥H1 ¸si respectiv ∥ · ∥H1 . • C ∞ (Ω)d este dens ˆın H1 ; • H1 ⊂ H inject¸ie compactˇa; • Existˇa o aplicat¸ie liniarˇa ¸si continuˇa γ : H1 → L2 (Γ)d astfel ˆıncˆat γ u = u |Γ pentru oricare u ∈ C 1 (Ω)d . Imaginea lui H1 prin operatorul urmˇa corespunzˇator cazului vectorial se noteazˇa cu HΓ . • γ : H1 → L2Γ nu este aplicat¸ie surjectivˇa ¸si nici injectivˇa; • HΓ ⊂ L2Γ inject¸ie continuˇa; • γ : H1 → HΓ este o aplicat¸ie liniarˇa, continuˇa ¸si surjectivˇa. • γ : H1 → L2Γ este un operator compact. ′

′

Desemnˇam prin HΓ dualul lui HΓ , ¸si (·, ·) produsul de dualitate ˆıntre HΓ ′ ¸si HΓ . Oricare σ ∈ H1 , existˇa un element notat σ ν ∈ HΓ astfel ˆıncˆat (σ ν, γv) = (σ, ε(v))H + (Div σ, v)H

∀v ∈ H1 .

ˆIn plus, dacˇa σ este suficient de regulat (de exemplu de clasˇa C 1 ), avem formula ∫ (σ ν, γv) = σ ν · γv da ∀v ∈ H1 . Γ

Deci, pentru σ suficient de regulat avem urmˇatoarea formulˇ a Green: ∫ (1.12) (σ, ε(v))H + (Div σ, v)H = σ ν · γv da ∀v ∈ H1 . Γ

Definim un subspat¸iu ˆınchis al lui H1 dupˇa cum urmeazˇa: V = { v ∈ H1 | γv = 0

a.p.t. pe

Γ1 }.

(1.13)

Amintim aici inegalitatea lui Korn. Teorema 1.2. [Inegalitatea lui Korn] Existˇa CK = CK (Ω, Γ1 ) > 0 astfel ˆıncˆat ∥ε(v)∥H ≥ CK ∥v∥H1

∀v ∈ V.

(1.14)

O demonstrat¸ie a acestei teoreme poate fi gˇasitˇa de exemplu ˆın [27]. Pe V considerˇam produsul scalar (u, v)V = (ε(u), ε(v))H

∀ u, v ∈ V

(1.15)

¸si fie ∥ · ∥V norma asociatˇa, ∥v∥V = ∥ε(v)∥H

∀ v ∈ V.

(1.16)

Utilizˆand inegalitatea lui Korn rezultˇa cˇa ∥ · ∥H1 ¸si ∥ · ∥V sunt norme echivalente pe V ¸si, prin urmare, (V, (·, ·)V ) este un spat¸iu Hilbert. Notˇam de asemenea cˇa existˇa o constantˇa c0 = c0 (Ω, Γ1 ) > 0 astfel ˆıncˆat ∥v∥L2 (Γ)d ≤ c0 ∥v∥V

∀ v ∈ V.

(1.17)

Pentru detalii privind aceastˇa sect¸iune recomandˇam lucrarile [35, 36].

1.3

Ecuat¸ii variat¸ionale eliptice; aproximarea solut¸iilor

Fie X un spat¸iu Hilbert ˆınzestrat cu produsul scalar (·, ·)X ¸si norma asociatˇa ∥ · ∥X . Fie de asemenea a : X × X → R o formˇa biliniarˇa astfel ˆıncˆat: i) existˇa M > 0 : |a(u, v)| ≤ M ∥u∥X ∥v∥X ii) existˇa m > 0 : a(u, u) ≥

m∥u∥2X

∀u, v ∈ X;

∀u ∈ X.

(1.18) (1.19)

Consider˘am urmˇatoarea problem˘a variat¸ional˘a elipticˇa. Problema 1.3. Sˇa se determine u ∈ X astfel ˆıncˆat, dat fiind L un element ˆın X, a(u, v) = (L, v)X , ∀v ∈ X. (1.20) Urmˇatorul enunt¸ furnizeazˇa existent¸a ¸si unicitatea solut¸iei Problemei 1.3. Teorema 1.4. [Lax-Milgram] Fie a : X × X → R o form˘a biliniar˘a ce verificˇ a (1.18) ¸si (1.19). Dat L ∈ X, exist˘ a un unic element u ∈ X astfel ˆıncˆ at (1.20) se verificˇa. Dac˘a a este ˆın plus formˇa simetric˘a, atunci u minimizeaz˘a funct¸ionala J :X→R

1 J(v) = a(v, v) − (L, v)X . 2

Demonstrat¸ia acestei teoreme se gˇase¸ste de exemplu ˆın [7]. Fie Xh subspat¸iu finit dimensional al lui X, de dimensiune N. Considerˇam urm˘atoarea problem˘a variat¸ionalˇa. Problema 1.5. Fie L ∈ X. Sˇa se determine uh ∈ Xh astfel ˆıncˆat a(uh , vh ) = (L, vh )X ,

∀vh ∈ Xh .

Problema 1.5 se nume¸ste problemˇ a discretˇ a asociatˇa Problemei 1.3.

Teorema 1.6. Admitem (1.18) ¸si (1.19). Atunci Problema 1.5 admite solut¸ie unic˘a. Demonstrat¸ie : Acest rezultat se poate obt¸ine printr-o aplicare direct˘a a Teoremei 1.4, deoarece subspat¸iul Xh este un spat¸iu Hilbert, fiind un subspat¸iu finit dimensional (deci ˆınchis) al unui spat¸iu Hilbert. Vom da ˆıns˘a ¸si o alt˘a demonstrat¸ie direct˘a, constructiv˘a, care va fi util˘a ˆın calculul efectiv al lui uh . Fie {η1 , η2 , ..., ηN } ⊂ Xh astfel ˆıncˆat ηj (xk , yk ) = δjk j, k ∈ 1, N . Vom c˘auta solut¸ia uh sub forma: uh =

N ∑

ξj ηj , ξj ∈ R, j ∈ 1, N .

j=1

Problema Ph se reduce la rezolvarea unui sistem liniar de N ecuat¸ii cu N necunoscute N ∑ a(ηj , ηi )ξj = (L, ηi )X , i ∈ 1, N . (1.21) j=1

Matricea sistemului este inversabil˘a deoarece sistemul N ∑

a(ηj , ηi )ξj = 0, i ∈ 1, N

j=1

are doar solut¸ia banal˘a. ˆIntr-adev˘ar, utilizˆand (1.19), obt¸inem 0=

N ∑ i,j=1

Rezult˘a de aici

N N N ∑ ∑ ∑ a(ηj , ηi )ξj ξi = a( ηj ξj , ηj ξj ) ≥ α∥ ηj ξj ∥2 . j=1

N ∑ j=1

j=1

ηj ξj = 0.

j=1

Cum {ηj }1≤j≤N este baz˘a, deducem cˇa ξj = 0, 1 ≤ j ≤ N. A¸sadar, din (1.21) obt¸inem ˆın mod unic coordonatele solut¸iei uh ˆın baza {η1 , η2 , ..., ηN } ⊂ Xh .  ˆIn cazul ˆın care forma bilinar˘a a(·, ·) este ¸si simetric˘a, rezolvarea Problemei 1.5 se reduce la rezolvarea unui sistem liniar a c˘arui matrice este simetric˘a ¸si pozitiv definit˘a. Analiza erorii produs˘a prin aproximarea solut¸iei u a Problemei 1.3 prin solut¸ia uh a Problemei 1.5, se bazeazˇa pe urmˇatorul rezultat. Teorema 1.7. (Lema lui C`ea) Admitem (1.18) ¸si (1.19). Fie u ∈ X unica solut¸ie a Problemei 1.3 ¸si uh ∈ Xh unica solut¸ie a Problemei 1.5. Atunci existˇ a C > 0 astfel ˆıncˆat ∥u − uh ∥X ≤ C inf ∥u − vh ∥X . vh ∈Xh

Demonstrat¸ie :

Pentru cˇa u ∈ X este solut¸ia Problemei 1.3, a(u, vh ) = (L, vh )X

∀vh ∈ Xh ,

¸si, pentru cˇa uh ∈ Xh este solut¸ia Problemei 1.5, a(uh , vh ) = (L, vh )X

∀vh ∈ Xh .

Utilizˆand liniaritatea lui a ˆın primul argument, putem scrie a(u − uh , vh ) = 0 ∀vh ∈ Xh . ˆIn particular, a(u − uh , uh ) = 0.

(1.22)

Deci, a(u − uh , vh − uh ) = 0 ∀vh ∈ Xh . Folosind (1.19), avem m∥u − uh ∥2X ≤ a(u − uh , u − uh ) = a(u − uh , u − vh ) + a(u − uh , vh − uh ) = a(u − uh , u − vh ) ≤ M ∥u − uh ∥X ∥u − vh ∥X . Prin urmare, M ∥u − vh ∥X ∀vh ∈ X m . ¸si de aici, trecˆand la inf, obt¸inem (1.22) cu C = M m ∥u − uh ∥X ≤



Pentru detalii privind aproximarea solut¸iilor problemelor variat¸ionale eliptice cititorul poate consulta de exemplu [17, 26].

Capitolul 2 Geometria deformat¸iei. Principii generale

2.1

Configurat¸ii. Coordonate. Derivata materialˇ a.

Prin mediu continuu ˆınt¸elegem un sistem material M care la fiecare moment umple complet o regiune a spat¸iului (mult¸ime deschisˇa, mˇarginitˇa, conexˇa din R3 ). Fie un mediu continuu ˆın mi¸scare, ce ocupˇa la momentul t = 0 domeniul Ω ⊂ R3 iar la momentul t > 0 domeniul Ωt ⊂ R3 . Domeniul Ω se nume¸ste configurat¸ie de referint¸ˇa sau configurat¸ie nedeformatˇ a a mediului. Pentru fiecare t > 0, Ωt se nume¸ste configurat¸ie curentˇ a sau deformatˇ a a mediului. Fie P o particulˇa a sistemului material M ¸si fie X = (Xi ) vectorul sˇau 15

de pozit¸ie la t = 0, iar x = (xi ) vectorul sˇau de pozit¸ie la un moment arbitrar t > 0, raportˆandu-ne la o bazˇa ortonormatˇa {e1 , e2 , e3 } din R3 . Componentele (Xi ) se numesc coordonate Lagrange sau coordonate materiale; componentele (xi ) se numesc coordonate Euler sau coordonate spat¸iale. Numim mi¸scare a unui mediu continuu o familie de aplicat¸ii (χ(·, t))t , χ(·, t) : Ω → χ(Ω, t) = Ωt x = χ(X, t),

t > 0.

(2.1)

Presupunem cˇa pentru oricare moment t > 0, aplicat¸ia χ(·, t) este o aplicat¸ie continuˇa ¸si bijectivˇa. Inversa sa se va nota cu χ−1 : χ−1 (·, t) : Ωt → Ω X = χ−1 (x, t) t > 0. Fie ϕ o funct¸ie scalarˇa sau vectorialˇa definitˇa pe Ω × R+ , ϕ = ϕ(X, t). Numim derivatˇa materialˇa, derivata ˆın raport cu t, pentru X constant. Vom ˙ nota aceastˇa derivatˇa prin dtd ϕ sau ϕ; d ∂ ˙ ϕ(X, t) = ϕ(X, t) = ϕ(X, t) dt ∂t

X ∈ Ω, t > 0.

ˆIn particular, viteza ¸si accelerat¸ia punctului material P la momentul t > 0 se definesc cu ajutorul derivatelor materiale astfel: v = x˙ =

d χ(X, t), dt

d2 χ(X, t). dt2 Pe componente, formulele (2.2) ¸si (2.3) se scriu astfel: ¨= a = v˙ = x

vi = x˙i =

d χi (X, t), dt

ai = v˙i = x¨i =

d2 χi (X, t). dt2

(2.2) (2.3)

Fie acum φ o funct¸ie scalarˇa definitˇa pe configurat¸ia deformatˇa, φ : Ωt × R+ → R. Pentru a calcula derivata materialˇa a funct¸iei φ, trebuie sˇa luˇam ˆın considerare dependent¸a lui x de X ¸si t. Astfel, utilizˆand (2.1) putem scrie, φ(x, ˙ t) =

∂ ∂ φ(x, t) vi . φ(x, t) + ∂t ∂ xi

(2.4)

Introducˆand notat¸ia ∇x φ pentru gradientul lui φ ˆın raport cu coordonatele spat¸iale, (2.4) devine ∂φ φ˙ = + ∇x φ · v. (2.5) ∂t Notˇam aici cˇa formula (2.5) este valabilˇa ¸si pentru funct¸ii vectoriale. ˆIn particular, din (2.3) ¸si (2.5) deducem pentru accelerat¸ie urmˇatoarea formulˇa vectorialˇa ∂v a= + (∇x v)v, ∂t formulˇa care pe componente se scrie astfel: ai =

∂ vi ∂ vi + vj . ∂t ∂ xj

Numim traiectorie pentru un punct material de coordonate Lagrange (X1 , X2 , X3 ), ˆıntre momentele 0 ¸si t1 , mult¸imea punctelor x = χ(X, t), pentru t ∈ [0, t1 ], T rajt∈[0,t1 ] (X) = {x | x = χ(X, t),

2.2

∀t ∈ [0, t1 ]}.

Vector deplasare. Tensorul deformat¸iilor infinitezimale.

Numim deformat¸ie aplicat¸ia χ(·, t) : Ω → Ωt , definitˇa prin (2.1) pentru t > 0 fixat. Gradientul sˇau ˆın raport cu coordonatele materiale Xi define¸ste

un cˆamp tensorial notat F , F = (Fij ),

Fij =

∂χi ∂Xj

i, j ∈ {1, 2, 3}.

Tensorul F se nume¸ste gradient al deformat¸iei. Pentru t = 0 avem F = I 3 deci J = detF = 1 ¸si cum J nu se anuleazˇa nicˇaieri (consecint¸ˇa a injectivitˇa¸tii aplicat¸iei χ(·, t) pentru orice t > 0) se obt¸ine J = detF > 0,

∀X ∈ Ω, t > 0.

(2.6)

ˆIn cele ce urmeazˇa vom nota prin ∇X ¸si ∇x operatorii gradient ˆın descriere materialˇa, respectiv spat¸ialˇa, i.e., (∇X φ)i =

∂φ , ∂Xi

(∇x φ)i =

∂φ ∂xi

pentru fiecare funct¸ie scalarˇa φ ¸si ∂ψ ∂ψ , (∇x ψ)ij = ∂Xj ∂xj

(∇X ψ)ij =

pentru fiecare funct¸ie vectorialˇa ψ. A¸sadar, F = ∇X χ. ˆIn plus, pentru orice funct¸ie scalarˇa φ(x, t) = φ(χ(X, t), t), putem scrie ∇X φ = (∇x φ)F . Notˇam cˇa aceastˇa formulˇa este de asemenea valabilˇa pentru funct¸iile vectoriale (verificare pe componente). Considerˇam ˆın aceastˇa sect¸iune tensorii C ¸si G definit¸i prin relat¸iile: C = FT F

(Cij = Fki Fkj ),

1 1 G = (C − I 3 ) (Gij = (Fki Fkj − δij )). 2 2 Tensorul C se nume¸ste tensor dilatare iar tensorul G se nume¸ste tensor deformat¸ie.

Tensorul dilatare C este un tensor simetric. Prin urmare, el posedˇa trei valori proprii reale CI , CII , CIII pe care le numim dilatat¸ii principale. Propozit¸ia 2.8. Valorile proprii CI , CII , CIII ale tensorului dilatare C sunt strict pozitive. Demonstrat¸ie : Dacˇa λ este o valoare proprie a tensorului dilatare C ¸si y este vector propriu asociat, avem Cij yj = λ yi ,

∀i ∈ {1, 2, 3}.

De aici rezultˇa cˇa λyi yi = Cij yj yi = Fki Fkj yj yi = Fki yi Fkj yj = (Fki yi )2 ≥ 0. Cum y este un vector nenul, deducem λ ≥ 0. Sˇa presupunem cˇa λ = 0. ˆIn acest caz Fij yj = 0,

∀i ∈ {1, 2, 3}.

Prin urmare det F = 0. Contradict¸ie cu (2.6).



Tensorul deformat¸ie G este simetric. Vectorii sˇai proprii se numesc direct¸ii principale de deformat¸ie iar valorile proprii se numesc deformat¸ii normale principale. Sˇa notˇam cu XJ , J ∈ {I, II, III} deformat¸iile normale principale. Ele depind de valorile proprii ale tensorului dilatare C dupˇa cum urmeazˇa 1 1 XJ = (CJ − 1) > − , 2 2

J ∈ {I, II, III}.

Semnificat¸ia fizicˇ a a tensorului dilatare C Cunoa¸sterea tensorului dilatare C permite calcularea lungimii vectorului dx (vector infinitezimal ˆın configurat¸ia deformatˇa, corespunzˇator vectorului infinitezimal dX din configurat¸ia init¸ialˇa). ˆIntr-adevˇar, dx = F dX

(dxi = Fij dXj ),

iar de aici |dx|2 = dx · dx = Fki Fkj dXi dXj . Prin urmare, |dx|2 = Cij dXi dXj .

Semnificat¸ia fizicˇ a a tensorului deformat¸ie G Cunoa¸sterea tensorului deformat¸ie G permite calcularea variat¸iei pˇatratului distant¸elor infinitezimale. ˆIntr-adevˇar, |dx|2 − |dX|2 = 2Gij dXi dXj .

(2.7)

Dacˇa mediul continuu se comportˇa ca un corp rigid (adicˇa ˆın timpul deformat¸iei distant¸ele relative dintre puncte sunt constante), din (2.7) deducem cˇa G = O 3 ; mai mult, C = I 3 . Reciproc, dacˇa C = I 3 , sau echivalent G = O 3 , atunci mediul continuu se comportˇa ca un rigid.

Vector deplasare Numim vector deplasare cˆampul vectorial u : Ω × R+ → R3 definit prin relat¸ia: x = χ(X, t) = X + u(X, t),

X ∈ Ω, t > 0.

(2.8)

A¸sadar, u(X, t) este deplasarea la momentul t a particulei reperatˇa prin X ˆın configurat¸ia de referint¸ˇa. ˆIn plus, ˙ v = u,

¨. a=u

Numim gradient al deplasˇarii tensorul H definit prin relat¸ia ( ∂ui ) H = ∇X u Hij = . ∂Xj Avem astfel, F = I3 + H

(Fij = δij + Hij ).

ˆIn plus, 1 G = (H + H T + H T H), 2

( Gij =

1 ( ∂ui ∂uj ∂uk ∂uk )) + + . 2 ∂Xj ∂Xi ∂Xi ∂Xj

−−→ Considerˇam un element material ˆın punctul M ∈ Ω, notat dM0 , −−→ dM0 = (dX1 , dX2 , dX3 ), −−→ de lungime dl0 ¸si de versor n0 la momentul t = 0, care devine dM la momentul t > 0, de lungime dl. Variat¸ia relativˇa a pˇatratului lungimii acestui element material se exprimˇa prin formula (dl)2 − (dl0 )2 = 2Gij n0i n0j . (dl0 )2 ˆIn plus,

( dl )2

= Cij n0i n0j . dl0 Dacˇa n0 este direct¸ia principalˇa de deformat¸ie asociatˇa valorii proprii λ a tensorului C, atunci Cij n0j = λn0i , de unde rezultˇa

( dl )2 dl0

= λ.

Valoarea proprie 1 µ = (λ − 1) 2 corespunzˇatoare tensorului G este legatˇa de variat¸ia relativˇa a patratului lungimilor (dl)2 − (dl0 )2 = 2µ. (dl0 )2 Deoarece ( dl )2 = 2µ + 1, dl0 deducem cˇa dl − dl0 √ = 1 + 2µ − 1. dl0 Astfel, pentru µ 0 aceste douˇa elemente devin dM ¸si δM, de versori n respectiv ν ¸si fac ˆıntre ele unghiul θ. −−→ dM = dl n = (dx1 , dx2 , dx3 ), −→ δM = δl ν = (δx1 , δx2 , δx3 ). Pe de o parte, −−→ −→ −−→ −−→ dM · δM − dM0 · δM0 = cosθ dl δ l − cosθ0 dl0 δl0 . Pe de altˇa parte, −−→ −→ −−→ −−→ dM · δM − dM0 · δM0 = dx · δx − dX · δX = (F T F − I 3 )dX · δX = 2Gαβ dl0 δl0 n0α ν0β .

Obt¸inem prin urmare, cos θ = √

cos θ0 + 2Gαβ n0α ν0β √ . 1 + 2Gαβ n0α n0β 1 + 2Gαβ ν0α ν0β

Presupunem cˇa n0 este direct¸ie principalˇa a deformat¸iei asociatˇa deformat¸iei normale principale µ, adicˇa: Gαβ n0β = µ n0α , ceea ce implicˇa

√ 1 + 2µ cos θ0 . cos θ = √ 1 + 2Gαβ ν0α ν0β

Notˇam cˇa dacˇa n0 ¸si ν 0 sunt vectori perpendiculari, atunci n ¸si ν sunt de asemenea vectori perpendiculari.

Tensorul deformat¸iilor infinitezimale ˆIn condit¸ii uzuale de utilizare, solidele nu sufer˘a decˆat mici deformat¸ii. Mai exact, vectorul deplasare u = u(X, t) variaz˘a lent ˆın raport cu X, derivatele ∂ui part¸iale fiind, ˆın consecint¸ˇa mici. ∂Xj Introducem urm˘atoarea ipotez˘a: Existˇa δ > 0 astfel ˆıncˆat: |H| ≤ δ, 2

δ ≪ 1,

(2.9)

2

¸si cantit˘a¸tile de ordin δ , O(δ ), sunt neglijabile, unde |H|, a se vedea (1.4), noteazˇa norma tensorului H. Pentru a simplifica scrierea, vom omite cantitˇa¸tile de ordin ˆIn aceast˘a ipotez˘a, ¸tinˆand cont c˘a H T H = O(δ 2 ) deducem

1 G = (H + H T ). 2

Aceasta justific˘a introducerea tensorului 1 ε = (H + H T ) 2

1 ∂ui ∂uj (εij = ( + )). 2 ∂Xj ∂Xi

(2.10)

Tensorul ε se nume¸ste tensorul deformat¸iilor infinitezimale. Pentru a marca dependent¸a lui ε de vectorul deplasare, uneori se va utiliza scrierea ε(u). Este important de semnalat c˘a ˆın ipoteza micilor deformat¸ii (2.9), dac˘a gradientul ∇x ϕ al unei funct¸ii scalare sau vectoriale ϕ este de ordin δ atunci ¸si ∇X ϕ este de ordin δ (¸si reciproc). ˆIn plus, ˆın acest caz ∇X ϕ = ∇x ϕ.

(2.11)

ˆIntr-adev˘ar, ∇X ϕ = (∇x ϕ)F = (∇x ϕ)(I 3 + H) = ∇x ϕ + (∇x ϕ)H. Deci, presupunˆand cˇa ∇x ϕ = O(δ), atunci (∇x ϕ)H = O(δ 2 ), de unde rezult˘a ∇X ϕ = ∇x ϕ. Reciproc, sˇa presupunem cˇa ∇X ϕ = O(δ). Cum ∇x ϕ = (∇X ϕ)(I 3 + H)−1 , admit¸ˆand (2.9), avem (I 3 + H)−1 = I 3 − H ¸si de aici rezultˇa c˘a ∇x ϕ = ∇X ϕ. ˆIn concluzie, admit¸ˆand (2.9), nu mai trebuie s˘a facem distinct¸ie ˆıntre coordonatele materiale ¸si coordonatele spat¸iale ˆın calculul derivatelor part¸iale ale lui ϕ. ˆIn particular, putem scrie: 1 ε = (∇x u + ∇Tx u) 2

1 ∂ui ∂uj (εij = ( + )). 2 ∂Xj ∂Xi

(2.12)

Semnificat¸ia fizicˇ a a tensorului ε Discutˇam ˆın cele ce urmeazˇa semnificat¸ia fizicˇa a componentelor tensorului ε. Definim alungirea relativ˘a Λ(n0 ) prin formula Λ(n0 ) =

dl − dl0 . dl0

Deducem c˘a Λ(n0 )(Λ(n0 ) + 2) = 2Gαβ n0α n0β iar de aici rezult˘a Λ(n0 ) = −1 +

√

1 + 2Gαβ n0α n0β .

ˆIn fat¸a radicalului am luat semnul plus c˘aci dacˇa punem semnul minus deducem cˇa Λ(n0 ) = −2 atunci cˆand G = 03 . Dar dacˇa G = 03 , atunci dl = dl0 ¸si prin urmare Λ(n0 ) = 0, deci contradict¸ie. Dac˘a presupunem acum c˘a n0 = (1, 0, 0), deducem √ Λ(e1 ) = −1 + 1 + 2G11 . Analog, Λ(e2 ) = −1 + ¸si Λ(e3 ) = −1 +

√

1 + 2G22

√ 1 + 2G33 .

Admit¸aˆnd (2.9), Λ(e1 ) = ε11 ,

Λ(e2 ) = ε22 ,

Λ(e3 ) = ε33 .

A¸sadar ε11 reprezint˘a alungirea relativ˘a ˆın direct¸ia axei OX1 , ε22 reprezint˘a alungirea relativ˘a ˆın direct¸ia axei OX2 , iar ε33 reprezint˘a alungirea relativ˘a ˆın direct¸ia axei OX3 , admit¸ˆand ipoteza micilor deformat¸ii. Fie acum n0 = (1, 0, 0) ¸si ν 0 = (0, 1, 0). Not˘am ˆın acest caz θ cu θ12 ; 2G12 √ cos θ12 = √ . 1 + 2G11 1 + 2G22

Folosind (2.9), rezult˘a cos θ12 = 2ε12 . Analog, cos θ13 = 2ε13 , cos θ23 = 2ε23 .

Mi¸scˇ ari izocore. Medii incompresibile ˆIn cele ce urmeazˇa analizˇam problema modificˇarii volumelor. Fie V ¸si Vt volumele ocupate de corp ˆın cofigurat¸ia nedeformat˘a, respectiv ˆın configurat¸ia deformat˘a. Avem ∫ ∫ ∫ V = dV, Vt = dv = J(X, t)dV. Ω

Ωt

Ω

Dac˘a J(X, t) = 1 X ∈ Ω, t > 0,

(2.13)

obt¸inem ∀t > 0,

Vt = V

adic˘a volumul corpului se conserv˘a ˆın timpul deformat¸iei. Reciproc, sˇa presupunem c˘a ˆın timpul deformat¸iei, la fiecare moment t > 0, volumul fiec˘arei pˇart¸i P a corpului se conserv˘a. Vom demonstra c˘a (2.13) se verific˘a. ˆIntr-adev˘ar: ∀t > 0

vt = v

(unde vt ¸si v reprezint˘a volumele ocupate de P la momentele t > 0 respectiv t = 0). Fie ωt ¸si ω domeniile din R3 ocupate de P la momentele t > 0 ¸si t = 0. Avem: ∫ ∫ vt = dv; v = dV ωt

ω

¸si deoarece ωt = χ(ω, t), utilizˆand aceea¸si schimbare de variabil˘a, egalitatea vt = v devine ∫ ∫ J(X, t)dV = ω

dV ω

∀t > 0.

(2.14)

Cum P este o parte arbitrar˘a a corpului (deci ω este o parte arbitrar˘a a lui Ω), (2.14) implic˘a (2.13). Mi¸sc˘arile χ pentru care se verificˇa (2.13) se numesc mi¸sc˘ ari izocore. Un mediu continuu care nu poate realiza decˆat mi¸scari izocore se nume¸ste mediu incompresibil. S¸tim cˇa 1+H H12 H13 11 J = det F = det(I 3 + H) = H21 1 + H22 H23 . H32 1 + H33 H31 Utilizˆand (2.9) se obt¸ine J = 1 + Hkk . Deducem astfel cˇa J = 1 + Div u.

(2.15)

ˆIn consecint¸˘a, o miscare χ este izocor˘ a dac˘a ¸si numai dac˘a Div u = 0; (am notat aici prin Div u divergent¸a cˆampului vectorial u ˆın raport cu coordonatele Lagrange). S˘a studiem variat¸ia relativ˘a a volumului. Pentru aceasta vom utiliza formula dv = JdV. Prin urmare,

dv − dV = J − 1. dV

Admit¸aˆnd (2.9), J = 1 + Div u, ¸si astfel,

dv − dV = Div u. dV

A¸sadar, ˆın ipoteza (2.9), Div u masoar˘a variat¸ia relativ˘a a elementului de volum.

Determinarea cˆ ampului deplasare Formula (2.10) ne arat˘a c˘a, definit fiind vectorul deplasare u pe Ω × R+ , componentele tensorului ε se obt¸in prin simple deriv˘ari. ˆIn aplicat¸ii apare ¸si problema invers˘a ¸si anume s˘a se determine vectorul depalsare u(X, t) cˆand se cunoa¸ste expresia tensorului ε(X, t), pe Ω × R+ . Problema revine la integrarea sistemului de ecuat¸ii cu derivate part¸iale de ordinul ˆıntˆai (2.10), εij fiind funct¸ii cunoscute pe Ω la momentul considerat t. S˘a introducem tensorul W de componente 1 ( ∂ui ∂uj ) − . 2 ∂Xj ∂Xi

Wij =

(2.16)

Presupunem c˘a ˆıntr-un punct P0 (X 0 ) se cunosc m˘arimile ui (X 0 , t) ≡ u0i ¸si Wij (X 0 , t) ≡ Wij0 . Admit¸ˆand c˘a ε este dat pe Ω la momentul dat t, ne propunem s˘a determinˇam u ˆıntr-un punct oarecare P1 (X 1 ). ∫ 1 0 ui (X , t) = ui (X , t) + dui = P0 P1

∫ =

u0i

+

(εij + Wij )dXj . P0 P1

Deoarece Wij este necunoscut pe Ω, urmeaz˘a s˘a-l elimin˘am ∫ ∫ Wij dXj = Wij d(Xj − Xj1 ) = P0 P1

P0 P1

∫ d(Wij (Xj − Xj1 )) − (Xj − Xj1 )Wij,k dXk

= P0 P1

∫ =

−Wij0 (Xj0

−

Xj1 )

−

(Xj − Xj1 )Wij,k dXk . P0 P1

Pe de altˇa parte, 2Wij,k = ui,jk − uj,ik = ui,jk + uk,ij − uk,ij − uj,ik = 2(εik,j − εjk,i ). A¸sadar, cumulˆand rezultatele obt¸inute, putem scrie ∫ 1 0 0 1 0 ui (X , t) = ui + Wij (Xj − Xj ) + Uik dXk (i ∈ {1, 2, 3}) P0 P1

unde Uik = εik + (Xj1 − Xj )(εik,j − εjk,i ). Vectorial, vom scrie ∫ u(X , t) = u + w × (X − X ) + 1

0

0

1

UdX

0

P0 P1

unde

1 Rot u. 2 Pe aceast˘a form˘a vectorial˘a se vede u¸sor c˘a partea integrat˘a reprezint˘a o rototranslat¸ie (mi¸scare de corp rigid); numai integrala propriu-zis˘a reprezint˘a o deformare. De altfel, dac˘a se calculeaz˘a εij cu w0 (W32 , W13 , W21 ) =

u(X, t) = u0 + w0 × (X 1 − X 0 ) se obt¸ine zero. Deducem a¸sadar c˘a solut¸ia sistemului εij = 0 (i, j ∈ {1, 2, 3})

este u(X 1 , t) = u0 + w0 × (X 1 − X 0 ). Not˘am de asemenea c˘a dou˘a st˘ari de deplasare c˘arora le corespund aceia¸si tensori ε difer˘a ˆıntre ele printr-o rototranslat¸ie.

Condit¸ii de compatibilitate Cum ˆın punctul P1 (X 1 ) trebuie s˘a obt¸inem deplas˘ari unice (nu putem avea dou˘a deplas˘ari) oricare ar fi drumul de integrare ˆıntre P0 ¸si P1 , pentru ca integralele curbilinii ∫ Uik dXk , i ∈ {1, 2, 3} P0 P1

s˘a nu depind˘a de drumul de integrare, este necesar ¸si suficient ca Uik dXk

(i ∈ {1, 2, 3})

s˘a fie diferent¸iale total exacte, adicˇa sˇa existe Ui astfel ˆıncˆat Uik dXk = dUi

i ∈ {1, 2, 3}.

(2.17)

Admitem c˘a Ω este simplu conex. A¸sadar, (2.17) are loc dac˘a ¸si numai dac˘a Uik,l = Uil,k (i ∈ {1, 2, 3}), (2.18) k, l ∈ {1, 2, 3}. Utilizˆand (2.2), relat¸iile (2.18) devin εik,l − δlj (εik,j − εjk,i ) + (Xj1 − Xj )(εik,jl − εjk,il ) = εil,k − δkj (εil,j − εjl,i ) + (Xj1 − Xj )(εil,jk − εjl,ik ), iar de aici deducem εik,jl + εjl,ik − εil,jk − εjk,il = 0 i, j, k, l ∈ {1, 2, 3}.

(2.19)

Avem 34 = 81 relat¸ii. Dat˘a fiind ˆıns˘a simetria lui ε ¸si posibilitatea schimb˘arii ordinii de derivare (admitem cˇa ui ∈ C 2 ) nu toate aceste relat¸ii sunt independente. Un simbol care se va dovedi util este simbolul lui Ricci, un simbol cu trei indici, ϵijk , definit dupˇa cum urmeazˇa:     1 dacˇa (i, j, k) = (1, 2, 3)     ϵijk = (2.20) −1 dacˇa (i, j, k) ̸= (1, 2, 3)        0 dacˇa cel put¸in doi indici sunt egali. prin (i, j, k) = (1, 2, 3) ˆınt¸elegˆand cˇa permutarea (i, j, k) are aceea¸si paritate cu permutarea (1, 2, 3) iar prin (i, j, k) ̸= (1, 2, 3) ˆınt¸elegˆand cˇa permutarea (i, j, k) nu are aceea¸si paritate cu permutarea (1, 2, 3). De exemplu, ϵ312 = 1, ϵ213 = −1, ϵ112 = 0. Relat¸iile (2.19) se pot rescrie cu ajutorul simbolului lui Ricci (2.20) astfel, ϵklm (εik,jl − εkj,il ) = 0 m ∈ {1, 2, 3} sau, echivalent ϵijn ϵklm εik,jl = 0 m, n ∈ {1, 2, 3}.

(2.21)

Observ˘am c˘a este vorba doar de 6 relat¸ii distincte, deoarece ˆın (2.21) avem o m˘arime care depinde doar de doi indici m ¸si n, ¸si care ˆın plus este simetric˘a. Relat¸iile (2.19) se numesc condit¸iile de compatibilitate ale lui Saint Venant (1860). Demonstrat¸ia dat˘a aici prin construct¸ia solut¸iei sistemului (2.10) se g˘ase¸ste ˆın lucr˘arile lui Cesaro (1906) ¸si Voltera (1907). Relat¸iile distincte sunt: 2ε12,12 = ε11,22 + ε22,11 (2.22) ε11,23 = (ε12,3 − ε23,1 + ε31,2 ),1 ¸si cele care se obt¸in din acestea prin permut˘ari circulare.

Existent¸a condit¸iilor de compatibilitate demonstreaz˘a faptul c˘a nu ne putem da arbitrar un tensor simetric ε pe Ω ¸si s˘a pretindem cˇa el reprezint˘a un tensor de deformat¸ie. Din punctul de vedere al integr˘arii sistemului (2.10) era natural s˘a ajungem la condit¸ii de compatibilitate deoarece ˆın (2.10) avem un sistem de 6 ecuat¸ii cu numai 3 necunoscute.

2.3

Principii generale

Principiul conservˇ arii masei Fie P o parte a unui corp M, care la momentul t = 0 ocup˘a ˆın R3 domeniul ω iar la momentul t > 0 ocup˘a domeniul ωt ⊂ R3 . ˆIn mecanica mediilor continue se presupune c˘a masa este funct¸ie absolut continu˘a de volumul suport. Pentru a exprima masa pˇart¸ii P vom utiliza Teorema Radon-Nikodym din teoria mˇasurii ¸si a integralei. Astfel, la momentul t = 0 masa p˘art¸ii P este ∫ m= ρ0 (X)dX, (2.23) ω

iar la momentul t > 0 masa p˘art¸ii P este ∫ mt = ρ(x, t)dx.

(2.24)

ωt

Aplicat¸ia ρ0 : Ω → R + se nume¸ste densitate de mas˘a ˆın configurat¸ia nedeformat˘ a. Aplicat¸ia ρ(·, t) : Ωt → R+

(2.25)

(2.26)

se nume¸ste densitate de mas˘a ˆın configurat¸ia deformat˘a. Presupunˆand c˘a masa p˘art¸ii P se conserv˘a ˆın timp, adicˇa m = mt

∀t > 0,

(2.27)

avem

∫

∫ ρ(x, t)dx ∀t > 0.

ρ0 (X)dX = ω

(2.28)

ωt

Deoarece ωt = χ(ω, t),

∫

∫ ρdx =

ρJdX.

ωt

(2.29)

ω

Din (2.28) si (2.29) rezultˇa ∫

∫ ρJdX =

ω

ρ0 (X)dX. ω

Cum ω este un domeniu arbitrar ˆın Ω (P fiind o parte oarecare a corpului), deducem ρ0 = ρJ (2.30) sau, explicit, ρ0 (X) = ρ(χ(X, t), t)J(X, t),

∀X ∈ Ω, t > 0.

Numim (2.30) principiul conserv˘ arii masei ˆın coordonate materiale. Derivˆand ˆın raport cu t ˆın (2.30), obt¸inem ˙ 0 = ρJ ˙ + ρJ. Sˇa demonstrˇam cˇa ˙ = div v dv. dv

(2.31)

Pentru aceasta, utiliz˘am formula dv = JdV care, prin derivare ˆın raport cu timpul se scrie ˙ = JdV. ˙ dv

(2.32)

Urmˇatoarea teoremˇa se va dovedi utilˇa ˆın cele ce urmeazˇa. Teorema 2.9. (Euler) J˙ = J div v.

(2.33)

Demonstrat¸ie : J = det F =

∂x1 ∂X1

∂x1 ∂X2

∂x1 ∂X3

∂x2 ∂X1

∂x2 ∂X2

∂x2 ∂X3

∂x3 ∂X1

∂x3 ∂X2

∂x3 ∂X3

Prin derivare obt¸inem x˙ 1,1 x˙ 1,2 x˙ 1,3 d ˙ J = J = x2,1 x2,2 x2,3 dt x3,1 x3,2 x3,3

x 1,1 x1,2 x1,3 = x 2,1 x2,2 x2,3 x3,1 x3,2 x3,3

x 1,1 x1,2 x1,3 + x˙ 2,1 x˙ 2,2 x˙ 2,3 x3,1 x3,2 x3,3

.

x 1,1 x1,2 x1,3 + x 2,1 x2,2 x2,3 x˙ 3,1 x˙ 3,2 x˙ 3,3



= J1 + J2 + J3 unde, am notat cu Ji determinantul obt¸inut din J prin ˆınlocuirea liniei i cu derivata ei ˆın raport cu timpul. Observˇam cˇa ∂x ∂x ∂x ∂v1 ∂v1 ∂v1 ∂v1 · j ∂v1 · j ∂v1 · j ∂X1 ∂X2 ∂X3 ∂xj ∂X1 ∂xj ∂X2 ∂xj ∂X3 = ∂x2 ∂x2 ∂x2 J1 = ∂x2 ∂x2 ∂x2 = ∂X1 ∂X2 ∂X3 ∂X1 ∂X2 ∂X3 ∂x3 ∂x3 ∂x3 ∂x3 ∂x3 ∂x3 ∂X1 ∂X2 ∂X3 ∂X1 ∂X2 ∂X3 =

∂v1 ∂x1

·

∂x1 ∂X1

∂v1 ∂x1

·

∂x1 ∂X2

∂v1 ∂x1

·

∂x1 ∂X3

∂x2 ∂X1

∂x2 ∂X2

∂x2 ∂X3

∂x3 ∂X1

∂x3 ∂X2

∂x3 ∂X3

Analog, deducem J2 =

∂v2 J, ∂x2

∂v1 = ∂x1 J.

J3 =

∂v3 J. ∂x3

A¸sadar, ∂v1 ∂v2 ∂v3 J˙ = J+ J+ J = J div v. ∂x1 ∂x2 ∂x3  Din (2.32) ¸si (2.33) rezult˘a (2.31). Deducem ρ˙ + ρ div v = 0,

(2.34)

sau explicit, ρ(x, ˙ t) + ρ(x, t)(

∂v1 ∂v2 ∂v3 + + ) = 0 ∀x ∈ Ωt , t > 0. ∂x1 ∂x2 ∂x3

Numim (2.34) principiul conserv˘ arii masei ˆın coordonate spat¸iale.

Principiul lui Cauchy ˆIn mecanic˘a act¸iunile se traduc prin fort¸e. Suntem interesat¸i aici de fort¸ele care act¸ioneaz˘a asupra unui subsistem material P. ˆIntr-o configurat¸ie deformatˇa, la un moment oarecare t > 0, avem • act¸iuni de contact (sau fort¸e de suprafat¸a˘) care sunt act¸iunile ce se exercit˘a pe ∂ωt din partea sistemului M\P (M fiind corpul ˆın intregime); • act¸iuni la distant¸˘a. Ca s˘a scriem ecuat¸ia de mi¸scare a lui P va trebui ca act¸iunile lui M\P asupra lui P, act¸iuni de contact pe suprafat¸a ∂ωt , s˘a le traducem prin fort¸e de contact pe ∂ωt . F˘acˆand aceast˘a operat¸ie, ˆın ecuat¸iile de mi¸scare ale lui P vom introduce aceste fort¸e ¸si vom face abstract¸ie de prezent¸a lui M\P. Principiul lui Cauchy: Exist˘a o distribut¸ie de fort¸e τ pe suprafat¸a ∂ωt a c˘aror act¸iune asupra lui P este echivalent˘a cu act¸iunea lui M\P. ˆIn plus,

vectorul τ este acela¸si pentru toate suprafet¸ele de aceea¸si orientare ¸si de acela¸si plan tangent ˆın x. Pe suprafat¸a ∂ωt , pentru un moment dat t > 0, τ va depinde de punctul curent. ˆIntr-un punct dat x ¸si pentru un moment dat t, τ va depinde de suprafat¸a ce trece prin acel punct. Dependent¸a de suprafat¸a˘ se traduce prin dependent¸a de versorul normal ν. ˆIn concluzie, putem scrie τ = τ (x, ν, t). Vectorul τ se nume¸ste vectorul tensiune Cauchy. Act¸iunea lui M\P asupra lui P este un cˆamp absolut continuu de suprafat¸ˇa. Astfel, aceast˘a act¸iune (total˘a) se va scrie ∫ τ da. ∂ωt

Act¸iunile la distant¸a˘ sunt funct¸iile absolut continue de mas˘a, ˆın timp ce fort¸ele electromagnetice vor fi funct¸ii absolut continue de volum. Act¸iunea total˘a care se exercit˘a asupra lui P ˆın configurat¸ia ωt va fi ∫ ρ b dv. ωt

unde ρ b reprezint˘a fort¸a pe unitatea de volum. Not˘am c˘a ˆın general b = b(x, t). Dac˘a b nu depinde explicit de t, cˆampul se va numi stat¸ionar, iar dac˘a b nu depinde de x, se va numi omogen. Un exemplu de cˆamp omogen este b = g (cˆampul gravitat¸ional). Pentru a simplifica scrierea vom nota f = ρ b. Rezumˆand, un corp M, care la momentul dat t > 0 se gˇase¸ste ˆın configurat¸ia deformat˘a Ωt , este supus la urm˘atoarele tipuri de fort¸e: • fort¸e volumice: de densitate f (·, t) : Ωt → R3 ; a: de densitate h(·, t) : ∂Ωt → R3 ; • fort¸e externe de suprafat¸˘ • fort¸e de suprafat¸˘a: de densitate τ (·, ·, t) : Ωt × S → R3 ,

unde S = {ν ∈ R3 ||ν| = 1}.

(2.35)

Dac˘a x ∈ ∂Ωt iar ν coincide cu vectorul normalei exterioare la Ωt ˆın x, atunci τ (x, ν, t) = h(x, t).

Principiul variat¸iei impulsului Impulsul subsistemului material P se define¸ste prin formula ∫ ∫ H(P) = v dm = ρ v dv. P

(2.36)

ωt

Principiul variat¸iei impulsului: pentru orice subsistem P, ˆın orice configurat¸ie a sa ωt , derivata impulsului ˆın raport cu timpul este egal˘a cu rezultanta fort¸elor care act¸ioneaz˘a asupra subsistemului, adic˘a ∫ ∫ ∫ d ρ vdv = f dv + τ da. (2.37) dt ωt ωt ∂ωt Not˘am c˘a (2.37) implic˘a ∫ ∫ ∫ ρ adv = f dv + ωt

ωt

τ da.

(2.38)

∂ωt

ˆIntr-adev˘ar, d dt

∫

∫ d ρ v dv = ρ vJdV dt ωt ∫ ω = (ρvJ ˙ + ρaJ + ρv div vJ)dV ∫ω = [(ρ˙ + ρ div v)vJ + ρaJ]dV. ω

Utilizˆand principiul conserv˘arii masei ˆın coordonate spat¸iale deducem c˘a ∫ ∫ ∫ d ρvdv = ρa JdV = ρa dv. dt ωt ω ωt

Principiul variat¸iei momentului cinetic Momentul cinetic al subsistemului P ˆın raport cu un pol O se define¸ste prin formula ∫ ∫ KO (P) = x × vdm = x × vρdv. P

ωt

Principiul variat¸iei momentului cinetic: pentru orice subsistem P, ˆın orice configurat¸ie a sa ωt , derivata momentului cinetic, calculat ˆın raport cu un pol O, este egalˇa cu momentul rezultant al fort¸elor care act¸ioneaza asupra subsistemului, calculat ˆın raport cu acela¸si pol O, adicˇa ∫ ∫ ∫ d x × ρ v dv = x × τ da + x × f dv. (2.39) dt ωt ∂ωt ωt Relat¸ia (2.39) implic˘a ∫ ∫ x × ρ a dv = ωt

∫ x × τ da +

x × f dv.

∂ωt

(2.40)

ωt

ˆIntr-adev˘ar: ∫ d x × ρ v dv dt ωt ∫ d x × ρ v JdV = dt ω ∫ ˙ + x × ρv div vJ)dV = (x˙ × ρvJ + x × ρvJ ˙ + x × ρvJ ∫ω = [v × ρvJ + x × (ρ˙ + ρ div v)vJ + x × ρaJ]dV ∫ω = x × ρaJdV ω ∫ = x × ρ a dv. ω

Pe componente, (2.40) se scrie cu ajutorul simbolului lui Ricci (2.20) astfel, ∫ ∫ ∫ ϵijk xj ρ ak dv = ϵijk xj τk da + ϵijk xj fk dv, (2.41) ωt

pentru fiecare i ∈ {1, 2, 3}.

∂ωt

ωt

Capitolul 3 Modelare prin ecuat¸ii cu derivate part¸iale

3.1

Ecuat¸iile Cauchy. Condit¸ii la limitˇ a

Obiectivul acestei sect¸iuni este scrierea ecuat¸iilor care modeleazˇa mi¸scarea unui solid deformabil. Prezentˇam ˆın cele ce urmeazˇa o teoremˇa fundamentalˇa. Teorema 3.10. (Teorema lui Cauchy) Pentru fiecare t > 0, admitem cˇa ρ(·, t) : Ωt → R+ , a(·, t) : Ωt → R3 ¸si f (·, t) : Ωt → R3 sunt funct¸ii continue iar τ (·, ·, t) : Ωt × S → R3 , τ = τ (x, ν, t) este funct¸ie de clas˘a C 1 ˆın raport cu x pentru oricare ν ∈ S ¸si continu˘ a ˆın raport cu ν pentru oricare x ∈ Ωt . Atunci exist˘a T (·, t) : Ωt → M3 de clas˘a C 1 (M3 fiind mult¸imea matricilor p˘atratice de ordin 3, cu componente reale) astfel ˆıncˆ at τ (x, ν, t) = T (x, t)ν 39

x ∈ Ωt , ν ∈ S, t > 0

(3.1)

S3 S S1

S2

Figura 3.1: Tetraedrul Tt ρ a = div T + f

x ∈ Ωt , t > 0

T (x, t) = T T (x, t) x ∈ Ωt , t > 0,

(3.2) (3.3)

unde S este mult¸imea definitˇa ˆın (2.35). Demonstrat¸ie : Fie t > 0, x ∈ Ωt ¸si ν ∈ S (de componente νi > 0). −→ Considerˇam tetraedrul Tt din Figura 3.1, SSi fiind vectori de versori ei , i ∈ {1, 2, 3}. Introducem urm˘atoarele notat¸ii: (SS2 S3 ) = F1 ; aria F1 = a1 , (SS3 S1 ) = F2 ; aria F2 = a2 , (SS1 S2 ) = F3 ; aria F3 = a3 , (S1 S2 S3 ) = F ; aria F = a. Not˘am c˘a ai = νi a ∀i ∈ {1, 2, 3}. ˆIntr-adev˘ar, cu argumente elementare de geometrie obt¸inem rapid, de exemplu c˘a a3 = ν3 a. Pentru aceasta scriem ′

a3 = aria(SS1 S2 ) =

|SS |·|S1 S2 | , 2

a = aria(S3 S1 S2 ) =

|S3 S |·|S1 S2 | , 2

′

unde S ′ este proiect¸ia lui S pe S1 S2 , ¸si de aici, a3 = e3 · ν a = ν3 a, e3 fiind −−→ versorul vectorului SS3 ¸si ˆın acela¸si timp versorul axei 0X3 . Analog pentru i ∈ {1, 2}. Utilizˆand (2.38) pentru ωt = Tt cu ∂Tt = F¯1 ∪ F¯2 ∪ F¯3 ∪ F¯ , vom avea, pe componente ∫ ∫ (ρ ai (x, t) − fi (x, t))dv = Tt

τi (x, θ, t)da

(3.4)

∂Tt

unde θ = −ej pentru Fj , j ∈ {1, 2, 3} ¸si θ = ν pentru F. Existˇa y i ∈ F astfel ˆıncˆat ∫ 1 τi (x, ν, t)da = τi (y i , ν, t) a F ¸si existˇa y ij ∈ Fj astfel ˆıncˆat ∫ 1 τi (x, −ej , t)da = τi (y ij , −ej , t). aj Fj Prin urmare, ∫

3 ∑ τi (x, θ, t)da = (aria F )τi (y i , ν, t) + (aria Fj )τi (y ij , −ej , t).

∂Tt

j=1

Utilizˆand (3.4) rezult˘a ∫ (ρai (x, t) − fi (x, t))dv = (τi (y i , ν, t) + Tt

3 ∑

νj τi (y ij , −ej , t))(aria F )

j=1

¸si de aici deducem |τi (y i , ν, t) +

∑3 j=1

νj τi (y ij , −ej , t)|(aria F ) ≤

≤ supz ∈Tt |ρ(z, t)ai (z, t) − fi (z, t)| vol(Tt ).

(3.5)

Reamintim aici c˘a ν este arbitrar fixat (cu νi > 0). Facem s˘a tind˘a vˆarfurile Sk la S (caz ˆın care y i ¸si y ij tind la x), vol(Tt ) = aria F

d(S,(S1 S2 S3 ))·aria F 3

aria F

=

d(S, (S1 S2 S3 )) → 0, cˆand Sk → S, k ∈ {1, 2, 3}. 3

Trecˆand la limit˘a Sk → S, k ∈ {1, 2, 3}, din (3.5) rezult˘a τi (x, ν, t) +

3 ∑

τi (x, −ej , t)νj = 0 ∀i ∈ {1, 2, 3}

j=1

sau, vectorial τ (x, ν, t) +

3 ∑

τ (x, −ej , t)νj = 0.

(3.6)

j=1

Din (3.6) rezult˘a imediat τ (x, ν, t) = −τ (x, −e1 , t)ν1 − τ (x, −e2 , t)ν2 − τ (x, −e3 , t)ν3 .

(3.7)

ˆIn (3.7) trecem la limit˘a ν → e1 ¸si obt¸inem τ (x, e1 , t) = −τ (x, −e1 , t). Analog, obt¸inem τ (x, e2 , t) = −τ (x, −e2 , t) ¸si τ (x, e3 , t) = −τ (x, −e3 , t). A¸sadar, τ (x, ej , t) = −τ (x, −ej , t) ∀j ∈ {1, 2, 3}.

(3.8)

Din (3.6) ¸si (3.8) rezult˘a τ (x, ν, t) = τ (x, ej , t)νj

(3.9)

τi (x, ν, t) = τi (x, ej , t)νj .

(3.10)

sau pe componente

Definim Tij (·, t) : Ωt → R astfel ˆıncˆat Tij (x, t) = τi (x, ej , t) ∀x ∈ Ωt , t > 0, ∀i, j ∈ {1, 2, 3}.

(3.11)

Utilizˆand (3.10) ¸si (3.11) rezult˘a τi (x, ν, t) = Tij (x, t)νj

∀i ∈ {1, 2, 3}.

(3.12)

Din (3.12) rezult˘a (3.1). Pentru a demonstra (3.2) vom utiliza (3.1) ¸si (2.38). ˆIntr-adev˘ar, din (2.38) deducem imediat ∫ ∫ (ρa − f )dv = τ da. ωt

∂ωt

Utilizˆand acum (3.1) ¸si formula Gauss-Ostrogradski obt¸inem ∫ ∫ ∫ (ρ a − f )dv = T (x, t)νda = div T dv. ωt

∂ωt

ωt

Cum ωt subdomeniu arbitrar ˆın Ωt deducem (3.2) care reprezintˇa ecuat¸iile de mi¸scare ˆın configurat¸ia deformat˘a. Pentru a demonstra simetria tensorului T vom utiliza (2.41). ˆIntradev˘ar, din (2.41) deducem imediat ∫ ∫ ϵijk xj (ρ ak − fk )dv = ϵijk xj τk da. ωt

∂ωt

ˆIn plus, utilizˆand (3.2), obt¸inem ∫ ∫ ϵijk xj Tkm,m dv = ωt

ϵijk xj Tkn νn da.

∂ωt

Utilizand formula Gauss-Ostrogradski putem scrie ∫ ∫ ϵijk xj Tkm,m dv = ϵijk (xj,n Tkn + xj Tkn,n )dv ωt

iar de aici

∫

ωt

∫ ϵijk xj Tkm,m dv =

ωt

∫ ϵijk δjn Tkn dv +

ωt

ϵijk xj Tkn,n dv. ωt

Deducem a¸sadar cˇa

∫ ϵijk Tkj dv = 0. ωt

Cum ωt subdomeniu arbitrar in Ωt , deducem imediat ϵijk Tkj = 0 ∀i ∈ {1, 2, 3}.

(3.13)

Pentru i = 1, ϵ123 T32 + ϵ132 T23 = 0 ⇒ T32 = T23 . Pentru i = 2, ϵ231 T13 + ϵ213 T31 = 0 ⇒ T13 = T31 . Pentru i = 3, ϵ312 T21 + ϵ321 T12 = 0 ⇒ T21 = T12 . A¸sadar, T este tensor simetric.



Teorema lui Cauchy stabile¸ste ˆın primul rˆand faptul c˘a vectorul tensiune al lui Cauchy depinde liniar de ν. Apoi, pune ˆın evident¸a˘ pentru orice x ∈ Ωt ¸si orice t > 0 tensorul simetric T (x, t) care se nume¸ste tensorul tensiune al lui Cauchy. Teorema lui Cauchy furnizeaz˘a de asemenea ecuat¸iile de mi¸scare ˆın configurat¸ie deformat˘a. Aceste ecuat¸ii nu sunt ˆıns˘a utilizabile pentru a calcula miscarea χ c˘aci ele se scriu ˆın variabile Euler x = (x1 , x2 , x3 ) care reprezint˘a de fapt necunoscuta problemei (x = χ(X, t)). Pentru aceasta, se impune scrierea ecuat¸iilor de mi¸scare ˆın configurat¸ia de referint¸˘a. Acest lucru se va realiza ˆın continuare, prin introducerea a doi tensori Piola - Kirchhoff.

Tensorii Piola-Kirchhoff Pornind de la tensorul tensiune al lui Chauchy T ¸si utilizˆand gradientul deformat¸iei F , construim urm˘atorul tensor numit primul tensor Piola -

Kirchhoff S = JT F −T (3.14) S(X, t) = J(X, t)T (χ(X, t), t)F

−T

(X, t),

unde J = det F . Deoarece F −T noteazˇa (F −1 )T , putem scrie (3.14) pe componente dupˇa cum urmeazˇa, Sij = JTik (F −1 )jk ∀i, j ∈ {1, 2, 3}. (3.15) Notˇam aici faptul c˘a Div S = J div T ,

(3.16)

sau, pe componente, ∂Sij ∂Tij =J ∂Xj ∂xj

∀i ∈ {1, 2, 3}.

(3.17)

Introducem notat¸ia urm˘atoare f 0 = Jf . f 0 (X, t) = J(X, t)f (χ(X, t), t). ˆInmult¸ind cu J ecuat¸iile de mi¸scare Cauchy ˆın configurat¸ia deformat˘a ¸si utilizˆand (3.16), precum ¸si principiul conserv˘arii masei ˆın configurat¸ie nedeformat˘a, (2.30), deducem Div S + f 0 = ρ0 a ˆın Ω × (0, ∞)

(3.18)

unde f 0 = Jf . Primul tensor Piola - Kirchhoff nu este simetric. Se impune introducerea unui nou tensor, al doilea tensor Piola - Kirchoff, care se va dovedi a fi simetric. Pornind de la primul tensor Piola - Kirchhoff S ¸si utilizˆand tensorul gradient al deformatiei F , construim acum tensorul numit al doilea tensor

Piola - Kirchhoff: Π = F −1 S; (3.19) Π(X, t) = F

−1

(X, t)S(X, t).

Pe componente, (3.19) se scrie astfel: Πij = (F −1 )ik Sjk ,

i, j ∈ {1, 2, 3}.

(3.20)

Al doilea tensor Piola- Kirchhoff este simetric. ˆIntr-adev˘ar, Π = F −1 S = F −1 JT F −T = JF −1 T F −T ¸si ΠT = (JF −1 T F −T )T = J(F −T )T T T (F −1 )T . Cum T = T T iar (F −T )T = F −1 , deducem Π = ΠT . Din (3.19) rezult˘a imediat S = F Π.

(3.21)

Combinˆand (3.18) ¸si (3.21) rezult˘a Div (F Π) + f 0 = ρ0 a ˆın Ω × (0, ∞),

(3.22)

care reprezint˘a ecuat¸iile de mi¸scare ˆın configurat¸ia nedeformat˘ a. Suntem interesat¸i acum s˘a scriem ecuat¸iile de mi¸scare ˆın ipoteza micilor deformat¸ii (2.9), caz ˆın care   H = O(δ),  

F = I 3 + O(δ), F

−T

= I 3 + O(δ),

F −1 = I 3 − H = I 3 + O(δ), J = 1 + O(δ).

(3.23)

Admitem ˆın plus ˆın cele ce urmeazˇa cˇa T = O(δ).

(3.24)

Utilizˆand acum (3.23) ¸si (3.14), obt¸inem S = JT F −T = (1 + O(δ))T (I 3 + O(δ)). A¸sadar, S = T (I 3 + O(δ) + O(δ 2 )).

(3.25)

Din (3.24) ¸si (3.25) deducem S = T.

(3.26)

De asemenea (3.19) ¸si (3.23) implic˘a Π = F −1 S = (I 3 + O(δ))S.

(3.27)

Cum din (3.26) ¸si (3.24) deducem S = O(δ), ¸si utilizˆand (3.27), rezult˘a imediat c˘a Π = S.

(3.28)

ˆIn concluzie, ˆın teoria liniarizat˘a nu se mai face distinctie ˆıntre tensorii T , S, Π (T = S = Π). Introducem urm˘atoarea notat¸ie: σ = S = Π = T; (3.29) σ(X, t) = S(X, t) = Π(X, t) = T (χ(X, t), t) ∀x ∈ Ω, ∀t > 0. Notˇam cˇa σ = σ T . Tensorul σ se va numi tensorul tensiune Cauchy ˆın teoria liniarizatˇa.

Ecuat¸iile Cauchy Fie T > 0. Suntem interesat¸i sˇa descriem evolut¸ia unui solid deformabil ˆın intervalul [0, T ]. Pentru aceasta, considerˇam ecuat¸ia vectorialˇa ¨ ˆın Ω × (0, T ), Div σ + f 0 = ρ0 u

(3.30)

ˆın care u : Ω × [0, T ] → R3 ¸si σ : Ω × [0, T ] → S3 reprezint˘a necunoscutele problemei, iar ρ0 : Ω → R+ ¸si

f 0 : Ω × [0, T ] → R3

reprezint˘a datele problemei. Ecuat¸ia (3.30) a fost scris˘a pe domeniul Ω care este fix ¸si cunoscut. Dac˘a dorim s˘a cunoa¸stem geometria structurii deformate vom utiliza relat¸ia de definit¸ie a vectorului deplasare x = X + u(X, t) ∀t > 0. Pentru a u¸sura scrierea vom nota ˆın continuare prin x = (x1 , x2 , x3 ) coordonatele unei particule ˆın configurat¸ia de referint¸a˘ ¸si vom omite indicele 0 la ρ0 ¸si f 0 . Astfel, ˆın cele ce urmeaz˘a (3.30) se va scrie astfel Div σ + f = ρ¨ u ˆın Ω × (0, T ).

(3.31)

Procesele de evolut¸ie modelate prin (3.31) se numesc procese dinamice. Daca v = 0, procesele modelate prin Div σ + f = 0

(3.32)

se numesc procese statice. Dac˘a u variaz˘a lent ˆın timp, astfel ˆıncˆat termenul ρ¨ u poate fi neglijat, procesele modelate prin (3.32) se numesc procese cvasistatice. Ecuat¸ia (3.31) se nume¸ste ecuat¸ia de mi¸scare Cauchy iar (3.32) se nume¸ste ecuat¸ie de echilibru Cauchy.

Γ2 x

3

Γ1

O

Ω

x2 x1

Figura 3.2: Domeniul Ω ¸si o partit¸ie a frontierei sale

Condit¸ii la limitˇ a Fie Γ frontiera lui Ω ¸si ν = ν(x) versorul normalei exterioare la Γ ˆın punctul ¯1 ∪ Γ ¯ 2 , Γ1 ∩ Γ2 = ∅, ca ˆın Figura x. Fie de asemenea o partit¸ie a lui Γ, Γ = Γ 3.2. Consider˘am urm˘atoarele condit¸ii: u = g

pe Γ1 × (0, T );

(3.33)

σν = h

pe Γ2 × (0, T ).

(3.34)

Condit¸ia (3.33) se nume¸ste condit¸ie la limit˘a ˆın deplas˘ari: cˆampul deplasare u are valoare impus˘a pe port¸iunea Γ1 a frontierei Γ, funct¸ia g : Γ1 × (0, T ) → R3 fiind dat˘a a problemei. Condit¸ia (3.34) se nume¸ste condit¸ie la limit˘a ˆın tract¸iuni: vectorul tensiune Cauchy σν este prescris pe port¸iunea Γ2 a frontierei Γ, funct¸ia h : Γ2 × (0, T ) → R3 fiind dat˘a a problemei. Dac˘a Γ1 = ∅ problema la limit˘a se va numi problem˘ a ˆın tract¸iuni (pur˘a) iar dac˘a Γ2 = ∅ problema la limit˘a se va numi problem˘ a ˆın deplas˘ari (pur˘a). Dac˘a Γ1 ̸= ∅ si Γ2 ̸= ∅, atunci problema considerat˘a este o problem˘a mixt˘a. Pe frontierˇa pot fi considerate ¸si alte tipuri de condit¸ii la limitˇa: condit¸ii de contact unilateral sau bilateral cu o fundat¸ie rigidˇa sau deformabilˇa, contactul putˆandu-se modela cu sau fˇarˇa frecare. Lucrarea de fat¸a se limiteazˇa la condit¸iile pe frontierˇa clasice, dar pentru cititorul interesat ¸si de alte tipuri de condit¸ii la limitˇa facem trimitere la monografiile [15, 38, 39].

F

S

l0

Figura 3.3: Barˇa supusˇa la ˆıncercˇari experimentale

3.2

Legi de comportament

Ecuat¸ia (3.31) este insuficient˘a pentru a descrie evolut¸ia unui solid deformabil. Pe de o parte, din punct de vedere matematic, este clar c˘a din 3 ecuat¸ii scalare nu putem determina 9 necunoscute: u1 , u2 , u3 , σ11 , σ12 , σ13 , σ22 , σ23 , σ33 . Pe de alt˘a parte, din punct de vedere fizic, remarc˘am faptul c˘a ecuat¸ia (3.31) este consecint¸a˘ a unor principii general valabile pentru mediile continue ¸si, prin urmare, dac˘a ecuat¸ia (3.31) ar fi suficient˘a pentru a descrie mi¸scarea unui anumit mediu continuu, ar ˆınsemna ca toate corpurile au acela¸si comportament indiferent de tipul de material din care sunt confect¸ionate. Trebuie a¸sadar, ca ecuat¸iei (3.31) s˘a i se al˘ature o ecuat¸ie care s˘a descrie ceea ce este propriu fiec˘arui tip de material din care este confect¸ionat corpul: legea de comportament. O lege de comportament ( sau lege constitutivˇa) este ˆın general o relat¸ie ˆıntre tensorul tensiune σ, tensorul deformat¸iilor infinitezimale ε ¸si derivatele ˙ lor temporale σ˙ ¸si ε. Legile de comportament se scriu ˆın general pe baza unor observat¸ii experimentale. O serie ˆıntreag˘a de ˆıncerc˘ari se realizeaz˘a pentru a se scrie o lege de comportament: ˆıncerc˘ari de ˆınc˘arcare-monoton˘a, ˆıncerc˘ari de relaxare, fluaj, ˆıncerc˘ari de ˆınc˘arcare-desc˘arcare, ˆıncerc˘ari pe care le vom prezenta pe scurt, ˆın ceea ce urmeaz˘a. Sˇa consider˘am o bar˘a de sect¸iune cu aria S, lungime l0 , c˘areia i se aplic˘a o fort¸˘a F = F (t) la una din extremit˘a¸ti, ˆın timp ce cealalt˘a extremitate r˘amˆane fix˘a; a se vedea Figura 3.3.

σ

σ

ε

ε σ

σ

ε

ε

Figura 3.4: ˆIncarcare monotona Se m˘asoar˘a alungirea bazei l = l(t) ¸si se definesc: ε(t) =

l(t) − l0 l0

¸si σ(t) =

F (t) . S

(3.35)

ˆ arcarea monoton˘a Se cre¸ste progresiv fort¸a ¸si se m˘asoar˘a l. Se • Inc˘ calculeaz˘a σ ¸si ε ¸si se deseneaz˘a curba σ = σ(ε). Pot apare situat¸iile din Figura 3.4. • Fluajul Se ˆıncepe cu o ˆıncercare de ˆınc˘arcare monoton˘a astfel ˆıncˆat la un moment t (pe care ˆın continuare ˆıl vom considera moment init¸ial) avem σ(0) = σ0 , ε(0) = ε0 , F (0) = F0 . De acum F se ment¸ine constant˘a (F = F0 ). Observˇam cˇa σ va r˘amˆane constant pentru t > 0. ˆIn general, ε cre¸ste ˆın timp: este fluajul. Pot apare situat¸iile din Figura 3.5. • Relaxarea Ca ¸si la fluaj se ˆıncepe cu o ˆıncercare de ˆınc˘arcare monoton˘a astfel ˆıncˆat la un moment t (pe care ˆıl vom presupune ˆın continuare moment init¸ial) avem: σ(0) = σ0 ,

ε(0) = ε0 .

σ

σ

σ

ε

ε

ε

Figura 3.5: Fluaj σ

σ

ε

σ

ε

ε

Figura 3.6: Relaxare Ment¸inem de aceast˘a dat˘a deformat¸ia constant˘a ˆın timp. ˆIn general se constat˘a c˘a tensiunile scad ˆın timp: este relaxarea; a se vedea Figura 3.6. ˆ arcare-desc˘arcare Se act¸ioneaz˘a cu o fort¸˘a F cresc˘atoare pˆan˘a la un • Inc˘ moment dat, apoi fort¸a F descre¸ste la zero. Aceast˘a ˆıncercare permite punerea ˆın evident¸˘a a comportamentului elastic sau plastic. Dac˘a curbele de ˆınc˘arcare-desc˘arcare σ = σ(ε) coincid, mediul se nume¸ste elastic; ˆın caz contrar el este plastic ¸si dup˘a desc˘arcare complet˘a se evident¸iaz˘a o deformat¸ie rezidual˘a. Prin urmare, materialele plastice sunt caracterizate printr-o descompunere a deformat¸iei ε ˆıntr-o parte elasticˇa (reversibilˇ a) εe ¸si o parte plasticˇa (ireversibilˇ a) εp . Altfel spus, putem scrie ε = εe + εp . Este posibil ca anelasticitatea materialului sˇa se manifeste de la un anumit prag numit limitˇ a de proport¸ionalitate; dacˇa acest prag este fix, mediul se nume¸ste perfect plastic; dacˇa acest prag variazˇa, se produce

σ

σ

ε

ε

Figura 3.7: ˆInc˘arcare-desc˘arcare ecruisaj; a se vedea Figura 3.7 Legile de comportament se scriu astfel ˆıncˆat, ˆın dimensiunea 1 s˘a fie descrise fenomenele puse ˆın evident¸˘a prin ˆıncerc˘ari experimentale cum ar fi cele patru ˆıncerc˘ari prezentate anterior. ˆIn prezenta lucrare sunt puse ˆın discut¸ie numai legile de comportament care corespund materialelor liniar elastice, caz ˆın care tensorul tensiune este o funct¸ie liniarˇa de tensorul deformat¸iilor infinitezimale ε, σ = F (ε).

(3.36)

ˆIn cazul unu-dimensional, legea (3.36) poate modela anumite proprietˇa¸ti puse ˆın evident¸aˇ prin ˆıncercˇari experimentale de ˆıncˇarcare-descˇarcare cum ar fi de exemplu, liniaritatea curbei σ = σ(ε). Relaxarea ¸si fluajul ˆınsˇa nu pot fi descrise de legea (3.36). ˆIntr-adevˇar, sˇa ne plasˇam ˆın cazul unu-dimensional; dacˇa de exemplu la momentul t = 0 avem ε(0) = ε0 ¸si ment¸inem deformat¸ia constantˇa ε(t) = ε0

∀t > 0

rezultˇa cˇa σ(t) = F (ε0 ) ∀t > 0. Prin urmare, modelul (3.36) nu poate descrie fenomenul de relaxare pus ˆın evident¸a prin ˆıncercˇari experimentale. Un comentariu similar se poate face pentru fenomenul de fluaj.

De asemenea, pentru ecuat¸ia (3.36), ˆın cazul unu-dimensional, curbele ˆıncˇarcare-descˇarcare σ = σ(ε) coincid. Deci, acest model nu poate descrie deformat¸iile reziduale. Cititorul interesat de alte legi de comportament decˆat cele elastice (3.36), poate consulta de exemplu [15, 38, 39]. ˆIn general, funct¸ia F din (3.36) depinde de punctul x ∈ Ω (Ω fiind domeniul din Rd ocupat de corp), σ(x, t) = F (x, ε(x, t)) ∀x ∈ Ω, t > 0. Dacˇa F nu depinde explicit de x, mediul se nume¸ste omogen. ˆIn caz contrar, el se nume¸ste neomogen. ˆIn elasticitatea liniarˇa σ este o funct¸ie liniarˇa de ε, adicˇa σ = Eε

(σij = Eijkh εkh )

(3.37)

unde E = (Eijkh ) este un tensor de ordin patru, adicˇa o aplicat¸ie liniarˇa de la spat¸iul vectorial al tensorilor cartezieni de ordinul doi ˆın el ˆınsu¸si. Componentele sale Eijkh , se numesc coeficient¸i elastici ¸si sunt independente de tensorul deformat¸ie. ˆIn cazul neomogen, Eijkh depind de x ∈ Ω ¸si ˆın cazul omogen, Eijkh sunt constante.

Legea lui Hooke ˆIn teoria clasicˇa a elasticitˇa¸tii nu sunt vizate materialele liniar elastice ˆın toatˇa generalitatea lor, considerˆandu-se cazul omogen ¸si izotrop. Un mediu continuu se nume¸ste izotrop dacˇa are acelea¸si propriet˘a¸ti ˆın orice direct¸ie ˆın jurul unui punct arbitrar fixat. Matematic, proprietatea de izotropie se traduce astfel: dac˘a deformat¸iei ε ˆıi corespunde tensiunea σ, σij = Eijkh εkh , atunci deformat¸iei ε′ = QεQT ,

unde Q este o matrice de schimbare de repere ortogonale QQT = QT Q = I 3 , ˆıi va corespunde tensiunea σ ′ = QσQT . Legat de proprietatea de izotropie a unui material vom demonstra urmˇatoarea teoremˇa. Teorema 3.11. Dac˘a un material este liniar elastic ¸si izotrop atunci funct¸ia scalar˘a w(ε) definit˘a prin w(ε) = σij εij = Eijkh εkh εij satisface w(QεQT ) = w(ε), Demonstrat¸ie :

∀Q matrice ortogonal˘ a.

Cum materialul este izotrop, putem scrie:

σ = Eε;

σ ′ = QσQT ;

ε′ = QεQT ;

σ ′ = Eε′ .

Rezult˘a a¸sadar w(QεQT ) = w(ε′ ) = tr(σ ′ ε′ ) = tr(QσQT QεQT ) = tr(QσεQT ) = Qik σim εmj Qkj = Qik Qkj σim εmj = δij σim εmj = σim εmi = σij εij = tr(σε) = w(ε).  Funct¸ia scalarˇa w = w(ε) definitˇa ˆın (3.38) se nume¸ste funct¸ie izotrop˘ a. Teorema 3.12. Dac˘a w(ε) este funct¸ie izotrop˘ a ¸si form˘a p˘atratic˘ a ˆın raport cu ε, atunci ea se scrie astfel w(ε) = A(tr(ε))2 + Btr(ε2 ),

(3.38)

A ¸si B fiind scalari independent¸i de ε. Demonstrat¸ie : w(ε) = w(QεQT ) ∀Q matrice ortogonal˘a Cum ε este un tensor simetric, el posed˘a trei valori proprii reale εI , εII , εIII ¸si trei direct¸ii principale νI , νII , νIII , care se pot alege ortogonale dou˘a cˆate dou˘a. ˆIn acest reper principal, ε se scrie   0  εI 0   εp =   0 εII 0   0 0 εIII

   .   

Fie Q matricea de trecere de la baza {e1 , e2 , e3 } la baza {ν I , ν II , ν III }; w(ε) = w(QεQT ) = w(εp ). Rezult˘a c˘a w(ε) depinde doar de εI , εII , εIII ¸si aceast˘a dependent¸˘a este simetric˘a ˆın raport cu cele 3 argumente. T ¸ inˆand cont ¸si de faptul c˘a w(ε) este form˘a p˘atratic˘a ˆın raport cu ε, rezult˘a c˘a w arat˘a astfel w(ε) = A′ (ε2I + ε2II + ε2III ) + 2B ′ (εI εII + εII εIII + εI εIII ) unde A′ si B ′ sunt constante. Avem de asemenea     ε2I + ε2II + ε2III = tr(ε2 );    2(εI εII + εII εIII + ε1 εIII ) = (trε)2 − tr(ε2 ), de unde rezult˘a imediat (3.38) schimbˆand notat¸iile constantelor. Presupunem cˇa Eijkh = Ekhij = Ejikh .



Exist˘a ˆın elasticitatea liniarizatˇa o energie intern˘a (de deformat¸ie) e :=

1 Eijkh εkh εij 2ρ0

astfel cˇa σij = ρ0

∂e . ∂εij

Cum w(ε) = e este funct¸ie izotrop˘a ¸si form˘a p˘atratic˘a ˆın raport cu ε, utilizˆand teorema anterioar˘a, deducem 1 Eijkh εkh εij = A(trε)2 + Btr(ε2 ) 2ρ0 sau echivalent σij εij = 2ρ0 [A(trε)δij εij + Bεij εij ]. Rezult˘a de aici urm˘atoarea expresie pentru componenta ij a tensorului tensiune σ : σij = ρ0 [2A(trε)δij + 2Bεij ], sau schimbˆand notat¸iile σ = λ(tr ε)I 3 + 2µε (σij = λ(trε)δij + 2µεij ).

(3.39)

Legea (3.39) este cunoscut˘a sub numele de legea lui Hooke. Coeficient¸ii λ ¸si µ se numesc coeficient¸ii lui Lam´ e. Dac˘a materialul este presupus omogen, ace¸sti coeficient¸i sunt independent¸i de x. Se constat˘a c˘a ˆın cazul izotrop avem numai doi coeficient¸i elastici independent¸i, ˆın loc de 21 de coeficient¸i Eijkh ˆın cazul anizotrop oarecare. Ace¸sti coeficient¸i Eijkh se exprim˘a ˆın funct¸ie de λ ¸si µ prin relat¸ia urm˘atoare Eijkh = λδij δkh + µ(δik δjh + δih δjk ),

i, j, k, h ∈ {1, 2, 3}.

a1111 = λ + 2µ

a2211 = λ

a3311 = λ

a1212 = µ

a1112 = 0

a2212 = 0

a3312 = 0

a1213 = 0

a1113 = 0

a2213 = 0

a3313 = 0

a1223 = 0

a1122 = 0

a2222 = λ + 2µ

a3322 = 0

a1323 = 0

a1123 = 0

a2223 = 0

a3323 = 0

a2323 = 0

a1133 = 0

a2233 = λ

a3333 = λ + 2µ a1313 = µ.

Legea lui Hooke inversat˘ a Formula (3.39) indic˘a dependent¸a tensorului tensiune de tensorul deformat¸iilor infinitezimale ε. Ne propunem acum sa exprim˘am tensorul deformat¸iilor infinitezimale ˆın funct¸ie de tensorul tensiune σ. Utilizˆand urma unui tensor, deducem din (3.39) trσ = 3λtrε + 2µtrε. De aici rezult˘a imediat trε =

1 trσ. 3µ + 2λ

(3.40)

Utilizˆand (3.39) ¸si (3.40) deducem acum formula urm˘atoare εij =

1 λ σij − δij , 2µ 2µ(3λ + 2µ)

ε=

1 λ σ− I 3. 2µ 2µ(3λ + 2µ)

sau tensorial, (3.41)

Relat¸ia deformat¸ie-tensiune (3.41) poartˇa numele de Legea lui Hooke inversatˇa.

Ecuat¸iile Navier-Lam´ e Dupˇa cum v˘azut, pentru materiale liniar elastice, omogene ¸si izotrope, tensorul tensiune al lui Cauchy se scrie pe componente σij = λ(trε)δij + 2µεij

∀i, j ∈ {1, 2, 3}.

Utilizˆand relat¸iile geometrice, 1 εij = (ui,j + uj,i ) 2

∀i, j ∈ {1, 2, 3},

deducem σij = λuk,k δij + µ(ui,j + uj,i ). Calculul divergent¸ei tensorului tensiune σ conduce la forma urm˘atoare a componentei i a vectorului Div σ : σij,j = λuk,kj δij + µ(ui,jj + uj,ij ) = λuk,ki + µuk,ik + µui,jj = (λ + µ)(uk,k ),i + µ△ui . A¸sadar, ecuat¸iile de echilibru se scriu (λ + µ)

∂ (div u) + µ△ui + fi = 0 ∂xi

∀i ∈ {1, 2, 3}.

Ecuat¸iile scalare anterioare, se scriu restrˆans sub forma vectorial˘a astfel (λ + µ) grad(div u) + µ△u + f = 0.

(3.42)

Ecuat¸iile de mi¸scare ale lui Cauchy se scriu ˆın acest caz (λ + µ)

∂ (div u) + µ△ui + fi = ρ¨ ui ∂xi

∀i ∈ {1, 2, 3},

iar vectorial, (λ + µ) grad(div u) + µ△u + f = ρ¨ u.

(3.43)

Ecuat¸ia vectorialˇa (3.42) se nume¸ste ecuat¸ia Navier iar ecuat¸ia vectorialˇa (3.43) se nume¸ste ecuat¸ia Lam´e.

Ecuat¸ia lui Navier este util˘a atunci cˆand se ¸stie c˘a vectorul deplasare c˘autat este cu rotat¸ional nul, adic˘a rot rot u = 0.

(3.44)

rot rot u = grad(div u) − △u.

(3.45)

Reamintim aici formula

Din (3.42) ¸si (3.45) deducem urm˘atoarea form˘a pentru ecuat¸ia lui Navier, (λ + 2µ) grad(div u) − µ rot rot u + f = 0.

(3.46)

Dac˘a vectorul deplasare u c˘autat este cu rotat¸ional nul, (3.46) se reduce la (λ + 2µ) grad(div u) + f = 0, care implic˘a rot f = 0 ¸si, prin urmare, f deriv˘a dintr-un potent¸ial, adic˘a f = − grad ϕ. Ecuat¸ia (3.47) devine ˆın acest caz grad((λ + 2µ) div u) − grad ϕ = 0, ecuat¸ie ce admite urm˘atoarea integral˘a prim˘a (λ + 2µ) div u − ϕ = const., ceea ce ˆın general simplific˘a c˘autarea lui u.

(3.47)

Ecuat¸iile Beltrami-Michell Contrar act¸iunii anterioare, de eliminare a tensiunilor ˆın vederea determin˘arii deplas˘arilor, ne propunem acum s˘a cautam direct tensorul tensiune. Tensorul tensiune trebuie s˘a satisfac˘a ecuat¸iile de mi¸scare (sau de echilibru) ¸si condit¸iile la limit˘a ˆın tensiuni. Unui astfel de tensor tensiune, (σij ), ˆıi putem asocia un tensor e de componente eij =

1+ν ν σij − (trσ)δij . E E

(3.48)

Acest tensor de ordinul 2, e = (eij ), este un tensor al deformat¸iilor liniarizate numai dac˘a ˆıi putem asocia un vector deplasare u = (ui ) astfel ˆıncˆat 1 ( ∂ui ∂uj ) eij = + = εij (u). 2 ∂xj ∂xi

(3.49)

Potrivit condit¸iilor de compatibilitate Saint-Venant, acest lucru este posibil dac˘a ¸si numai dac˘a ϵipq ϵjrs epr,qs = 0

∀i, j ∈ {1, 2, 3}.

(3.50)

Ne propunem acum s˘a exprim˘am aceste condit¸ii de compatibilitate ˆın funct¸ie de tensorul tesiune. Este suficient s˘a ˆınlocuim ˆın (3.50) eij prin expresia (3.48) obt¸inˆand, ϵipq ϵjrs epr,qs =

1+ν ν ϵipq ϵjrs σpr,qs − ϵipq ϵjrs σmn,qs δpr . E E

(3.51)

Deducem astfel ecuat¸iile Beltrami-Michell (ecuat¸ii de compatibilitate ˆın tensiuni), (1 + ν)ϵipq ϵjrs σpr,qs − νϵipq ϵjrs σmm,qs δqr = 0. Observˇam cˇa (1 + ν)(σij,kl + σkl,ij − σil,jk − σjk,il ) = ν[σmm,kl δij + σmm,ij δkl − σmm,jk δil − σmm,il δjk ].

Dacˇa k = l, atunci, (1 + ν)(σij,kk + σkk,ij − σik,jk − σjk,ik ) = ν[σmm,kk δij + σmm,ij δkk − σmm,jk δik − σmm,ik δjk ] iar de aici, deducem (1 + ν)[σij,kk + σkk,ij − σik,jk − σjk,ik ] = ν[σmm,kk δij + σmm,ij ]. Utilizˆand ecuat¸iile de echilibru, obt¸inem σij,kk +

1 ν σkk,ij + fi,j + fj,i = σmm,kk δij 1+ν 1+ν

sau, echivalent △σij +

3.3

1 ν σkk,ij − △σmm δij = −(fi,j + fj,i ). 1+ν 1+ν

Modele simplificate

Modelul deformat¸iilor plane Dac˘a ˆıntr-un corp elastic, omogen ¸si izotrop, de form˘a cilindric˘a, a se vedea Figura 3.8, cu generatoarele paralele cu axa Ox3 , vectorul deplasare are componentele de urmˇatoarea formˇa, u1 = u1 (x1 , x2 ),

u2 = u2 (x1 , x2 ),

u3 = 0,

atunci, tensorul deformat¸iilor infinitezimale are urmˇatoarele componente: ∂u1 ∂u2 , ε22 = , ε33 = 0, ∂x1 ∂x2 1 ∂u1 ∂u2 + ), ε23 = ε31 = 0. = ( 2 ∂x2 ∂x1

ε11 = ε12

(3.52)

x1 0 x2

x3

Figura 3.8: Modelul plan Tensorul ε de componente (3.52) se nume¸ste tensorul deformat¸iilor plane. Pentru materiale elastice, omogene ¸si izotrope, tensorul tensiune Cauchy are urmˇatoarea formˇa,    σ11 σ21 0   σ=  σ21 σ22 0   0 0 σ33

   ,   

unde σαβ = λ(trε)δαβ + 2µεαβ , α, β ∈ {1, 2}. Cum 0 = ε33 =

1+ν ν σ33 − (trσ)δ33 , E E

deducem (1 + ν)σ33 − ν(σ11 + σ22 + σ33 ) = 0, iar de aici σ33 = ν(σ11 + σ22 ).

(3.53)

Fie Ω o sect¸iune dreapt˘a oarecare a lui B. Problema ˆın acest caz revine la a rezolva o problem˘a de elasticitate bidimensionalˇa     σαβ,β + fα = 0 ˆın Ω ⊂ R2 α, β ∈ {1, 2}     (Pd ) : σαβ = λεγγ δαβ + 2µεαβ ˆın Ω α, β ∈ {1, 2}        condit¸ii la limit˘a pe ∂Ω.

Necunoscutele acestei probleme sunt: σ11 , σ12 , σ21 , σ22 , ε11 , ε12 , ε21 , ε22 , u1 , u2 .

(3.54)

Dup˘a rezolvarea problemei Pd , σ33 se determin˘a, folosind (3.53).

Modelul tensiunilor plane Prin definit¸ie, un tensor cartezian de ordin doi se nume¸ste cˆ amp de tensiuni plane dac˘a are forma urm˘atoare    σ11 (x1 , x2 ) σ12 (x1 , x2 ) 0      . σ= σ (x , x ) σ (x , x ) 0  12 1 2  22 1 2     0 0 0 Admit¸ˆand c˘a lucr˘am cu un corp liniar elastic ¸si izotrop, utilizˆand legea lui Hooke, σαβ = λ(ε11 + ε22 + ε33 )δαβ + 2µεαβ ,

∀α, β ∈ {1, 2}

(3.55)

0 = ε23 = ε13 ,

(3.56)

0 = λ(ε11 + ε22 + ε33 ) + 2µε33 .

(3.57)

Din (3.57) rezult˘a ε33 = −

λ (ε11 + ε22 ). λ + 2µ

(3.58)

Utilizˆand (3.58) putem rescrie (3.55) astfel σαβ = λ∗ εγγ δαβ + 2µεαβ , unde λ∗ =

2λµ . λ + 2µ

(3.59)

ˆIn acest caz problema la limitˇa se scrie astfel.     σαβ,β + fα = 0 ˆın Ω     (Pt ) : σαβ = λ∗ εγγ δαβ + 2µεαβ ˆın Ω        condit¸ii la limit˘a pe ∂Ω iar ε33 se determin˘a cu (3.58). Remarca 3.13. Din punct de vedere matematic, problemele de tensiuni plane ¸si cele de deformat¸ii plane sunt de aceea¸si natur˘a. De aceea, ˆın Capitolul 4, vom ˆınt¸elege prin problemˇ a planˇa urmˇatoarea problemˇa la limtˇa     Div σ + f = 0 ˆın Ω     (P ) : σ = Cε(u) ˆın Ω        condit¸ii la limit˘a pe ∂Ω, fˇarˇa a mai face distinct¸ie ˆıntre cazul deformat¸iilor plane ¸si cel al tensiunilor plane. Odatˇa determinatˇa solut¸ia problemei (P), (u1 , u2 , σ11 , σ12 , σ21 , σ22 ), ˆın cazul deformat¸iilor plane σ33 se determinˇa cu (3.53) iar ˆın cazul tensiunilor plane ε33 se determinˇa cu (3.58.)

Modelul antiplan Fie un corp care ocupˇa un domeniu cilindric B ⊂ R3 , raportat la un sistem cartezian ortogonal Ox1 x2 x3 ; a se vedea Figura 3.9. Cosiderˇam cˇa generatoarele sunt paralele cu axa Ox3 ¸si vom nota prin Ω domeniul mˇarginit al lui R2 care reprezintˇa sect¸iunea transversalˇa a cilindrului. Presupunem cˇa cilindrul este suficient de lung pentru a putea neglija efectul fort¸elor aplicate la capete; astfel putem considera B = Ω × (−∞, ∞). Fie Γ = ∂Ω. Presupunem cˇa frontiera Γ este netedˇa ¸si considerˇam o partit¸ie ˆın trei zone mˇasurabile

x1

Ω

0 x2

x3

Figura 3.9: Modelul antiplan Γ1 , Γ2 ¸si Γ3 , astfel ˆıncˆat mes Γ1 > 0. Notˇam cu ν versorul normalei exterioare la Γ. Studiem starea de deformat¸ie a corpului material datoratˇa fort¸elor volumice de densitate f 0 ¸si fort¸elor de tract¸iune de densitate f 2 , ¸stiind cˇa pe Γ1 × (−∞, ∞) cilindrul este ˆıncastrat ˆıntr-o structurˇa fixˇa. Admitem ipoteza micilor deformat¸ii ¸si presupunem cˇa densitatea fort¸elor volumice ¸si densitatea tract¸iunilor au urmˇatoarea formˇa particularˇa. f 0 = (0, 0, f0 )

cu f0 = f0 (x1 , x2 ) : Ω → R,

(3.60)

f 2 = (0, 0, f2 )

cu f2 = f2 (x1 , x2 ) : Γ2 → R.

(3.61)

Admit¸ˆand aceastˇa ˆıncˇarcare particularˇa, vectorul deplasare este de forma u = (0, 0, u)

cu u = u(x1 , x2 , t) : Ω × [0, T ] → R,

(3.62)

astfel cˇa, tensorul deformat¸iilor infinitezimale are urmˇatoarea expresie   1 u 0  0 2 ,1        1 ε(u) =  0 (3.63) 0 u . 2 ,2       1 1 u u 0 2 ,1 2 ,2 Folosind legea constitutivˇa a lui Hooke σ = λ tr(ε)I + 2µε, ecuat¸iile de mi¸scare ale lui Cauchy se reduc la urmˇatoarea ecuat¸ie scalarˇa, div(µ(x)∇u(x)) + f0 (x) = ρ0 (x)¨ u(x)

ˆın Ω.

(3.64)

Dacˇa materialul este omogen, atunci µ este o constantˇa strict pozitivˇa ¸si ecuat¸ia (3.64) se rescrie ˆın cazul stat¸ionar astfel, µ△u(x)) + f0 (x) = ρ0 (x)¨ u(x)

ˆın Ω.

Deoarece cilindrul este ˆıncastrat pe Γ1 × (−∞, +∞), vectorul deplasare se anuleazˇa aici. De asemenea, (3.62) implicˇa u=0

pe Γ1 .

(3.65)

Cˆat prive¸ste versorul normalei exterioare, ˆın contextul antiplan el este de urmˇatoarea formˇa, ν = (ν1 , ν2 , 0) cu νi = νi (x1 , x2 ) : Γ → R, i = 1, 2.

(3.66)

Deducem cˇa pentru materiale liniar elastice ¸si izotrope, σν = 0 ¸si σν = (0, 0, µ unde

∂u ), ∂ν

∂u = ∇u · ν. Prin urmare, condit¸ia la limitˇa in tensiuni se reduce la ∂ν µ(x)

∂u (x) = f2 (x) ∂ν

pe Γ2 .

Capitolul 4

Probleme clasice ˆın elastostaticˇ a

4.1

Comprimarea uniform˘ a

Cadrul fizic (Solid elastic supus la presiune hidrostatic˘ a) Fie un corp liniar elastic, omogen ¸si izotrop, care ocup˘a o regiune Ω, raportat la un sistem cartezian ortogonal Ox1 x2 x3 , a se vedea Figura 4.1. Corpul este introdus ˆıntr-un gaz avˆand presiunea constant˘a p. Admitem c˘a fort¸ele volumice sunt neglijabile. Suntem interesat¸i s˘a determin˘am vectorul deplasare, tensorul tensiune ¸si variat¸ia relativ˘a a volumului. Modelare matematic˘a prin E.D.P. Deoarece presiunea gazului este constant˘a, perturbat¸iile exterioare care act¸ioneaz˘a asupra solidului sunt independente de timp. Astfel, deplas˘arile, deformat¸iile ¸si tensiunile sunt funct¸ii 69

−p n

x3

−p n

Ω O x2 x

1

−p n

−p n

Figura 4.1: Comprimarea uniformˇa independente de timp: u = u(x),

ε = ε(x),

σ = σ(x).

ˆIn acest caz (elastostatic), componentele vectorului deplasare ui , componentele tensorului deformat¸iilor infinitezimale εij ¸si componentele tensorului tensiune σij , ui = ui (x1 , x2 , x3 ) = ui (xm ) = ui (x)

∀i ∈ {1, 2, 3}

εij = εij (x1 , x2 , x3 ) = εij (xm ) = εij (x)

∀i, j ∈ {1, 2, 3}

σij = σij (x1 , x2 , x3 ) = σij (xm ) = σij (x)

∀i, j ∈ {1, 2, 3}

satisfac ecuat¸iile urm˘atoare: • ecuat¸iile de echilibru ale lui Cauchy σij,j = 0

∀i ∈ {1, 2, 3} ˆın Ω,

• ecuat¸iile geometrice 1 εij = (ui,j + uj,i ) 2 unde ui,j =

∂ui ; ∂xj

∀i, j ∈ {1, 2, 3} ˆın Ω,

• legea de comportament pentru materiale liniar elastice, omogene ¸si izotrope (legea lui Hooke) σij = λεkk δij + 2µεij

∀i, j ∈ {1, 2, 3} ˆın Ω,

unde λ > 0 si µ > 0 sunt doua constante de material; • condit¸ii la limit˘a ˆın tract¸iuni σij νj = −pνi

∀i ∈ {1, 2, 3} pe ∂Ω.

Rezumˆand, avem urm˘atoarea problem˘a mecanic˘a: Problema 4.14.

¯ → R3 ¸si σ : Ω ¯ → R3 astfel ˆıncˆat: Sˇa se determine u : Ω Div σ = 0,

ˆın Ω

ε = 12 (∇u + ∇uT )

ˆın Ω (4.1)

σ = λ(trε)I 3 + 2µε, ˆın Ω σn = −pn,

pe ∂Ω.

Fie σ tensorul de componente σij = −pδij ˆın Ω.

(4.2)

Evident, pentru σ dat ˆın (4.2), condit¸iile la limit˘a descrise de relat¸ia a patra a lui (4.1) sunt verificate. De asemenea, prima relat¸ie din (4.1) se verific˘a. Din a treia relat¸ie a lui (4.1), deducem εij = −

p δij , 3λ + 2µ

deoarece trσ = −3p. A¸sadar, tensorul ε este un tensor sferic ¸si constant, ε=−

p I 3, 3λ + 2µ

unde I 3 este tensorul unitate. Evident, relat¸iile de compatibilitate SaintVenant sunt satisfacute ˆın aceastˇa situat¸ie. Componentele vectorului deplasare sunt: ui = −

p xi 3λ + 2µ

∀i ∈ {1, 2, 3},

(4.3)

neglijˆand deplasarea rigid˘a infinitezimal˘a. ˆIntr-adev˘ar:   p  p  u = − ⇒ u = − x1 + f1 (x2 , x3 );  1 1,1 3λ+2µ   3λ + 2µ   p p u x2 + f2 (x1 , x3 ); 2,2 = − 3λ+2µ ⇒ u2 = −  3λ + 2µ     p  p  x3 + f3 (x1 , x2 ), ⇒ u3 = −  u3,3 = − 3λ+2µ 3λ + 2µ ¸si de aici,  ∂f1 (x2 , x3 ) ∂f2 (x1 , x3 )    u1,2 + u2,1 = 0 ⇒ + = 0;   ∂x2 ∂x1   ∂f1 (x2 , x3 ) ∂f3 (x1 , x2 ) u1,3 + u3,1 = 0 ⇒ + = 0;  ∂x3 ∂x1       u2,3 + u3,2 = 0 ⇒ ∂f2 (x1 , x3 ) + ∂f3 (x1 , x2 ) = 0. ∂x3 ∂x2 Mai mult, ∂ 2 f1 (x2 , x3 ) = 0 ⇒ f1 (x2 , x3 ) = a1 (x3 )x2 + b1 (x3 ); ∂x22 ∂ 2 f3 (x1 , x2 ) = 0 ⇒ f2 (x1 , x2 ) = a2 (x2 )x1 + b2 (x2 ); ∂x21 ∂ 2 f2 (x1 , x3 ) = 0 ⇒ f3 (x1 , x3 ) = a3 (x1 )x3 + b3 (x1 ); ∂x23 ∂ 2 f1 (x2 , x3 ) = 0 ⇒ f1 (x2 , x3 ) = c1 (x2 )x3 + d1 (x2 ); ∂x23

(4.4)

∂ 2 f2 (x1 , x3 ) = 0 ⇒ f2 (x1 , x3 ) = c2 (x3 )x1 + d2 (x3 ); ∂x21 ∂ 2 f3 (x1 , x2 ) = 0 ⇒ f3 (x1 , x2 ) = c3 (x1 )x2 + d3 (x1 ); ∂x22 ∂f1 (x2 , x3 ) = a′1 (x3 )x2 + b′1 (x3 ); ∂x3 ∂ 2 f1 (x2 , x3 ) 0= = a′′1 (x3 )x2 + b′′1 (x3 ). 2 ∂x3 A¸sadar, pentru oricare x2 , a′′1 (x3 )x2 + b′′1 (x3 ) = 0.

(4.5)

Din (4.5) deducem a′′1 (x3 ) = 0;

b′′1 (x3 ) = 0

astfel cˇa a1 (x3 ) = α1 x3 + β1 ,

b1 (x3 ) = γ1 x3 + δ1 .

Deducem a¸sadar, u1 (x1 , x2 , x3 ) = −

p x1 + α1 x2 x3 + β1 x2 + γ1 x3 + δ1 . 3λ + 2µ

(4.6)

Prin permut˘ari circulare ale indicilor g˘asim: u2 (x1 , x2 , x3 ) = −

p x2 + α2 x3 x1 + β2 x3 + γ2 x1 + δ2 , 3λ + 2µ

(4.7)

u3 (x1 , x2 , x3 ) = −

p x3 + α3 x1 x2 + β3 x1 + γ3 x2 + δ3 . 3λ + 2µ

(4.8)

ˆIn aceste relat¸ii, α1 , α2 , α3 , β1 , β2 , β3 , γ1 , γ2 , γ3 , δ1 , δ2 , δ3 sunt constante arbitrare.

Funct¸iile (4.6)-(4.8) trebuie s˘a satisfac˘a (4.4). Rezult˘a a¸sadar α1 x3 + β1 + α2 x3 + γ2 = 0

∀x3 ,

α2 x1 + β2 + α3 x1 + γ3 = 0

∀x1 ,

α3 x2 + β3 + α1 x2 + γ1 = 0

∀x2 ,

¸si de aici,

    α1 + α2 = 0    β1 + γ2 = 0,     α2 + α3 = 0    β2 + γ3 = 0,     α3 + α1 = 0    β3 + γ1 = 0.

Prin urmare, α1 = α2 = α3 = 0. Introducem notat¸iile γ1 = ω2 ,

γ2 = ω3 ,

γ 3 = ω1 .

Cu aceste notat¸ii, u1 =

−p x1 + b1 + ω2 x3 − ω3 x2 3λ + 2µ

−p x2 + b2 + ω3 x1 − ω1 x3 3λ + 2µ −p u3 = x3 + b3 + ω1 x2 − ω2 x1 , 3λ + 2µ u2 =

sau, echivalent u=

−p x + b + ω × x. 3λ + 2µ

Neglijˆand deplasarea rigid˘a infinitezimal˘a (care nu modific˘a nici tensorul deformat¸iilor infinitezimale ε, nici tensorul tensiune σ) obt¸inem formula anunt¸at˘a pe componente ˆın (4.3). Sˇa introducem notat¸ia 3K = 3λ + 2µ

(4.9)

unde constanta K se nume¸ste modulul de rigiditate la comprimare. Semnificat¸ia fizic˘a pentru modulul de rigiditate la comprimare K apare imediat dac˘a calcul˘am variat¸ia relativ˘a a volumului, V − V0 3p p = div u = uk,k = trε = − =− . V0 3λ + 2µ K

(4.10)

Relat¸ia (4.10) arat˘a c˘a pentru o presiune p dat˘a, variat¸ia relativ˘a a volumului este invers proport¸ional˘a cu modulul de rigiditate la comprimare K. ˆIn plus, remarc˘am c˘a modulul K are dimensiunea unei presiuni. Experimental se poate observa c˘a o presiune extern˘a (p > 0) asupra unui solid deformabil antreneaz˘a o mic¸sorare a volumului corpului, ceea ce matematic implic˘a 3λ + 2µ > 0. 1 Compresibilitatea materialului se m˘asoar˘a prin : materialul este cu atˆat K mai compresibil cu cˆat K este mai mic. Invers, un material este mai put¸in compresibil cu cˆat K este mai mare. Un material incompresibil este un material limit˘a, corespunzator cazului K → ∞. Observat¸iile experimentale au condus la urm˘atoarele inegalit˘a¸ti ce indic˘a pozitivitatea constantelor lui Lam´e λ > 0,

µ > 0.

(4.11)

De aici rezult˘a o consecint¸a˘ important˘a pentru energia de deformat¸ie: inegalit˘a¸tile (4.11) implic˘a faptul c˘a energia de deformat¸ie e este o form˘a p˘atratic˘a pozitiv definit˘a (ˆın raport cu componentele tensorului deformat¸iilor infinitezimale ε). ˆIntr-adevˇar, cum energia de deformat¸ie este e=

1 σij εij , 2ρ0

utilizˆand legea lui Hooke, σij = λ(trε)δij + 2µεij , deducem e=

1 [λ(trε)2 + 2µtr(ε2 )]. 2ρ0

De aici rezult˘a imediat

µ εij εij , ρ0 unde µ, ρ0 > 0, deci e este form˘a p˘atratic˘a pozitiv definit˘a. Pornind de la acest rezultat bazat pe observat¸ii experimentale vom formula o ipotez˘a pentru materiale liniar elastice, dar nu neap˘arat omogene sau izotrope. Considerˆand un material liniar elastic, al c˘arui comportament este modelat prin legea constituitiv˘a e≥

σij = Eijkh εkh ,

i, j ∈ {1, 2, 3},

energia de deformat¸ie are expresia urm˘atoare 1 ρ0 e = Eijkh εkh εij . 2 ˆIn teoria clasicˇa a elasticitˇa¸tii se admite ipoteza c˘a energia de deformat¸ie e este form˘a p˘atratic˘a pozitiv definit˘a, i.e., ∃α0 > 0 a.ˆı. Eijkh εkh εij ≥ α0 εij εij

4.2

∀ε ∈ S3 .

Tract¸iunea simpl˘ a

Cadrul fizic. Fie o grind˘a cilindric˘a de lungime L limitat˘a de sect¸iunile drepte Γ0 ¸si Γ1 ca ˆın Figura 4.2. Raport˘am grinda la un sistem cartezian ortogonal Ox1 x2 x3 astfel ˆıncˆat Γ0 este ˆın planul x1 = 0 iar Γ1 este ˆın planul x1 = L. Grinda este supus˘a unor fort¸e de tract¸iune uniform repartizate,

x2 Γ0

Γ1

0

−F

x1

F

x3 L

Figura 4.2: Tract¸iunea simplˇa (F, 0, 0) pe Γ1 ¸si (−F, 0, 0) pe Γ0 , paralele cu axa cilindrului. Fort¸ele volumice sunt neglijabile iar suprafat¸a lateral˘a Γl a cilindrului nu este supus˘a nici unei fort¸e. Suntem interesat¸i s˘a determinam vectorul deplasare ¸si tensorul tensiune. Modelare matematic˘a prin E.D.P. Perturbat¸iile exterioare fiind independente de timp, problema mecanic˘a se traduce matematic prin urm˘atoarea problem˘a elastostatic˘a, admit¸aˆnd c˘a grinda este liniar elastic˘a, omogen˘a ¸si izotrop˘a. ¯ ⊂ R3 ¸si σ : Ω ¯ → S3 astfel ˆıncˆat Problema 4.15. Sˇa se determine u : Ω Div σ = 0 σ = λ(trε)I 3 + 2µε 1 ε = (∇u + ∇uT ) 2 σν = 0

in Ω (4.12) in Ω (4.13) in Ω (4.14) pe Γl (4.15)

(σν)1 = F, (σν)2 = 0, (σν)3 = 0

pe Γ1 (4.16)

(σν)1 = −F, (σν)2 = 0, (σν)3 = 0

pe Γ0 .(4.17)

Deoarece ν(x) = (1, 0, 0) pe Γ1 , din (4.16) deducem σ11 = F,

σ21 = σ31 = 0 pe Γ1 .

(4.18)

De asemenea, ν(x) = (−1, 0, 0) pe Γ0 , ¸si prin urmare, (4.17) implicˇa σ11 = F,

σ21 = σ31 = 0 pe Γ0 .

(4.19)

ˆIn plus, pe Γl , ν(x) = (0, ν2 (x1 , x2 ), ν3 (x2 , x3 )). Prin urmare, din (4.15), deducem σi2 ν2 + σi3 ν3 = 0 ∀i ∈ {1, 2, 3}.

(4.20)

Ecuat¸ia (4.12) reprezint˘a ecuat¸ia de echilibru ˆın ipoteza c˘a fort¸ele volumice sunt neglijabile; relat¸ia (4.13) este legea de comportament pentru materiale liniar elastice, omogene ¸si izotrope; ecuat¸ia vectorial˘a (4.14) este ecuat¸ia geometric˘a iar relat¸iile (4.15)-(4.17) sunt condit¸iile la limit˘a. Fie tensorul σ avˆand urmˇatoarele componente, σ11 = F

σ21 = σ31 = σ22 = σ33 = σ23 = 0.

(4.21)

Evident acest tensor verific˘a (4.12) ¸si (4.15)-(4.17). Utilizˆand (4.21), (4.13) se scrie astfel: F = λtrε + 2µε11 , (4.22) 0 = λtrε + 2µε22 ,

(4.23)

0 = λtrε + 2µε33 ,

(4.24)

ε23 = ε31 = ε12 = 0.

(4.25)

Din (4.22)-(4.24) deducem F = 3λtrε + 2µtrε, adic˘a trε =

F . 3λ + 2µ

Din (4.26) ¸si (4.22) obt¸inem F( λ ) F (λ + µ) ε11 = 1− = . 2µ 3λ + 2µ µ(3λ + 2µ)

(4.26)

(4.27)

Utilizˆand (4.23) ¸si (4.24) avem ε22 = ε33 = −

λF . 2µ(3λ + 2µ)

(4.28)

A¸sadar, tensorul deformat¸iilor infinitezimale este de forma urmˇatoare      ε=   

(λ+µ) F µ(3λ+2µ)

0

0

0

−λ F 2µ(3λ+2µ)

0

0

0

−λ F 2µ(3λ+2µ)

   .   

(4.29)

Not˘am pe de o parte faptul c˘a ε este ˆın acest caz un tensor diagonal ¸si prin urmare direct¸iile principale vor fi cele ale axelor de coordonate. Pe de alt˘a parte, ε este ˆın acest caz un tensor constant, ¸si astfel condit¸iile de compatibilitate Saint-Venant Voltera vor fi evident satisfacute. Pentru a deduce componentele vectorului deplasare vom integra sistemul urm˘ator de ecuat¸ii cu derivate part¸iale   λ+µ   u1,1 = F    µ(3λ + 2µ)     λ   u2,2 = u3,3 = − F 2µ(3λ + 2µ) (4.30)     u1,2 + u2,1 = 0         u1,3 + u3,1 = 0. Neglijˆand deplas˘arile rigide infinitezimale obt¸inem u1 = ε11 x1 ,

u2 = ε22 x2 ,

u3 = ε33 x3 .

(4.31)

Observˇam cˇa planul x1 = 0 r˘amˆane invariant ca ¸si punctele axei Ox1 . ˆIn consecint¸˘a am reu¸sit s˘a verific˘am Problema 4.15 cu ajutorul • cˆampului deplasare dat ˆın (4.31); • tensorului deformat¸iilor infinitezimale dat ˆın (4.29); • tensorului tensiune Cauchy dat ˆın (4.21).

Analiza solut¸iei obt¸inute. Modulul lui Young. Coeficientul lui Poisson. Alungirea ∆L a barei sub efectul fort¸elor de tract¸iune este dat˘a prin deplasarea punctului (L, 0, 0). ∆L = ε11 L =

λ+µ F L, µ(3λ + 2µ)

(4.32)

(x1 − X1 = u1 (X, t) = ∆L = ε11 X1 = ε11 L), astfel cˇa alungirea relativ˘a este ∆L λ+µ = F = ε11 . L µ(3λ + 2µ)

(4.33)

Introducem notat¸ia E=

µ(3λ + 2µ) λ+µ

(4.34)

unde E este modulul de rigiditate la alungire, numit modulul lui Young. Pentru F dat, alungirea relativ˘a ∆L F = ε11 = L E

(4.35)

este cu atˆat mai mic˘a cu cˆat E este mai mare. Modulul lui Young este un coeficient de material care se m˘asoar˘a ˆın laborator pe ma¸sina de tract¸iune ¸si care are dimensiunea unei fort¸e pe unitatea de suprafat¸a˘ (i.e., a unei tensiuni sau a unei presiuni), N/m2 . Experimental s-a observat c˘a • E > 0, • dac˘a bara se alunge¸ste (F > 0), dimensiunile sale transversale scad x2 − X2 = u2 (X, t) = −

λ F x2 , 2µ(3λ + 2µ)

(4.36)

x3 − X3 = u3 (X, t) = −

λ F x3 . 2µ(3λ + 2µ)

(4.37)

Fie ν astfel ˆıncˆat ε22 = ε33 = −νε11 . Deducem cˇa ν=

(4.38)

λ . 2(λ + µ)

(4.39)

Constanta ν este un coeficient de material cunoscut sub numele de coeficientul lui Poisson. Dac˘a l este diametrul sect¸iunii transversale a grinzii ˆınainte de deformat¸ie iar l + ∆l este diametrul sect¸iunii transversale a grinzii dup˘a deformat¸ie, putem scrie ∆l ∆L = −ν , (4.40) l L de unde putem specifica semnificat¸ia fizic˘a a coeficientului lui Poisson: este un coeficient f˘ar˘a dimensiune iar experimental s-a dovedit cˇa ν > 0.

(4.41)

Au loc a¸sadar urmˇatoarele inegalitˇa¸ti: λ > 0, λ+µ

λ+µ > 0. µ(3λ + 2µ)

S¸tiam ˆıns˘a c˘a 3K = 3λ + 2µ > 0. Prin urmare, λ > 0,

µ>0

¸si, deducem astfel 0 0). λ+µ

(4.42)

Putem exprima ˆın acest moment modulul de rigiditate la compresiune precum ¸si constantele lui Lam´e, λ ¸si µ, ˆın funct¸ie de constantele de material determinate experimental: modulul lui Young E ¸si coeficientul lui Poisson ν, astfel, 3K = 3λ + 2µ =

E , 1 − 2ν

λ=

νE , (1 − 2ν)(1 + ν)

µ=

E 2(1 + ν)

Legea lui Hooke inversat˘a, poate fi scris˘a cu ajutorul lui E ¸si ν astfel,

ε=

1+ν ν σ − (trε)I 3 , E E

sau pe componente

εij =

4.3

1+ν ν σij − σkk δij E E

∀i, j ∈ {1, 2, 3}.

Torsiunea unui arbore cilindric

Cadrul fizic Consider˘am un arbore cilindric care ocup˘a un domeniu Ω ⊂ R3 . Cilindrul este de lungime L ¸si vom nota prin ω o sect¸iune transversal˘a arbitrar˘a a cilindrului ( ω ⊂ R2 domeniu simplu conex); a se vedea Figura 4.3. Presupunem c˘a baza Γ0 este fix˘a ¸si impunem bazei Γ1 o torsiune αL, unde α > 0. Presupunem c˘a densitatea fort¸elor masice este neglijabil˘a, suprafat¸a lateral˘a Γlat este liber˘a de constrˆangeri, iar pe bazele Γ0 ¸si Γ1 impunem eforturi normale nule. Modelare matematic˘a prin E.D.P. Admit¸ˆand cˇa cilindrul este liniar elastic, omogen ¸si izotrop, problema torsiunii se formuleazˇa astfel: ¯ → R, i ∈ {1, 2, 3}¸si σij : Ω ¯ → R, Problema 4.16. Sˇa se determine ui : Ω

Γ1 Ω

x3

L

ω

O x2 x

1

1111111 0000000 0000000 1111111 Γ0 0000000 1111111 0000000 1111111

Figura 4.3: Cadrul fizic i, j ∈ {1, 2, 3}, astfel ˆıncˆat, σij,j = 0,

ˆın Ω

σij = λ(εkk )δij + 2µεij , ˆın Ω σij νj = 0,

pe Γlat

σ3j νj = 0,

pe Γ0

u1 = u2 = 0,

pe Γ0

u1 = −αLx2 ,

pe Γ1

u2 = αLx1 ,

pe Γ1 .

Consider vectorul u de componente     u1 = −αx3 x2     u2 = αx3 x1        u3 = αφ(x1 , x2 ).

∪

Γ1

(4.43)

ˆIn acest caz,  0 0    ε= 0 0    α α ∂φ ∂φ ) ) (−x2 + (x1 + 2 ∂x1 2 ∂x2

∂φ + ) ∂x1 ∂φ α (x1 + ) 2 ∂x2

α (−x2 2

    .   

0

Cilindrul fiind un solid omogen, elastic ¸si izotrop putem g˘asi cu u¸surint¸a˘ pentru tensorul σ urm˘atoarea expresie:   ∂φ 0 0 αµ(−x2 + )  ∂x1      ∂φ  σ= (4.44) )  0 0 αµ(x1 + . ∂x2     ∂φ ∂φ αµ(−x2 + ) αµ(x1 + ) 0 ∂x1 ∂x2 Ecuat¸iile de echilibru conduc ˆın acest caz la urm˘atoarea ecuat¸ie scalar˘a ∆φ = 0 ˆın ω. Pe Γlat avem ν = (ν1 (x1 , x2 ), ν2 (x1 , x2 ), 0). Pentru i = 1, σ11 ν1 + σ12 ν2 + σ13 ν3 = 0, pentru i = 2, σ21 ν1 + σ22 ν2 + σ23 ν3 = 0, pentru i = 3, σ31 ν1 + σ32 ν2 + σ33 ν3 = 0. Deci, µα(−x2 +

∂φ ∂φ )ν1 + µα(x1 + )ν2 = 0, ∂x1 ∂x2

pe ∂ω

¸si, µα(x1 ν2 − x2 ν1 ) + µα(

∂φ ∂φ ν1 + ν2 ) = 0 pe ∂ω. ∂x1 ∂x2

Notˇam cˇa

∂φ ∂φ ∂φ ν1 + ν2 = ∇φ · (ν1 , ν2 ) = . ∂x1 ∂x2 ∂ν

Astfel, ∂φ = x2 ν1 − x1 ν2 , pe ∂ω. ∂ν Deoarece pe Γ0 avem ν = (0, 0, −1), atunci σ31 ν1 + σ32 ν2 + σ33 ν3 = 0 echivalent cu 0 = 0. De asemenea, pe Γ1 avem ν = (0, 0, 1). Atunci, σ31 ν1 + σ32 ν2 + σ33 ν3 = 0 echivalent cu 0 = 0. Prin urmare, funct¸ia φ este solut¸ie a problemei urmˇatoare. Problema 4.17. Sˇa se determine φ : ω ¯ → R astfel ˆıncˆat     ∆φ = 0 ˆın ω (N )    ∂φ = x2 ν1 − x1 ν2 pe ∂ω. ∂ν Aceastˇa problemˇa are solut¸ie unicˇa pˆanˇa la o constantˇa aditiv˘a. A¸sadar, solut¸ia Problemei 4.16 este datˇa prin (4.43) ¸si (4.44), unde φ este solut¸ia unicˇa pˆanˇa la o constantˇa aditivˇa a Problemei 4.17. ˆIn cele ce urmeazˇa suntem interesat¸i sˇa determinˇam zona primelor puncte ale cilindrului ˆın care materialul ˆınceteazˇa sˇa mai fie elastic. Acest studiu necesitˇa cˆateva preliminarii.

Tensiuni normale principale. Direct¸ii principale Fie σ tensorul tensiune al lui Cauchy. Cum σ este tensor de ordinul doi simetric, el are trei valori proprii reale, σI ≤ σII ≤ σIII .

Valorile σj , j ∈ {I, II, III}, se numesc tensiuni normale principale. Fie ν j , j ∈ {I, II, III}, vectorii proprii corespunz˘atori, ortogonali doi cˆate doi. Ace¸stia se numesc direct¸ii principale de tensiune. ˆIntr-un reper local cu axe de direct¸ii ν j , σ are urm˘atoarea reprezentare 

 0  σI 0   σ= 0  0 σII   0 0 σIII

   .   

Cercurile lui Mohr. Fort¸e de forfecare Fie n = n1 e1 +n2 e2 +n3 e3 astfel ˆıncˆat |n| = R3 fiind baza canonic˘a a lui R3 .

√ n21 + n22 + n23 = 1, {e1 , e2 , e3 } ⊆

Sˇa considerˇam vectorul F = (F1 , F2 , F3 ), unde F1 = σ I n1 , F2 = σ II n2 , F3 = σIII n3 . Acest vector se descompune astfel : F = FN · n + F T , unde FN = F · n. Avem |F |2 = F · F = (FN · n + F T ) · (FN · n + F T ) = |FN |2 + |F T |2 . Prin urmare, |F T |2 = |F |2 − |FN |2 . ˆIn cele ce urmeaz˘a studiem cum variaz˘a |F T | ˆın funct¸ie de FN , atunci cˆand n variaz˘a. Pentru ˆınceput putem scrie     n21 + n22 + n23 = 1;     σI n21 + σII n22 + σIII n23 = FN ;       2 2 2  σI2 n21 + σII n23 = |F T |2 + FN2 . n2 + σIII

Fie polinomul P (x) = x2 + ax + b. Avem 2 2 σI2 n21 + aσI n21 + bn21 + σII n2 + aσII n22 + bn22 2 +σIII n23 + aσIII n23 + bn23

= |F T |2 + FN2 + aFN + b, ¸si de aici n21 P (σI ) + n22 P (σII ) + n23 P (σIII ) = |F T |2 + FN2 + aFN + b. Fie P (x) = (x − σI )(x − σII ). Atunci, P (σI ) = 0 ¸si P (σII ) = 0 ¸si n23 (σIII − σI )(σIII − σII ) = |F T |2 + FN2 − (σI + σII )FN + σI σII . Fie P (x) = (x − σII )(x − σIII ). ˆIn acest caz, P (σII ) = 0 ¸si P (σIII ) = 0 ¸si n21 (σI − σII )(σI − σIII ) = |F T |2 + FN2 − (σII + σIII )FN + σII σIII . Fie P (x) = (x − σI )(x − σIII ). Avem acum P (σI ) = 0 ¸si P (σIII ) = 0 ¸si n22 (σIII − σI )(σIII − σII ) = |F T |2 + FN2 − (σI + σIII )FN + σI σIII . A¸sadar, n23 (σIII − σI )(σIII − σII ) = |F T |2 + (FN − σI )(FN − σII );

n21 (σI − σIII )(σI − σIII ) = |F T |2 + (FN − σII )(FN − σIII ); n22 (σII − σI )(σII − σIII ) = |F T |2 + (FN − σI )(FN − σIII ). Cazul 1. Dac˘a σI < σII < σIII , atunci n23 =

|F T |2 + (FN − σI )(FN − σII ) , (σIII − σI )(σIII − σI )

n21 =

|F T |2 + (FN − σII )(FN − σIII ) , (σI − σII )(σI − σIII )

n22 =

|F T |2 + (FN − σI )(FN − σIII ) , (σII − σI )(σII − σIII )

¸si de aici, |F T |2 + (FN − σI )(FN − σII ) ≥ 0, |F T |2 + (FN − σII )(FN − σIII ) ≥ 0, |F T |2 + (FN − σI )(FN − σIII ) ≥ 0. Dac˘a not˘am |F T | = Y ¸si FN = X atunci (

(σ − σ ) σI + σII )2 II I ) +Y2 ≥ ) ; 2 2 ( (σ − σ ) σII + σIII )2 III I X− ) +Y2 ≥ ) ; 2 2 ( (σ − σ ) σI + σIII )2 III I X− ) +Y2 ≥ ) . 2 2 X−

Din Figura 4.3 se observˇa u¸sor cˇa Y are valoare maximˇa pentru X = σIII − σI . 2 σIII − σI |F T |max := 2 se nume¸ste forfecare maximal˘a. Forfecarea maximal˘a este valoarea la care se produc rupturi ˆın material. Cazul 2. Dac˘a σI = σII ≤ σIII , atunci n23 =

|F T |2 + (FN − σI )2 . (σIII − σI )2

Prin urmare, ˆın acest caz Y 2 + (X − σI )2 ≥ 0, adicˇa (X, Y ) apart¸ine σIII − σI , 0). Deducem imediat cˇa Y pˇart¸ii pozitive a discului centrat ˆın ( 2 σIII − σI are valoare maximˇa pentru X = . 2 Cazul 3. σI = σII = σIII ˆIn acest caz avem σI2 (n21 +n22 +n23 ) = FN2 +|F T |2 ¸si σI (n21 + n22 + n23 ) = FN , astfel cˇa σI2 = FN2 + |F T |2 ¸si σI2 = FN2 . Din ultimele dou˘a relat¸ii rezult˘a c˘a |F T | = 0. Deci, X = σI ¸si Y = 0. Valoarea maximˇa σIII − σI . a lui Y este 0 = 2

Calculul forfecˇ arii maximale Revenim acum asupra Problemei 4.16. Cum Div σ = 0, ∂σ13 ∂σ23 + = 0 ˆın ω. ∂x1 ∂x2 Deci, ∂σ23 ∂σ13 =− ∂x1 ∂x2

ˆın ω.

Rezult˘a c˘a exist˘a θ = θ(x1 , x2 ) astfel ˆıncˆat σ13 =

∂θ , ∂x2

¸si σ23 = −

∂θ . ∂x1

A¸sadar, µα(−x2 +

∂φ ∂θ )= ∂x1 ∂x2

¸si µα(x1 +

∂φ ∂θ )=− . ∂x2 ∂x1

Deducem cˇa, ∂ 2θ ∂ 2θ ∂ 2φ ∂ 2φ + = −µα + µα − µα − µα . ∂x22 ∂x21 ∂x1 ∂x2 ∂x2 ∂x1 Deci, pentru φ funct¸ie suficient de regulatˇa, ∆θ = −2µα,

ˆın ω.

Pe de altˇa parte, σ31 ν1 + σ32 ν2 = 0,

pe ∂ω

sau, echivalent, ∂θ ∂θ ν1 − ν2 = 0, ∂x2 ∂x1 Fie τ = (τ1 , τ2 ) astfel ˆıncˆat τ1 = −ν2

¸si

pe

τ2 = ν1 .

Atunci, pe ∂ω avem ∂θ ∂θ τ1 + τ2 = 0, ∂x1 ∂x2 sau echivalent, ∂θ =0 ∂τ

pe

∂ω.

∂ω.

Deducem cˇa ∇θ = 0 pe ∂ω ¸si prin urmare θ este constant pe ∂ω. Alegem θ=0

pe

∂ω

¸si considerˇam urmˇatoarea problemˇa. Problema 4.18. Sˇa se determine θ : ω ¯ → R astfel ˆıncˆat     ∆θ = −2µα ˆın ω (D)    θ=0 pe ∂ω. Problema 4.18 are solut¸ie unicˇa. Vom vedea cˇa fort¸a de forfecare maximalˇa νIII − νI |F T |max = 2 poate fi exprimatˇa prin intermediul solut¸iei Problemei 4.18. Sˇa evaluˇam tensiunile normale principale νI ¸si νIII . √ √ 2 2 2 2 + ν23 νI = − ν13 + ν23 ; νIII = ν13 Prin urmare, fort¸a de forfecare maximalˇa are expresia √ ∂θ 2 ∂θ 2 |F T |max = + = | grad θ|, ∂x2 ∂x1 unde θ este unica solut¸ie a Problemei 4.18. Se ¸stie cˇa | grad θ| este funct¸ie subarmonic˘a, a se vedea de exemplu [30]. Conform principiului de maxim pentru funct¸ii subarmonice, deducem c˘a aceast˘a funct¸ie ˆı¸si atinge maximul pe ∂ω. De aceea, primele puncte ˆın care materialul va ˆınceta s˘a fie elastic sunt situate pe Γlat .

Capitolul 5 Solut¸ii slabe. Aproximarea solut¸iilor slabe

5.1

Probleme antiplane

Din punct de vedere matematic, problemele antiplane ˆın deplasˇari-tract¸iuni sunt problemˇa la limitˇa Neumann-Dirichlet. Pentru materiale liniar elastice, omogene ¸si izotrope, avem urmˇatoarea problemˇa la limitˇa. Problema 5.19. Fie f0 : Ω → R, f2 : Γ2 → R ¸si µ > 0. Sˇa se determine ¯ → R astfel ˆıncˆat u:Ω µ∆u(x, t) + f0 = 0

in Ω

(5.1)

u=0 ∂u µ = f2 ∂ν

pe Γ1

(5.2)

pe Γ2 .

(5.3)

93

¯ → R sufiPrin solut¸ie tare a Problemei 5.19, ˆınt¸elegem o funct¸ie u : Ω ¯ care verificˇa Procient de regulatˇa ( de exemplu de clasˇa C 2 (Ω) ∩ C 1 (Ω)) blema 5.19. Dacˇa datele f0 ¸si f2 au expresii complicate, este uneori dificil, alteori chiar imposibil sˇa determinˇam explicit solut¸ia tare. De aceea este util sˇa definim pentru Problema 5.19 not¸iunea de solut¸ie slabˇa, interesˆandu-ne apoi de o bunˇa aproximare a acesteia. Admitem ˆın aceastˇa sect¸iune urmˇatoarea ipotezˇa asupra datelor f0 ¸si f2 , f0 ∈ L2 (Ω),

f2 ∈ L2 (Γ2 ).

(5.4)

Sˇa considerˇam pentru ˆınceput cˇa u este o funct¸ie suficient de regulatˇa. ¯ ¸si utilizˆand prima formulˇa a lui Multiplicˆand (5.1) cu o funct¸ie v ∈ C ∞ (Ω) Green, obt¸inem ∫ ∫ ∫ ∂u µ∇u · ∇v dx = f0 v dx + µ v dx. Ω Ω Γ ∂ν ¯ Folosind acum (5.3), putem scrie pentru orice v ∈ C ∞ (Ω), ∫ ∫ ∫ ∫ ∂u µ∇u · ∇v dx = f0 v dx + µ v dx + f2 v dx. Ω Ω Γ1 ∂ν Γ2 ¯ = H 1 (Ω), Deoarece C ∞ (Ω) ∫ ∫ ∫ ∫ ∂u µ∇u · ∇v dx = f0 v dx + µ γv dx + f2 γv dx ∀v ∈ H 1 (Ω). ∂ν Ω Ω Γ1 Γ2 Cum V este un subspat¸iu al lui H 1 (Ω), iar ∫ ∂u µ γv dx = 0 ∀v ∈ V, Γ1 ∂ν obt¸inem

∫

∫ µ∇u · ∇v dx = Ω

∫ f2 γv dx ∀v ∈ V.

f0 v dx + Ω

Γ2

Observˇam cˇa dacˇa u ∈ V atunci ∇u ∈ L2 (Ω)2 . Cum ∇v ∈ L2 (Ω)2 ¸si µ > 0, µ∇u · ∇v ∈ L1 (Ω). Notˇam ˆın plus cˇa ipoteza (5.4) combinatˇa cu v ∈ V

implicˇa f0 v ∈ L1 (Ω) ¸si f2 γv ∈ L1 (Γ2 ), deci integralele din membrul drept sunt de asemenea bine definite. Aceasta ne permite sˇa definim ∫ ∫ ϕ : V → R ϕ(v) = f0 v dx + f2 γv dx. Ω

Γ2

Deoarece γ este operator liniar, utilizˆand ¸si liniaritatea integralelor, deducem cˇa ϕ este aplicat¸ie liniarˇa. ˆIn plus, ϕ este aplicat¸ie continuˇa. Pentru a justifica continuitatea, arˇat cˇa existˇa M > 0 astfel ˆıncˆat |ϕ(v)| ≤ M ∥v∥V

∀v ∈ V.

ˆIntr-adevˇar, |ϕ(v)| ≤ ∥f0 ∥L2 (Ω) ∥v∥L2 (Ω) + ∥f2 ∥L2 (Γ2 ) ∥γv∥L2 (Γ2 ) ≤ ∥f0 ∥L2 (Ω) ∥v∥H 1 (Ω) + ∥f2 ∥L2 (Γ2 ) ∥γv∥L2 (Γ) . Cum γ este operator liniar ¸si continuu, existˇa c0 > 0 astfel ˆıncˆat ∥γv∥L2 (Γ) ≤ c0 ∥v∥H 1 (Ω) . Notˇam de asemenea cˇa existˇa c1 , c2 > 0 astfel ˆıncˆat c1 ∥v∥V ≤ ∥v∥H 1 (Ω) ≤ c2 ∥v∥V

∀v ∈ V.

Prin urmare, putem considera M = c2 (∥f0 ∥L2 (Ω) + c0 ∥f2 ∥L2 (Γ2 ) ). Utilizˆand acum teorema de reprezentare a lui Riesz, putem defini f ∈ V, dupˇa cum urmeazˇa, (f, v)V := (f0 , v)L2 (Ω) + (f2 , γv)L2 (Γ2 ) . ˆIn plus, definim a : V × V → R forma biliniarˇa, astfel ˆıncˆat ∫ a(u, v) = µ∇u · ∇v dx u, v ∈ V. Ω

Suntem condu¸si astfel la urmˇatoarea formulare variat¸ionalˇa a Problemei 5.19. Problema 5.20. Dat f ∈ V, sˇa se determine u ∈ V astfel ˆıncˆat a(u, v) = (f, v)V

∀v ∈ V.

(5.5)

Definit¸ia 5.21. Problemei 5.20.

Numim solut¸ie slabˇa a Problemei 5.19 orice solut¸ie a

Teorema 5.22. astfel ˆıncˆat

ˆ plus, existˇa C > 0 Problema 5.20 are solut¸ie unicˇa. In ∥u∥V ≤ C∥f ∥V

∀v ∈ V

(5.6)

¸si ∥u1 − u2 ∥V ≤ C∥f 1 − f 2 ∥V

(5.7)

unde u1 , u2 ∈ V sunt solut¸iile corespunzˇ atoare datelor f 1 , f 2 ∈ V. Demonstrat¸ie : Pentru a justifica existent¸a ¸si unicitatea solut¸iei Problemei 5.20, verificˇam ipotezele Teoremei Lax-Milgram, a se vedea Teorema 1.4, plasˆandu-ne ˆın spat¸iul Hilbert V. |a(u, v)| ≤ µ∥∇u∥L2 (Ω)2 ∥∇v∥L2 (Ω)2 = µ∥u∥V ∥u∥V . Observˇam cˇa (1.18) se verificˇa pentru M = µ. ˆIn plus, a(u, u) = µ∥∇u∥2L2 (Γ2 ) = µ∥u∥2V . Deci, (1.19) se verificˇa pentru m = µ. Considerˇam ˆın (5.5) v = u astfel cˇa a(u, u) = (f, u)V . Cum a este formˇa V -elipticˇa, µ∥u∥2V ≤ |(f, u)V |.

Aplicˇam acum inegalitatea Cauchy-Schwarz obt¸inˆand, µ∥u∥2V ≤ ∥f ∥V ∥u∥V . Dacˇa ∥u∥V ̸= 0, ˆımpˇart¸ind inegalitatea precedentˇa prin ∥u∥V obt¸inem (5.6) cu C = µ1 . Dacˇa ∥u∥V = 0, (5.6) este de asemenea adevˇaratˇa. Fie f 1 , f 2 ∈ V. a(u1 , v) = (f 1 , v)V

∀v ∈ V,

a(u2 , v) = (f 2 , v)V

∀v ∈ V.

Alegem v = u2 − u1 ¸si sumˇam cele douˇa identitˇa¸ti membru cu membru obt¸inˆand a(u1 − u2 , u2 − u1 ) = (f 1 − f 2 , u2 − u1 )V ¸si de aici, a(u1 − u2 , u1 − u2 ) = (f 1 − f 2 , u1 − u2 )V . Utilizˆand V −elipticitatea formei a ¸si inegalitatea Cauchy-Schwarz obt¸inem µ∥u1 − u2 ∥2V ≤ ∥f 1 − f 2 ∥V ∥u1 − u2 ∥V care ne conduce imediat la (5.7). Putem considera de exemplu C = µ1 .



Fie Vh subspat¸iu de dimensiune finitˇa al spat¸iului V. Considerˇam problema Problema 5.23. Dat f ∈ V, sˇa se determine uh ∈ Vh astfel ˆıncˆat a(uh , vh ) = (f, vh )V

∀uh ∈ Vh .

Problema 5.23 are solut¸ie unicˇa, aceasta rezultˆand tot pe baza Teoremei Lax-Milgram 1.4. Solut¸ia Problemei 5.23 aproximeazˇa solut¸ia Problemei 5.20. Ne propunem ˆın cele ce urmeazˇa sˇa aproximˇam solut¸ia slabˇa a Problemei ¯1 ∪ Γ ¯ 2, 5.19 ˆın urmˇatorul caz particular. Fie Ω = (1, 2) × (0, 1) ¸si ∂Ω = Γ = Γ Γ1 ∩ Γ2 = ∅, unde Γ1 = (AB] ∪ (BC),

Γ2 = [DC] ∪ [AD),

f0 : Ω → R f0 (x, y) = 2(x − 2)2 + 2y 2 ,     2y 2 pe [AD) f2 : Γ2 → R f2 (x, y) =    2(x − 2)2 pe [DC] ¸si µ = 1. Notˇam cˇa f0 ∈ L2 (Ω) ¸si f2 ∈ L2 (Γ2 ). Pentru acest exemplu se poate verifica prin calcul direct cˇa funct¸ia ¯ → R u(x, y) = (x − 2)2 y 2 u:Ω este solut¸ie exactˇa. Trasarea graficului solut¸iei exacte se face ˆın Matlab astfel: x=1:0.09:2; y=0:0.09:1; [X,Y]=meshgrid(x,y); U1=Y.^2; U2=(X-2).^2; U=U1.*U2; surf(X,Y,U). Ne propunem sˇa trasˇam graficul solut¸iei exacte ¸si sa-l vizualizˇam ˆıntr-o fereastrˇa paralelˇa cu fereastra ce va prezenta graficul solut¸iei slabe aproximative. Pentru a aproxima solut¸ia slabˇa vom utiliza metoda elementului finit, iar implementarea se va realiza ˆın Matlab. Pentru aceasta, discretizˇam domeniul utilizˆand 32 elemente finite (elemente triunghiulare cu trei noduri), ca ˆın Figura 5.1. Pentru detalii privind triangulat¸ia unui domeniu, a se consulta de exemplu [26]. Numˇarul total de noduri este 25. Matlab-ul va citi datele problemei ˆın format ASCII prin fi¸siere cu extensia .dat. ˆIn Sect¸iunea 6.2 sunt prezentate tablourile de date utilizate: coordinates.dat, elementes.dat, dirichlet.dat ¸si neumann.dat. Calculul solut¸iei aproximative se reduce la rezolvarea urmˇatorului sistem liniar, Au = b,

x2

C(2,1)

D(1,1) 21

22

32

31 19

18

24

23

22 18

17

16

19

20 14

15

14

13

10

11

9

8

8

7

6

1

2 2

7

15

12

6

5 3

3

16

21

13

9

1

28 17

12

10

25

29 27

20

11

24

30 26

25

O

23

4 4

5

x1

B(2,0)

A(1,0)

Figura 5.1: Discretizarea domeniului Ω = (1, 2) × (0, 1) ∫

unde A = (Aij );

∇ηi · ∇ηj dx i, j ∈ 1, 25

Aij = Ω

∫

¸si b = (bi );

bi =

∫ f2 ηi ds i ∈ 1, 25.

f0 ηi dx + Ω

Γ2

Deoarece discretizarea s-a fˇacut prin intermediul a 25 de noduri, spat¸iul Vh are dimensiunea 25. Drept bazˇa considerˇam mult¸imea elementelor ηj , j ∈ 1, 25, de forma urmˇatoare ηj (x, y) =

1 [xk ym − xm yk + (yk − ym )x + (xm − xk )y] 2|T |

pentru oricare (x, y) ∈ T, unde T este un triunghi de vˆarfuri (xj , yj ), (xk , yk ) ¸si (xm , ym ) iar |T | este aria triunghiului. Notˇam cˇa ηj (xj , yj ) = 1 ˆın timp ce ηj (xk , yk ) = 0 ¸si ηj (xm , ym ) = 0. Fie T un triunghi arbitrar din triangulat¸ia domeniului Ω. Definim matricea M de componente ∫ Mmn = ∇ηm · ∇ηn dx m, n ∈ 1, 3. T

solutia exacta

1

1

0.8

0.8

0.6

0.6

axa 0z

axa 0z

solutia aproximativa

0.4

0.4

0.2

0.2

0 1

0 1 2 0.5

axa 0y

1.5 0 1

axa 0x

2 0.5 axa 0y

Figura 5.2: Exemplul 1

1.5 0 1

axa 0x

Matricea pˇatraticˇa de ordinul 3, M = (Mmn ), se nume¸ste matrice elementalˇa. Matricea A se nume¸ste matrice de rigiditate ¸si se calculeazˇa prin intermediul matricilor elementale, printr-o procedurˇa de asamblare dupˇa cum urmeazˇa. • Init¸ializˇam cu zero toate componentele matricei A; • Asamblˇam matricea elementalˇa corespunzˇatoare primului element finit (primul element finit este triunghiul care are drept vˆarfuri nodurile 1, 2 ¸si 9 din ret¸ea; a se vedea Figura 5.1; A11 = A11 + M11 , A12 = A12 + M12 , A19 = A19 + M13 , A21 = A21 + M21 , A22 = A22 + M22 , A29 = A29 + M23 , A91 = A91 + M31 , A92 = A92 + M32 , A99 = A99 + M33 , restul elementelor matricei A rˇamˆanˆand neschimbate; • Asamblˇam matricea elementalˇa corespunzˇ atoare celui de al doilea element finit (al doilea element finit este triunghiul care are drept vˆarfuri nodurile 2, 3 ¸si 8 din ret¸ea; vezi Figura 5.1; A22 = A22 + M11 , A23 = A23 + M12 , A28 = A28 + M13 , A32 = A32 + M21 , A33 = A33 + M22 , A38 = A38 + M23 , A82 = A82 + M31 , A83 = A83 + M32 , A88 = A88 + M33 , restul elementelor rˇamˆanˆand neschimbate, ¸s.a.m.d. ˆIn MATLAB, aceastˇa asamblare se face dupˇa cum urmeazˇa: for j = 1:size(elementes,1) A(elementes(j,:),elementes(j,:)) = ... A(elementes(j,:),elementes(j,:))... + melementala(coordinates(elementes(j,:),:)); end

Calculul fort¸elor volumice elementale (vectori de 3 componente) se face astfel, ∫ ∫ f0 ηm ds ≈ f0 (xM , yM )ηm (xM , yM ) ds = f0 (xM , yM )ηm (xM , yM )|T | T

T

pentru fiecare m ∈ 1, 3, unde (xM , yM ) este centrul de greutate al triunghiului T. Asamblarea fort¸elor volumice elementale se face ˆın MATLAB astfel: for j=1:size(elementes,1) b(elementes(j,:))=b(elementes(j,:)) + ... det([1 1 1 ; coordinates(elementes(j,:),:)’])*... f_0(sum(coordinates(elementes(j,:),:))/3)/6; end Fie E o laturˇa a unui triunghi al triangulat¸iei astfel ˆıncˆat E ⊂ Γ2 . Calculul tract¸iunilor elementale (vectori cu douˇa componente) se face astfel, ∫ ∫ f2 ηm ds ≈ f2 (xM , yM )ηm (xM , yM ) ds = f2 (xM , yM )ηm (xM , yM )|E|. E

E

pentru fiecare m ∈ 1, 2 (ˆın numerotare localˇa). Asamblarea tract¸iunilor elementale ˆıntr-un vector cu 7 componente (deoarece 7 noduri sunt implicate ˆın condit¸ia Neumann, celelalte 9 noduri fiind noduri implicate ˆın condit¸ia Dirichlet) se realizeazˇa astfel: for j = 1 : size(neumann,1) b(neumann(j,:)) = b(neumann(j,:)) + ... norm(coordinates(neumann(j,1),:)-... coordinates(neumann(j,2),:))*... f_2(sum(coordinates(neumann(j,:),:))/2)/2; end Potrivit condit¸iei Dirichlet, deja ¸stim valoarea funct¸iei u ˆın cele 9 noduri de pe Γ1 . Celelalte 16 noduri le vom introduce ˆıntr-un vector numit

solutia aproximativa pentru f,f/100,f/1000 0.8 0.7

1 0.8

0.6

0.6

0.5

0.4

0.4

0.2

0.3

0 1

0.2 2 0.1

0.5

1.5 0

1

Figura 5.3: Test Exemplul 1

FreeNodes. Prin urmare vom rezolva un sistem de numai 16 ecuat¸ii cu 16 necunoscute. Scrierea ˆın MATLAB se va realiza astfel: u = sparse(size(coordinates,1),1); u(unique(dirichlet))= ... u_d(coordinates(unique(dirichlet),:)); b = b-A*u; u(FreeNodes) = A(FreeNodes,FreeNodes)\b(FreeNodes) Codul matlab complet care are drept rezultat trasarea ˆın ferestre paralele a graficului solut¸iei exacte ¸si a graficului solut¸iei aproximative poate fi gˇasit ˆın Anexˇa, Sect¸iunea 6.1.

Putem compara graficul solut¸iei exacte cu graficul solut¸iei aproximative, testˆand astfel eficient¸a algoritmului propus; a se vedea Figura 5.2. Din (5.6) deducem cˇa dacˇa f = 0V atunci solut¸ia Problemei 5.20 este solut¸ia nulˇa. Mai mult, din (5.7) deducem cˇa dacˇa f se diminueazˇa(ˆın normˇa) atunci ¸si deplasarea se diminueazˇa. Testˇam numeric ˆın cele ce urmeazˇa aceastˇa ultimˇa afirmat¸ie. ˆIn Figura 5.3 reprezentˇam ˆın acela¸si sistem de axe solut¸iile aproximative corespunzˇatoare datelor: • f0 ¸si f2 ; • f0 /100 ¸si f2 /100; • f0 /1000 ¸si f2 /1000. Sesizˇam apropierea graficului solut¸iei aproximative de planul Ox1 x2 pe mˇasurˇa ce fort¸ele se diminueazˇa.

5.2

Probleme plane

ˆIn aceastˇa sect¸iune sunt studiate probleme plane ˆın deplasˇari-tract¸iuni, pentru materiale liniar elastice, omogene ¸si izotrope. Problema 5.24. Sˇa se determine u : Ω ⊂ R2 → R2 ¸si σ : Ω → S2 astfel ˆıncˆat Div σ + f 0 = 0

ˆın Ω,

(5.8)

σ = Cε(u)

ˆın Ω,

(5.9)

u = w

pe Γ1 ,

(5.10)

σ n = f2

pe Γ2 ,

(5.11)

unde λ, µ > 0 ¸si Cε(u) = λ(trε(u))I 3 + 2µε(u).

¯ → Prin solut¸ie tare a Problemei 5.24 ˆınt¸elegem o pereche (u, σ), u : Ω ¯ → S2 , u ∈ C 2 (Ω) ∩ C 1 (Ω), ¯ σ ∈ C 1 (Ω), care verificˇa (5.8)-(5.11). R2 , σ : Ω Admitem ˆın cele ce urmeazˇa cˇa densitatea fort¸elor volumice ¸si densitatea tract¸iunilor verificˇa ipoteza f 0 ∈ L2 (Ω)2 ,

f 2 ∈ L2 (Γ2 )2 .

(5.12)

Fˇarˇa a restrˆange generalitatea, vom presupune ˆın cele ce urmeazˇa cˇa ˜ ∈ H1 astfel ˆıncˆat w = γ w ˜ ¸si apoi w = 0 (altfel, presupunem cˇa existˇa w ˜ facem translat¸ia u := u − w). Ne propunem sˇa introducem not¸iunea de solut¸ie slabˇa pentru Problema 5.24. Fie (u, σ) solut¸ie tare a Problemei 5.24. Utilizˆand (5.8), deducem (Div σ, v)H + (f 0 , v)H = 0 ∀v ∈ H1 . Folosind formula lui Green (1.12), putem scrie ∫ ∫ (σ, ε(v))H = (f 0 , v)H + σn · γv dΓ + σn · γv dΓ ∀v ∈ H1 . Γ1

Γ2

Fie V spat¸iul de funct¸ii vectoriale definit ˆın (1.13). Utilizˆand (5.11) obt¸inem ∫ (σ, ε(v))H = (f 0 , v)H + σn · γv dΓ ∀v ∈ V. Γ2

Sˇa folosim acum (5.9). Avem astfel,

∫ f 2 · γv dΓ.

(Cε(u), ε(v))H = (f 0 , v)H + Γ2

Utilizˆand (5.12) definim

∫

Φ : V → R Φ(v) =

∫ f 0 · v dx +

Ω

f 2 · γv dx. Γ2

Deoarece γ este operator liniar, utilizˆand ¸si liniaritatea integralelor, deducem cˇa Φ este aplicat¸ie liniarˇa. ˆIn plus, utilizˆand continuitatea operatorului urmˇa γ, deducem cˇa Φ este aplicat¸ie continuˇa. Aceste proprietˇa¸ti ale

funct¸iei Φ ne permit sˇa aplicˇam teorema de reprezentare a lui Riesz care ne furnizeazˇa un unic element f ∈ V astfel ˆıncˆat ∫ ∫ (f , v)V = f 0 · v dx + f 2 · γv da ∀v ∈ V. (5.13) Ω

Γ2

Definim acum forma bilinearˇa ∫ a : V × V → R a(u, v) =

Cε(u) · ε(v) dx,

∀v ∈ V.

(5.14)

Ω

Astfel, Problema 5.24 are urmˇatoarea formulare variat¸ionalˇa. Problema 5.25. Dat f ∈ V, sˇa se determine u ∈ V astfel ˆıncˆat a(u, v) = (f , v)V Definit¸ia 5.26. Problemei 5.25.

∀v ∈ V.

(5.15)

Numim solut¸ie slabˇa a Problemei 5.24 orice solut¸ie a

Teorema 5.27. Admitem (5.12). Atunci, Problema 5.25 are solut¸ie unicˇa. ˆ plus, existˇa C > 0 astfel ˆıncˆ In at ∥u∥V ≤ C∥f ∥V

∀v ∈ V

(5.16)

¸si ∥u1 − u2 ∥V ≤ C∥f 1 − f 2 ∥V

(5.17)

unde u1 , u2 ∈ V sunt solut¸iile corespunzˇ atoare datelor f 1 , f 2 ∈ V. Demonstrat¸ie :

Folosind inegalitatea Cauchy-Schwarz, putem scrie |a(u, v)| ≤ ∥Cε(u)∥H ∥ε(v)∥H .

Deoarece se poate verifica u¸sor cˇa existˇa C > 0 astfel ˆıncˆat ∥Cε(u)∥2H ≤ C∥ε(u)∥2H ,

folosind definit¸ia produsului scalar pe V, a se vedea (1.15), deducem cˇa forma a verificˇa (1.18) pentru oricare u, v ∈ V. Sˇa evaluˇam acum a(u, u) pentru u ∈ V. a(u, u) ≥ 2µ∥ε(u)∥2H = 2µ∥u∥2V . Prin urmare, forma biliniarˇa a verificˇa ¸si (1.19). Pentru a demonstra existent¸a ¸si unicitatea solut¸iei Problemei 5.25 aplicˇam Teorema 1.4 plasˆandu-ne ˆın spat¸iul X = V. Considerˇam ˆın (5.15) v = u astfel cˇa a(u, u) = (f , u)V . Inegalitatea (5.22) rezultˇa de aici utilizˆand inegalitatea Cauchy-Schwartz ¸si V−elipticitatea formei a. Fie f 1 , f 2 ∈ V. a(u1 , v) = (f 1 , v)V

∀v ∈ V,

a(u2 , v) = (f 2 , v)V

∀v ∈ V.

Alegem v = u2 − u1 ¸si sumˇam cele douˇa identitˇa¸ti membru cu membru obt¸inˆand a(u1 − u2 , u1 − u2 ) = (f 1 − f 2 , u1 − u2 )V , iar de aici, µ∥u1 − u2 ∥2V ≤ ∥f 1 − f 2 ∥V ∥u∥V . Considerˆand de exemplu C = µ1 , obt¸inem (5.23).



ˆIn cazul problemelor plane putem scrie: ( µ2 1) 2 − σ11 σ22 ) σ D |2 = + (σ11 + σ22 )2 + 2(σ12 6(µ + λ) 2 unde |σ D | este norma deviatorului tensorului σ, a se vedea (1.6). Notˇam cˇa |σ D |2 4µ

3.5

3

2.5

2

1.5

1

0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

Figura 5.4: Exemplul 2: discretizarea domeniului

−3

x 10

4

4

3.5

3.5

3

3 2.5 2.5 2 2 1.5 1.5 1 1 0.5 0

0.5 0

0.5

1

1.5

2

2.5

3

Figura 5.5: Exemplul 2

3.5

−3

x 10

3.5

3.5

3

3

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0

0.5

0

0.5

1

1.5

2

2.5

3

3.5

Figura 5.6: Test Exemplul 2

are semnificat¸ie de densitate de energie elasticˇa. Fie Vh un subspat¸iu de dimensiune 2N al spat¸iului V. Considerˇam baza {η 1 , ..., η 2N } = {η1 e1 , η1 e2 , ..., ηN e1 , ηN e2 } unde N este numˇarul nodurilor din ret¸eaua care discretizeazˇa domeniul. Considerˇam problema discretˇa ata¸satˇa Problemei 5.25. Problema 5.28. Dat f ∈ V, sˇa se determine uh ∈ Vh astfel ˆıncˆat a(uh , v h )V = (f , v h )V

∀uh ∈ Vh .

Pe baza Teoremei 1.4 (Lax-Milgram), Problema 5.28 are o unicˇa solut¸ie care se poate calcula rezolvˆand sistemul algebric liniar Ax = b, unde

∫ A = (Ai,j );

ε(η k ) · Cε(η l )dx k, l ∈ 1, 2N

Akl = Ω

(5.18)

∫

¸si b = (bk );

∫ f 2 · η k dx +

bk = Ω

g · η k ds k ∈ 1, 2N .

(5.19)

ΓN

Sˇa presupunem cˇa domeniul Ω are frontierˇa poligonalˇa, astfel cˇa ¯ = ∪T ∈T T Ω unde T desemneazˇa o triangulat¸ie iar T un element finit al triangulat¸iei. Vom considera T ca fiind un triunghi cu trei noduri sau un pˇatrat. Integralele din (5.18) ¸si (5.19) pot fi calculate pentru fiecare j, k ∈ 1, ..., 2N astfel: ∑∫ Ajk = ε(η j ) · Cε(η k )dx j, k ∈ 1, 2N T ∈T

¸si bj =

∑∫ T ∈T

T

T

f 2 · η j dx +

∑ ∫ E⊂ΓN

g · η j ds j ∈ 1, 2N . E

Deci, putem calcula matricea A ¸si vectorul b prin asamblarea matricelor elementale ¸si a vectorilor elementali. Un calcul laborios conduce la urmˇatoarea expresie pentru matricea elementalˇa, aceea¸si pentru fiecare element triunghiular cu trei noduri: M = |T |RT CR

(5.20)

unde T este un element finit triunghiular din T , |T | aria triunghiului T,   λ 0   λ + 2µ      C= λ λ + 2µ 0       0 0 µ iar

 0 ηk2 ,x 0 ηk3 ,x 0  ηk1 ,x   R= ηk1 ,y 0 ηk2 ,x 0 ηk3 ,x  0   ηk1 ,y ηk1 ,x ηk2 ,y ηk2 ,x ηk3 ,y ηk3 ,x

    ,   

unde k1 , k2 , k3 sunt nodurile din ret¸ea care reprezintˇa vˆarfurile triunghilui T. Implementˇam ˆın Matlab (5.20) dupˇa cum urmeazˇa: function M=melementala3(vertices,lambda,mu) PhiGrad=[1,1,1;vertices’]\[zeros(1,2);eye(2)]; R=zeros(3,6); R([1,3],[1,3,5])=PhiGrad’; R([3,2],[2,4,6])=PhiGrad’; C=mu*[2,0,0;0,2,0;0,0,1] +lambda*[1,1,0;1,1,0;0,0,0]; M=det([1,1,1;vertices’])/2*R’*C*R; Asemˇanˇator, dar cu o scriere mult mai complicatˇa se poate calcula matricea elementalˇa corespunzˇatoare unui element finit cu patru noduri, a se vedea [3]. Matricea de rigiditate A se calculeazˇa ˆın Matlab astfel: for j=1:size(elements3,1) I=2*elements3(j,[1,1,2,2,3,3]) -[1,0,1,0,1,0]; A(I,I)=A(I,I)+... melementala3(coordinates(elements3(j,:),:),lambda,mu); end for j=1:size(elements4,1) I=2*elements4(j,[1,1,2,2,3,3,4,4])-[1,0,1,0,1,0,1,0]; A(I,I)=A(I,I)+... melementala4(coordinates(elements4(j,:),:),lambda,mu); end Considerˇam ˆın aceastˇa sect¸iune urmˇatorul exemplu: o membranˇa elas√ ticˇa (E=2900, ν = 0.4) ocupˇa domeniul Ω = (−3, 3) × (−3, 3) \ B(0, 2). Dacˇa y = 3 sau y = −3 atunci f 2 = n. ˆIn rest, frontiera este liberˇa de tract¸iuni. Cum problema este simetricˇa, vom discretiza numai un sfert √ al domeniului Ω. Avem a¸sadar, u2 = 0 pe [ 2, 3] × {0} ¸si −u1 = 0 pe √ {0} × [ 2, 3].

Figura 5.4 prezintˇa me¸sa nedeformatˇa; Figura 5.5 reprezintˇa me¸sa deformatˇa cu un factor egal cu 200 iar prin codul culorilor, trecˆand de la albastru, care corespunde valorii minime, ¸si ajungˆand la ro¸su, care corespunde valorii maxime, este reprezentatˇa o medie nodalˇa a energiei elastice. ˆIn Sect¸iunea 6.1 este prezentat codul matlab complet pentru rezolvarea acestu exemplu. Acest cod a fost realizat de J. Alberty, C. Carstensen, S. A. Funken and R. Klose ¸si a fost publicat ˆın [3]. Putem testa ¸si de aceastˇa datˇa cˇa, diminuˆand f se diminueazˇa vectorul deplasare, astfel cˇa ˆın Figura 5.6 me¸sa este mai put¸in deformatˇa decˆat ˆın Figura 5.5.

5.3

Probleme 3D

ˆIn aceastˇa sect¸iune este pusˇa ˆın discut¸ie rezolvarea ˆın sens slab a problemelor ˆın deplasˇari-tract¸iuni, pentru materiale liniar elastice, omogene ¸si izotrope, admit¸and cˇa d = 3. Problema 5.29. Sˇa se determine u : Ω ⊂ R3 → R3 ¸si σ : Ω → S3 astfel ˆıncˆat Div σ + f 0 = 0

ˆın Ω,

σ = Cε(u)

ˆın Ω,

u = w

pe Γ1 ,

σ n = f2

pe Γ2 .

¯ → Prin solut¸ie tare a Problemei 5.29 ˆınt¸elegem o pereche (u, σ), u : Ω ¯ → S3 , u ∈ C 2 (Ω) ∩ C 1 (Ω), ¯ σ ∈ C 1 (Ω), care verificˇa (5.21)-(5.21). R3 , σ : Ω Admitem ˆın cele ce urmeazˇa cˇa densitatea fort¸elor volumice ¸si densitatea tract¸iunilor verificˇa ipoteza f 0 ∈ L2 (Ω)3 ,

f 2 ∈ L2 (Γ2 )3 .

(5.21)

Ca ¸si ˆın sect¸iunea precedentˇa admitem fˇarˇa a restrˆange generalitatea cˇa w = 0. Problema 5.29 poate fi formulatˇa variat¸ional dupˇa cum urmeazˇa. Problema 5.30. Dat f ∈ V, sˇa se determine u ∈ V astfel ˆıncˆat a(u, v) = (f , v)V

∀v ∈ V,

unde V este spat¸iul funct¸ional definit ˆın (1.13) pentru Ω ⊂ R3 . Definit¸ia 5.31. Problemei 5.30.

Numim solut¸ie slabˇa a Problemei 5.29 orice solut¸ie a

Reluˆand demersul realizat ˆın sect¸iunea precedentˇa, se poate demonstra urmˇatoarea teoremˇa. Teorema 5.32. Admitem (5.21). Atunci, Problema 5.30 are solut¸ie unicˇa. ˆ In plus, existˇa C > 0 astfel ˆıncˆat ∥u∥V ≤ C∥f ∥V

(5.22)

∥u1 − u2 ∥V ≤ C∥f 1 − f 2 ∥V

(5.23)

¸si unde u1 , u2 ∈ V sunt solut¸iile corespunzˇ atoare datelor f 1 , f 2 ∈ V. Fie Vh un subspat¸iu de dimensiune 3N al spat¸iului V. Considerˇam baza {η 1 , ..., η 3N } = {η1 e1 , η1 e2 , η1 e3 ..., ηN e1 , ηN e2 , ηN e3 } unde N este numˇarul nodurilor din ret¸eaua care discretizeazˇa domeniul. Considerˇam problema discretˇa ata¸satˇa Problemei 5.25. Problema 5.33. Dat f ∈ V, sˇa se determine uh ∈ Vh astfel ˆıncˆat a(uh , v h )V = (f , v h )V

∀uh ∈ Vh .

z

40

20

0 20 0 60

−20 40

−40 −60

20 −80

y

0 −100

x

Figura 5.7: Exemplul 3: discretizarea domeniului

Pe baza Teoremei 1.4 (Lax-Milgram), Problema 5.33 are solut¸ie unicˇa. ˆIn plus, calculul acestei unice solut¸ii se face rezolvˆand sistemul liniar Ax = b, ∫

unde A = (Ai,j );

ε(η k ) · Cε(η l )dx k, l ∈ 1, 3N

Akl = Ω

∫

¸si b = (bk );

∫ f 2 · η k dx +

bk = Ω

g · η k ds k ∈ 1, 3N . ΓN

ˆIn aceastˇa sect¸iune considerˇam un exemplu ˆın care domeniul Ω este domeniul ocupat ˆın R3 de o piesˇa cu o geometrie complexˇa, a se vedea Figura 5.7. Materialul este fixat ˆın douˇa gˇauri de centre (17, −1.5, 21) ¸si (48, −1.5, 21), astfel cˇa pe port¸iunile corespunzˇatoare de frontierˇa avem condit¸ii Dirichlet omogene. Port¸iunea de frontierˇa corespunzˇatoare unei

50 45

z

40 40

35

20

30 25

0 20

20 0 60

−20 40

−40 −60 y

15 10

20 −80

0 −100

5 x

Figura 5.8: Exemplul 3

gˇauri centrate ˆın punctul de coordonate (20, −77, 11.5) este supusˇa unei condit¸ii Dirichlet neomogene u_d=(0,0,0.1). Restul frontierei este liberˇa de tract¸iuni. ˆIn Figura 5.7 este prezentatˇa me¸sa nedeformatˇa. ˆIn Figura 5.8 este prezentatˇa me¸sa deformatˇa pentru 6512 grade de libertate iar deplasarea este amplificatˇa de 100 de ori. Prin intermediul codului culorilor este prezentatˇa mˇarimea energiei elastice. Pentru discretizarea domeniului s-au folosit 2172 noduri ¸si 8257 elemente finite (tetraedre). Codul matlab corespunzˇator acestui exemplu este prezentat in Sect¸iunea 6.1. El a fost realizat de J. Alberty, C. Carstensen, S. A. Funken and R. Klose ¸si a fost publicat ˆın [3]. Dificultatea principalˇa ˆın rezolvarea problemelor 3D constˇa ˆın realizarea geometriei domeniilor. Pentru detalii privind implementarea ˆın Matlab a problemelor elastos-

tatice recomandˇam lucrˇarile [2, 3, 4] ¸si referint¸ele acestora.

Capitolul 6 Anexˇ a

6.1

Coduri MATLAB

ˆIn scrierea acestei sect¸iuni am utilizat lucrˇarile [2, 3, 4].

Codul MATLAB pentru probleme antiplane main.m Initializari load coordinates.dat; coordinates(:,1)=[]; load elementes.dat; elementes(:,1)=[]; load dirichlet.dat; dirichlet(:,1) = []; load neumann.dat; neumann(:,1) = []; FreeNodes=setdiff(1:size(coordinates,1),unique(dirichlet)); A = sparse(size(coordinates,1),size(coordinates,1)); b=sparse(size(coordinates,1),1); Asamblarea matricei de rigiditate globala, din matrici elementale 117

for j = 1:size(elementes,1) A(elementes(j,:),elementes(j,:)) = A(elementes(j,:),elementes(j,:)) ... + melementala(coordinates(elementes(j,:),:)); end Calculul fortelor volumice for j=1:size(elementes,1) b(elementes(j,:))=b(elementes(j,:)) + ... det([1 1 1 ; coordinates(elementes(j,:),:)’])*... f(sum(coordinates(elementes(j,:),:))/3)/6; end Calculul tractiunilor for j = 1 : size(neumann,1) b(neumann(j,:)) = b(neumann(j,:)) + ... norm(coordinates(neumann(j,1),:)-... coordinates(neumann(j,2),:))*... g(sum(coordinates(neumann(j,:),:))/2)/2; end Conditiile Dirichlet u = sparse(size(coordinates,1),1); u(unique(dirichlet))=... u_d(coordinates(unique(dirichlet),:)); b = b-A*u; Calculul solutiei aproximative u(FreeNodes) = A(FreeNodes,FreeNodes)\b(FreeNodes) Reprezentarea grafica a solutiei aproximative subplot(1,2,1) show(elementes,coordinates,full(u)); title(’solutia aproximativa’) xlabel(’axa 0x’) ylabel(’axa 0y’) zlabel(’axa 0z’)

Reprezentarea grafica a solutiei exacte x=1:0.09:2; y=0:0.09:1; [X,Y]=meshgrid(x,y); U1=Y.^2; U2=(X-2).^2; U=U1.*U2; subplot(1,2,2) surf(X,Y,U) xlabel(’axa 0x’) ylabel(’axa 0y’) zlabel(’axa 0z’) title(’solutia exacta’) f_0.m function VolumeForce = f_0(v); VolumeForce=sparse(size(v,1),1); x=v(1); y=v(2); VolumeForce=VolumeForce-2*y^2-2*(x-2)^2; f_2.m function Stress = f_2(v); Stress=sparse(size(v,1),1); x=v(1); y=v(2); if x==1 Stress=2*y^2; else Stress=2*(x-2)^2; end show.m function show(elementes,coordinates,u) trisurf(elementes,coordinates(:,1),coordinates(:,2),...

u’,’facecolor’,’interp’) elementala.m function M = melementala(vertices) d=size(vertices,2); G=[ones(1,d+1);vertices’]\[zeros(1,d);eye(d)]; M =det([ones(1,d+1);vertices’])*G*G’/prod(1:d); end u_d.m function DirichletBoundaryValue = u_d(x) DirichletBoundaryValue = zeros(size(x,1),1);

Codul MATLAB pentru rezolvarea problemelor plane main.m Initializari E = 2900; nu = 0.4; mu = E/(2*(1+nu)); lambda = E*nu/((1+nu)*(1-2*nu)); load coordinates.dat; eval(’load elements3.dat;’,’elements3 = [];’); eval(’load elements4.dat;’,’elements4 = [];’); eval(’load neumann.dat;’,’neumann = [];’); load dirichlet.dat; Asamblarea matricei de rigiditate globala, din matrici elementale A = sparse(2*size(coordinates,1),2*size(coordinates,1)); b = zeros(2*size(coordinates,1),1); for j = 1:size(elements3,1) I = 2*elements3(j,[1,1,2,2,3,3]) -[1,0,1,0,1,0]; A(I,I) = A(I,I)+...

melementala3(coordinates(elements3(j,:),:),lambda,mu); end for j = 1:size(elements4,1) I = 2*elements4(j,[1,1,2,2,3,3,4,4]) -[1,0,1,0,1,0,1,0]; A(I,I)=A(I,I)+... melementala4(coordinates(elements4(j,:),:),lambda,mu); end Calculul fortelor volumice for j = 1:size(elements3,1) I = 2*elements3(j,[1,1,2,2,3,3]) -[1,0,1,0,1,0]; fs = f(sum(coordinates(elements3(j,:),:))/3)’; b(I)=b(I)+... det([1,1,1;coordinates(elements3(j,:),:)’])*[fs;fs;fs]/6; end for j = 1:size(elements4,1) I = 2*elements4(j,[1,1,2,2,3,3,4,4]) -[1,0,1,0,1,0,1,0]; fs = f(sum(coordinates(elements4(j,:),:))/4)’; b(I)=b(I)+... det([1,1,1;coordinates(elements4(j,1:3),:)’])*[fs;fs;fs;fs]/4; end Calculul tractiunilor if ~isempty(neumann) n = (coordinates(neumann(:,2),:)-... coordinates(neumann(:,1),:))*[0,-1;1,0]; for j = 1:size(neumann,1); I = 2*neumann(j,[1,1,2,2]) -[1,0,1,0]; gm=g(sum(coordinates(neumann(j,:),:))/2,n(j,:)/norm(n(j,:)))’; b(I) = b(I) +norm(n(j,:))*[gm;gm]/2; end end Conditiile Dirichlet DirichletNodes = unique(dirichlet);

[W,M] = u_d(coordinates(DirichletNodes,:)); B = sparse(size(W,1),2*size(coordinates,1)); for k = 0:1 for l = 0:1 B(1+l:2:size(M,1),2*DirichletNodes-1+k) =... diag(M(1+l:2:size(M,1),1+k)); end end mask = find(sum(abs(B)’)); A = [A, B(mask,:)’; B(mask,:), sparse(length(mask),length(mask))]; b = [b;W(mask,:)]; Calculul solutiei x = A\b; u = x(1:2*size(coordinates,1)); Reprezentarea grafica a solutiei aproximative [AvE,Eps3,Eps4,AvS,Sigma3,Sigma4] =... avmatrix(coordinates,elements3,elements4,u,lambda,mu); show(elements3,elements4,coordinates,AvS,u,lambda,mu); estimate =... aposteriori(coordinates,elements3,elements4,AvE,Eps3,Eps4,... AvS,Sigma3,Sigma4,u,lambda,mu) f.m function volforce = f(x); volforce = zeros(size(x,1),2); g.m function sforce = g(x,n) sforce = zeros(size(x,1),2); sforce(find(n(:,2)==1),2) = 1; show.m function show(elements3,elements4,...

coordinates,AvS,u,lambda,mu) for i=1:size(coordinates,1) AvC(i)=(mu/(24*(mu+lambda)^2)+1/(8*mu))*(AvS(i,1)+... AvS(i,4))^2+1/(2*mu)*(AvS(2)^2-AvS(1)*AvS(4)); end factor=20; colormap(1-gray) trisurf(elements3,factor*u(1:2:size(u,1))+... coordinates(:,1),... factor*u(2:2:size(u,1))+coordinates(:,2), ... zeros(size(coordinates,1),1), AvC, ’facecolor’,’interp’); hold on trisurf(elements4,factor*u(1:2:size(u,1))+... coordinates(:,1), ... factor*u(2:2:size(u,1))+coordinates(:,2), ... zeros(size(coordinates,1),1), AvC,’facecolor’,’interp’); view(0,90) hold off colorbar(’vert’) melementala3.m function M=melementala3(vertices,lambda,mu) PhiGrad = [1,1,1;vertices’]\[zeros(1,2);eye(2)]; R = zeros(3,6); R([1,3],[1,3,5]) = PhiGrad’; R([3,2],[2,4,6]) = PhiGrad’; C = mu*[2,0,0;0,2,0;0,0,1] +lambda*[1,1,0;1,1,0;0,0,0]; M = det([1,1,1;vertices’])/2*R’*C*R; melementala4 function M= melementala4(vertices,lambda,mu) R_11 = [2 -2 -1 1; -2 2 1 -1; -1 1 2 -2;

1 -1 -2 2]/6;

R_12 = [1 1 -1 -1; -1 -1 1 1; -1 -1 1 1; 1 1 -1-1]/4; R_22 = [2 1 -1 -2; 1 2 -2 -1; -1 -2 2 1; -2 -1 1 2]/6; F = inv([vertices(2,:)-vertices(1,:); vertices(4,:)-vertices(1,:)]); L = [lambda+2*mu,lambda,mu]; M= zeros(8,8); E = F’*[L(1),0;0,L(3)]*F; M(1:2:8,1:2:8) =E(1,1)*R_11 +... E(1,2)*R_12 +E(2,1)*R_12’ +E(2,2)*R_22; E =F’*[L(3),0;0,L(1)]*F; M(2:2:8,2:2:8) = E(1,1)*R_11 +... E(1,2)*R_12 +E(2,1)*R_12’ +E(2,2)*R_22; E = F’*[0,L(3);L(2),0]*F; M(2:2:8,1:2:8) = E(1,1)*R_11 +E(1,2)*R_12 +... E(2,1)*R_12’ +E(2,2)*R_22; M(1:2:8,2:2:8) = M(2:2:8,1:2:8)’; M = M/det(F); u_d.m function [W,M] = u_d(x) M = zeros(2*size(x,1),2); W = zeros(2*size(x,1),1); temp = find(x(:,1)>0 & x(:,2)==0); M(2*temp-1,2) = 1; temp = find(x(:,2)>0 & x(:,1)==0); M(2*temp-1,1) = 1; aposteriori.m function estimate=... aposteriori(coordinates,elements3,elements4, ... AvE,Eps3,Eps4,AvS,Sigma3,Sigma4,u,lambda,mu); eta3 = zeros(size(elements3,1),1);

eta4 = zeros(size(elements4,1),1); e3 = zeros(size(elements3,1),1); e4 = zeros(size(elements4,1),1); for j=1:size(elements3,1) Area3=det([1,1,1;coordinates(elements3(j,:),:)’])/2; for i=1:4 eta3(j) = eta3(j) + Area3 * (Sigma3(j,i)*Eps3(j,i) ... +AvS(elements3(j,:),i)’*[2,1,1;1,2,1;1,1,2]*... AvE(elements3(j,:),i)/12 ... - AvS(elements3(j,:),i)’*[1;1;1]*Eps3(j,i)/3 ... - AvE(elements3(j,:),i)’*[1;1;1]*Sigma3(j,i)/3); e3(j) = e3(j) + Area3 * ... AvS(elements3(j,:),i)’*[2,1,1;1,2,1;1,1,2]*... AvE(elements3(j,:),i)/12; end end for j=1:size(elements4,1) Area4=det([1,1,1;coordinates(elements4(j,1:3),:)’]); for i=1:4 eta4(j) = eta4(j) + Area4 * (Sigma4(j,i)*Eps4(j,i) ... +AvS(elements4(j,:),i)’*[4,2,1,2;2,4,2,1;1,2,4,2;2,1,2,4]*... AvE(elements4(j,:),i)/36 ... - AvS(elements4(j,:),i)’*[1;1;1;1]*Eps4(j,i)/4 ... - AvE(elements4(j,:),i)’*[1;1;1;1]*Sigma4(j,i)/4); e4(j) = e4(j) + Area4 * ... AvS(elements4(j,:),i)’*... [4,2,1,2;2,4,2,1;1,2,4,2;2,1,2,4]... AvE(elements4(j,:),i)/36; end end estimate = sqrt(sum(eta3)+sum(eta4))... /sqrt(sum(e3)+sum(e4));

avmatix.m function [AvE,Eps3,Eps4,AvS,Sigma3,Sigma4] = ... avmatrix(coordinates,elements3,elements4,u,... lambda,mu); Eps3 = zeros(size(elements3,1),4); Sigma3 = zeros(size(elements3,1),4); Eps4 = zeros(size(elements4,1),4); Sigma4 = zeros(size(elements4,1),4); AreaOmega = zeros(size(coordinates,1),1); AvS = zeros(size(coordinates,1),4); AvE = zeros(size(coordinates,1),4); for j = 1:size(elements3,1) area3 = det([1,1,1;coordinates(elements3(j,:),:)’])/2; AreaOmega(elements3(j,:))=... AreaOmega(elements3(j,:))+area3; PhiGrad=[1,1,1;coordinates(elements3(j,:),:)’]\... [zeros(1,2);eye(2)]; U_Grad =u([1;1]*2*elements3(j,:)-[1;0]*... [1,1,1])*PhiGrad; Eps3(j,:) = reshape((U_Grad+U_Grad’)/2,1,4); Sigma3(j,:) = reshape(lambda*trace(U_Grad)*... eye(2)+2*mu*(U_Grad+U_Grad’)/2,1,4); AvE(elements3(j,:),:) = AvE(elements3(j,:),:)... +area3*[1;1;1]*Eps3(j,:); AvS(elements3(j,:),:) = AvS(elements3(j,:),:)... +area3*[1;1;1]*Sigma3(j,:); end; for j = 1:size(elements4,1) area4 = det([1,1,1;coordinates(elements4(j,1:3),:)’]); AreaOmega(elements4(j,:),:)=... AreaOmega(elements4(j,:))+area4; Vertices = coordinates(elements4(j,:),:);

U_Grad=inv([Vertices(2,:)... -Vertices(1,:);Vertices(4,:)-... Vertices(1,:)])*[-1,1,1,-1;-1,-1,1,1]*... [u(2*elements4(j,:)-1),u(2*elements4(j,:))]/2; Eps4(j,:) = reshape((U_Grad+U_Grad’)/2,1,4); Sigma4(j,:) = reshape(lambda*... trace(U_Grad)*eye(2)+2*mu*(U_Grad+U_Grad’)/2,1,4); AvE(elements4(j,:),:) = AvE(elements4(j,:),:)... +area4*[1;1;1;1]*Eps4(j,:); AvS(elements4(j,:),:) = AvS(elements4(j,:),:)... +area4*[1;1;1;1]*Sigma4(j,:); end AvE = AvE./(AreaOmega*[1,1,1,1]); AvS =AvS./(AreaOmega*[1,1,1,1]); Variabilele Sigma3 ¸si Eps3 notezˇa tensorul tensiune ¸si tensorul deformat¸ie pentru elemente finite triunghiulare. Sigma4 ¸si Eps4 notezˇa tensorul tensiune ¸si tensorul deformat¸ie pentru elemente finite quadrilaterale. Variabilele AvS ¸si AvE sunt matrici cu 4 coloane, avˆand numˇarul liniilor egal cu numˇarul nodurilor din ret¸ea. Componentele fiecˇarei linii noteazˇa tensorul tensiune medie respectiv, tensorul deformat¸ie medie, pentru nodul notat cu numˇarul liniei curente. Discretizarea s-a fˇacut folosind 561 noduri ¸si 640 elemente finite (256 elemente finite cu trei noduri ¸si 384 elemente finite cu patru noduri). Pentru mai multe detalii a se consulta [3].

Codul MATLAB pentru rezolvarea problemelor 3D main.m Init¸ializˇari E=100000;nu=0.3; mu=E/(2*(1+nu)); lambda=E*nu/((1+nu)*(1-2*nu));

load coordinates.dat; load elements.dat; eval(’load neumann.dat;’,’neumann = [];’); load dirichlet.dat; A = sparse(3*size(coordinates,1),3*size(coordinates,1)); b = zeros(3*size(coordinates,1),1); Asamblarea matricei de rigiditate din matrici elementale for j = 1:size(elements,1) I=3*elements(j,[1,1,1,2,2,2,3,3,3,4,4,4])... -[2,1,0,2,1,0,2,1,0,2,1,0]; A(I,I) = A(I,I)+... stima(coordinates(elements(j,:),:),lambda,mu); end Calculul fortelor volumice for j = 1:size(elements,1) I=3*elements(j,[1,1,1,2,2,2,3,3,3,4,4,4])-... [2,1,0,2,1,0,2,1,0,2,1,0]; fs = f(sum(coordinates(elements(j,:),:))/4)’; b(I) = b(I)+... det([1,1,1,1;coordinates(elements(j,:),:)’])*... [fs;fs;fs;fs]/24; end Calculul tractiunilor if ~isempty(neumann) for j = 1:size(neumann,1) n = cross( coordinates(neumann(j,2),:)-... coordinates(neumann(j,1),:), ... coordinates(neumann(j,3),:)-coordinates(neumann(j,1),:)); I =3*neumann(j,[1,1,1,2,2,2,3,3,3])-[2,1,0,2,1,0,2,1,0]; gm = g(sum(coordinates(neumann(j,:),:))/3,n/norm(n))’; b(I) = b(I) +norm(n)*[gm;gm;gm]/6; end

end Conditii Dirichlet dirichletnodes = unique(dirichlet); [W,M] = u_d(coordinates(dirichletnodes,:)); B = sparse(size(W,1),3*size(coordinates,1)); for k = 0:2 for l = 0:2 B(1+l:3:size(M,1),3*dirichletnodes-2+k) = ... diag(M(1+l:3:size(M,1),1+k)); end end mask=find(sum(abs(B)’)); A=[A,B(mask,:)’;B(mask,:),sparse(length(mask),length(mask))]; b = [b;W(mask,:)]; Calculul solutiei aproximative x = A \ b; u = x(1:3*size(coordinates,1)); Reprezentare grafica a solutiei aproximative show(elements,dirichlet,neumann,coordinates,u,lambda,mu); f_0.m function volforce = f_0(x) volforce=zeros(size(x,1),3); f_2.m function sforce = f_2(x,n) sforce=zeros(size(x,1),3); show.m function show(elements,dirichlet,neumann,... coordinates,u,lambda,mu) E=zeros(4*size(elements,1),3);

C=zeros(size(elements,1),1); for j=1:size(elements,1) PhiGrad=[1,1,1,1;coordinates(elements(j,:),:)’]\... [zeros(1,3);eye(3)]; U_Grad=u([1;1;1]*3*elements(j,:)-... [2;1;0]*[1,1,1,1])*PhiGrad; SIGMA = lambda * trace(U_Grad)*eye(3)+... mu*(U_Grad +U_Grad’); C(j) = sqrt(sum(eig(SIGMA).^2)); end Area = zeros(size(elements,1),1); AreaOmega = zeros(max(max(elements)),1); AvC = zeros(max(max(elements)),1); for j=1:size(elements,1) Area = det([1,1,1,1;coordinates(elements(j,:),:)’])/6; AreaOmega(elements(j,:))=... AreaOmega(elements(j,:))+Area; AvC(elements(j,:)) =... AvC(elements(j,:))+Area*[1;1;1;1]*C(j); end AvC = AvC./AreaOmega; E=[dirichlet;neumann]; factor=100.0; colormap(1-gray); trisurf(E,factor*u(1:3:size(u,1))+coordinates(:,1), ... factor*u(2:3:size(u,1))+coordinates(:,2), ... factor*u(3:3:size(u,1))+coordinates(:,3),AvC,... ’facecolor’,’interp’); view(-50,35) axis equal; axis([-10 70 -100 20 0 40]) xlabel(’x’); ylabel(’y’);

zlabel(’z’); stima.m function stima = stima(vertices,lambda,mu) PhiGrad = [1,1,1,1;vertices’]\[zeros(1,3);eye(3)]; R = zeros(6,12); R([1,4,5],1:3:10) = PhiGrad’; R([4,2,6],2:3:11) = PhiGrad’; R([5,6,3],3:3:12) = PhiGrad’; C([1:3],[1:3]) = lambda*ones(3,3)+2*mu*eye(3); C([4:6],[4:6]) = mu*eye(3); stima = det([1,1,1,1;vertices’])/6*R’*C*R; u_d.m function [W,M] = u_d(x) M = zeros(3*size(x,1),3); W = zeros(3*size(x,1),1); M(1:3:end,1) = 1; M(2:3:end,2) = 1; M(3:3:end,3) = 1; W(3*find(x(:,2)