Modelowanie wspomagające projektowanie maszyn PROJEKT	Politechnika Poznańska Instytut Mechaniki Stosowanej Zakład Mechaniki Technicznej	2014/2015
Imię i nazwisko	Wydział	Kierunek
Kamil Jakubczak Bartosz Góra	BMiZ	MiBM
	Semestr/Stopień	Grupa
	2 / II	KMiU
Temat projektu	Nr projektu	Prowadzący
Wyznaczanie wykresów ugięć	54	dr inż. P. Fritzkowski
oraz częstości i postaci drgań własnych	Data oddania	Ocena
dla ramy płaskiej
Opis projektu Dla ramy statycznie niewyznaczalnej przedstawionej na rysunku wykonać: 1. W rozwiązaniu analitycznym: Uwolnić układ od więzów i zaznaczyć lokalne układy współrzędnych {x, y} Napisać równania statyki Napisać warunki brzegowe i warunki zgodności Wyznaczyć reakcje Wykonać wykresy momentów gnących Dobrać przekrój poprzeczny dla ramy na podstawie warunku wytrzymałościowego Wykonać wykresy ugięć dla całej ramy 2. W rozwiązaniu metodą sztywnych elementów skończonych: Wprowadzić podział na części Wyznaczyć sztywności dla wszystkich części Wprowadzić EST i SES Przyjąć stopnie swobody, współrzędne uogólnione i prędkości uogólnione Napisać energię kinetyczną i potencjalną Wyznaczyć macierze bezwładności i sztywności Wyznaczyć ugięcie ramy, wynikające z zadanego obciążenia Wykonać wykresy ugięć dla całej ramy Wyznaczyć częstości i postaci drgań własnych układu Wykonać wykresy postaci drgań własnych dla kilku pierwszych częstości drgań własnych 3. Dodatkowo przeprowadzić analizę wpływu podziału (dyskretyzacji) ramy na wyniki obliczeń

Rozwiązanie analityczne:

1.1. Uwolnienie ramy od więzów:

W rozpatrywanej belce są 4 przedziały nieciągłości

Podstawienie danych:

$$\text{dane} = \{ a \rightarrow 1.1,b \rightarrow 1.4,F \rightarrow 4000,q \rightarrow 2000,\text{EE} \rightarrow 2.1*10\hat{}11,\text{ro} \rightarrow 7850\};$$

1.2. Równania statyki:

HA − HD + F = =0,

$$\text{VA} + \text{VD} - q*a - \frac{3}{2}q*b = = 0,$$

$$\text{MA} - \frac{1}{2}q*a^{2} - F*\frac{1}{2}a - (a + \frac{1}{2}b)*\frac{3}{2}*q*b + \text{VD}*(a + b) + \text{HD}*a = = 0$$

Momenty gnące i linie ugięcia belek:

Przedział I:

$$Mg1\lbrack x1\_\rbrack: = - (\text{VA}*x1 - \text{MA} - \frac{1}{2}q*{x1}^{2})$$

$$\theta 1\lbrack x1\_\rbrack: = - \frac{1}{\text{EE}*\text{Jz}}*(\text{VA}*\frac{{x1}^{2}}{2} - \text{MA}*x1 - \frac{1}{6}q*{x1}^{3}) + C1$$

$$y1\lbrack x1\_\rbrack: = - \frac{1}{\text{EE}*\text{Jz}}*(\text{VA}*\frac{{x1}^{3}}{6} - \text{MA}*\frac{{x1}^{2}}{2} - \frac{1}{24}q*{x1}^{4}) + C1*x1 + D1$$

Przedział II:

$$Mg2\lbrack x2\_\rbrack: = (\text{HA}*x2 + q*\frac{1}{2}a^{2} - \text{VA}*a + \text{MA})$$

$$\theta 2\lbrack x2\_\rbrack: = \frac{1}{\text{EE}*\text{Jz}}*( - \text{HA}*\frac{{x2}^{2}}{2} - (q*\frac{1}{2}a^{2} + \text{VA}*a + \text{MA})*x2) + C2$$

$$y2\lbrack x2\_\rbrack: = \frac{1}{\text{EE}*\text{Jz}}*( - \text{HA}*\frac{{x2}^{3}}{6} - (q*\frac{1}{2}a^{2} + \text{VA}*a + \text{MA})*\frac{1}{2}{x2}^{2}) + C2*x2 + D2$$

Przedział III:

$$Mg3\lbrack x3\_\rbrack: = HD*x3 + VD*b - \frac{1}{2}*\frac{3}{2}*q*b^{2}$$

$$\theta 3\lbrack x3\_\rbrack: = - \frac{1}{EE*Jz}*(HD*\frac{{x3}^{2}}{2} + (VD*b - \frac{1}{2}*\frac{3}{2}*q*b^{2})*x3) + C3$$

$$y3\left\lbrack x3_{} \right\rbrack - \frac{1}{\text{EE}*\text{Jz}}*\left( \text{HD}*\frac{{x3}^{3}}{6} + \left( \text{VD}*b - \frac{1}{2}*\frac{3}{2}*q*b^{2} \right)*\frac{{x3}^{2}}{2} \right) + C3*x3 + D3$$

Przedział IV:

$$Mg4\lbrack x4\_\rbrack: = \text{VD}*x4 - \frac{1}{2}*\frac{3}{2}*q*{x4}^{2}$$

$$\theta 4\lbrack x4\_\rbrack: = - \frac{1}{\text{EE}*\text{Jz}}*(\text{VD}*\frac{{x4}^{2}}{2} - \frac{1}{2}*\frac{3}{2}*q*\frac{{x4}^{3}}{3}) + C4$$

$$y4\left\lbrack x4_{\_} \right\rbrack - \frac{1}{\text{EE}*\text{Jz}}*\left( \text{VD}*\frac{{x4}^{3}}{6} - \frac{1}{2}*\frac{3}{2}*q*\frac{{x4}^{4}}{12} \right) + C4*x4 + D4$$

3. Warunki brzegowe i warunki zgodności

y1[0] = =0,

θ1[0] = =0,

y4[0] = =0,

$$y2\left\lbrack \frac{a}{2} \right\rbrack = = - y3\left\lbrack \frac{a}{2} \right\rbrack,$$

$$\theta 2\left\lbrack \frac{a}{2} \right\rbrack = = \theta 3\left\lbrack \frac{a}{2} \right\rbrack,$$

θ1[a] = =θ2[0],

θ3[0] = =θ4[b],

y4[b] = = − y1[a],

y3[0] = =0,

y2[0]= = 0}

1.4. Wyznaczenie reakcji i stałych całkowania:

F + HA − HD = =0,

$$- \text{aq} - \frac{3\text{bq}}{2} + VA + VD = = 0,$$

$$- \frac{\text{aF}}{2} + aHD + MA - \frac{a^{2}q}{2} - \frac{3}{2}\left( a + \frac{b}{2} \right)\text{bq} + \left( a + b \right)VD = = 0,$$

D1 = =0, C1 = =0, D4 = =0,

$$\frac{aC2}{2} + D2 + \frac{- \frac{a^{3}\text{HA}}{48} - \frac{1}{8}a^{2}\left( MA + \frac{a^{2}q}{2} + a\text{VA} \right)}{\text{EEJz}} = = - \frac{aC3}{2} - D3 + \frac{\frac{a^{3}\text{HD}}{48} + \frac{1}{8}a^{2}\left( - \frac{3b^{2}q}{4} + b\text{VD} \right)}{\text{EEJz}},$$

$$C2 + \frac{- \frac{a^{2}\text{HA}}{8} - \frac{1}{2}a\left( MA + \frac{a^{2}q}{2} + a\text{VA} \right)}{\text{EEJz}} = = C3 - \frac{\frac{a^{2}\text{HD}}{8} + \frac{1}{2}a\left( - \frac{3b^{2}q}{4} + b\text{VD} \right)}{\text{EEJz}},$$

$$C1 - \frac{- aMA - \frac{a^{3}q}{6} + \frac{a^{2}\text{VA}}{2}}{\text{EEJz}} = = C2,$$

$$C3 = = C4 - \frac{- \frac{b^{3}q}{4} + \frac{b^{2}\text{VD}}{2}}{\text{EEJz}},$$

$$bC4 + D4 - \frac{- \frac{b^{4}q}{16} + \frac{b^{3}\text{VD}}{6}}{\text{EEJz}} = = - aC1 - D1 + \frac{- \frac{a^{2}\text{MA}}{2} - \frac{a^{4}q}{24} + \frac{a^{3}\text{VA}}{6}}{\text{EEJz}},$$

D3 = =0, D2 = =0

Rozwiązanie:

HA → −3788.72, VA → 5246.91, MA → 7854.88, HD → 211.28, VD → 1153.09, C1 → 0,

$$\mathbf{C2} \rightarrow \frac{2.81412 \times 10^{- 8}}{\text{Jz}},\mathbf{C3} \rightarrow - \frac{1.13073 \times 10^{- 8}}{\text{Jz}},\mathbf{C4} \rightarrow - \frac{1.24596 \times 10^{- 8}}{\text{Jz}},\mathbf{D1} \rightarrow 0,\mathbf{D2} \rightarrow 0,\mathbf{D3} \rightarrow 0,\mathbf{D4} \rightarrow 0$$

5. Wykresy momentów gnących:

- Wykres momentów dla przedziału I

- Wykres momentów dla przedziału I

- Wykres momentów dla przedziału I

- Wykres momentów dla przedziału I

5. Wykresy ugięć dla całej ramy :

1. Rozwiązanie metodą Sztywnych Elementów Skończonych

   1. Podział na części:

1. Współczynniki sztywności każdego z elementów:

CN1=(EE*Ap)/ΔL1 – współczynnik wytrzymałości pierwszych elementów na rozciąganie

CT1=(GG*Ap)/(kt*ΔL1) – współczynnik wytrzymałości pierwszych elementów na ścinanie

CG1=(EE*Jz)/ ΔL1 – współczynnik wytrzymałości pierwszych elementów na zginanie

CN2=(EE*Ap)/ ΔL2 – współczynnik wytrzymałości drugich elementów na rozciąganie

CT2=(GG*Ap)/(kt*ΔL2) – współczynnik wytrzymałości drugich elementów na ścinanie

CG2=(EE*Jz)/ ΔL2 – współczynnik wytrzymałości drugich elementów na zginanie

CN3=(EE*Ap)/ ΔL3 – współczynnik wytrzymałości trzecich elementów na rozciąganie

CT3=(GG*Ap)/(kt*ΔL3) – współczynnik wytrzymałości trzecich elementów na ścinanie

CG3=(EE*Jz)/ ΔL3 – współczynnik wytrzymałości trzecich elementów na zginanie

CN4=(EE*Ap)/ ΔL4 – współczynnik wytrzymałości czwartych elementów na rozciąganie

CT4=(GG*Ap)/(kt*ΔL4) – współczynnik wytrzymałości czwartych elementów na ścinanie

CG4=(EE*Jz)/ ΔL4 – współczynnik wytrzymałości czwartych elementów na zginanie

CN1u=5*CN1; – współczynnik wytrzymałości utwierdzenia w początku ramy na rozciąganie

CT1u=5*CT1; – współczynnik wytrzymałości utwierdzenia w początku ramy na ścinanie

CG1u=5*CG1; – współczynnik wytrzymałości utwierdzenia w początku ramy na zginanie

CN4p=5*CN4; – współczynnik wytrzymałości podpory na końcu ramy na rozciąganie

CT4p=5*CT4; – współczynnik wytrzymałości podpory na końcu ramy na ścinanie

CG4p=0; – współczynnik wytrzymałości podpory na końcu ramy na zginanie

1. Wprowadzenie EST i SES:

W narożnikach nie powinna być L ??

1. Liczba stopni swobody:

Nss = 3 * Ne = 3 * (N1+N2+N3+N4+1) = 3 * (4+5+7+5+1) = 3 * 22 = 66 

Nss = 66

1. Kod generujący wektor współrzędnych uogólnionych:

Do[

{wq[3 * k − 2]=qx[k],

wq[3*k−1] = qy[k],

wq[3*k] = φ[k]},

{k, 1, Ne}]

Postać wektora współrzędnych uogólnionych – przemieszczenia w osiach X i Y oraz obrotów w płaszczyźnie XY:

{qx[1],qy[1],φ[1],qx[2],qy[2],φ[2],qx[3],qy[3],φ[3],qx[4],qy[4],φ[4],qx[5],qy[5],φ[5],

qx[6], qy[6], φ[6], qx[7], qy[7], φ[7], qx[8], qy[8], φ[8], qx[9], qy[9], φ[9], qx[10], qy[10], φ[10],

qx[11], qy[11], φ[11], qx[12], qy[12], φ[12], qx[13], qy[13], φ[13], qx[14], qy[14], φ[14],

qx[15], qy[15], φ[15], qx[16], qy[16], φ[16], qx[17], qy[17], φ[17], qx[18], qy[18], φ[18],

qx[19],qy[19],φ[19],qx[20],qy[20],φ[20],qx[21],qy[21],φ[21],qx[22],qy[22],φ[22]}

1. Kod generujący wektor sił uogólnionych:

Do[

 {wp[3 * k − 2]=fx[k],

wp[3*k−1] = fy[k],

wp[3 * k]=mφ[k]}, 

{k, 1, Ne}]

Postać wektora sił uogólnionych dla podanego wcześniej podziału ramy –
– kolejne wiersze odpowiadają wierszom wektora współrzędnych uogólnionych:

$$\{ 0, - \frac{\text{aq}}{8},0,0, - \frac{\text{aq}}{4},0,0, - \frac{\text{aq}}{4},0,0, - \frac{\text{aq}}{4},0,0, - \frac{\text{aq}}{8},0,$$

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, −F, 0,

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

$$0,0,0,0,0,0,0, - \frac{\text{bq}}{10},0,0, - \frac{\text{bq}}{5},0,$$

$$0, - \frac{\text{bq}}{5},0,0, - \frac{\text{bq}}{5},0,0, - \frac{\text{bq}}{5},0,0, - \frac{\text{bq}}{5},0\}$$

1. Zapis energii potencjalnej sprężystości:

$$\text{Ep} = \frac{1}{2}*(\sum_{k = 1}^{N1}\left( CN1*\left( \text{xA}\left\lbrack k + 1 \right\rbrack - \text{xB}\left\lbrack k \right\rbrack \right)^{2} + CT1*\left( \text{yA}\left\lbrack k + 1 \right\rbrack - \text{yB}\left\lbrack k \right\rbrack \right)^{2} + CG1*\left( \varphi\left\lbrack k + 1 \right\rbrack - \varphi\left\lbrack k \right\rbrack \right)^{2} \right)$$

$$+ \sum_{k = N1 + 2}^{N1 + N2}\left( CN2*\left( \text{xA}\left\lbrack k + 1 \right\rbrack - \text{xB}\left\lbrack k \right\rbrack \right)^{2} + CT2*\left( \text{yA}\left\lbrack k + 1 \right\rbrack - \text{yB}\left\lbrack k \right\rbrack \right)^{2} + CG2*\left( \varphi\left\lbrack k + 1 \right\rbrack - \varphi\left\lbrack k \right\rbrack \right)^{2} \right)$$

$$+ \sum_{k = N1 + N2 + 1}^{N1 + N2 + N3 - 1}\left( CN3*\left( \text{xA}\left\lbrack k + 1 \right\rbrack - \text{xB}\left\lbrack k \right\rbrack \right)^{2} + CT3*\left( \text{yA}\left\lbrack k + 1 \right\rbrack - \text{yB}\left\lbrack k \right\rbrack \right)^{2} + CG3*\left( \varphi\left\lbrack k + 1 \right\rbrack - \varphi\left\lbrack k \right\rbrack \right)^{2} \right)$$

$$+ \sum_{k = N1 + N2 + N3 + 1}^{N1 + N2 + N3 + N4}\left( CN4*\left( \text{xA}\left\lbrack k + 1 \right\rbrack - \text{xB}\left\lbrack k \right\rbrack \right)^{2} + CT4*\left( \text{yA}\left\lbrack k + 1 \right\rbrack - \text{yB}\left\lbrack k \right\rbrack \right)^{2} + CG4*\left( \varphi\left\lbrack k + 1 \right\rbrack - \varphi\left\lbrack k \right\rbrack \right)^{2} \right)$$

+(CN2*(xA[N1+2]−yB[N1+1])²+CT2*(yA[N1+2]+xB[N1+1])²+CG2*(φ[N1+2]−φ[N1+1])²)

+(CN4 * (yA[N1+N2+N3+1]−xB[N1+N2+N3])² + CT4 * (xA[N1+N2+N3+1]+yB[N1+N2+N3])²

+CG4 * (φ[N1+N2+N3+1]−φ[N1+N2+N3])²)+(CN1u * xA[1]² + CT1u * yA[1]² + CG1u * φ[1]²)

+(CN4p * xB[Ne]² + CT4p * yB[Ne]² + CG4p * φ[Ne]²))

Obliczenie pochodnej energii potencjalnej po przemieszczeniach uogólnionych:

pEp = ∂_{WQ}Ep

1. Kod generowania macierzy sztywności:

Do[

{Do[

{mC[kw,kk]=Coefficient[pEp[[kw]],WQ[[kk]]]},

{kk, 1, Nss}]},

{kw, 1, Nss}]

Zapis wektora bezwładności:

k=1

wA[3*k-2]=m1k

wA[3*k-1]=m1k

wA[3*k]=J1k

Do[

{wA[3*k-2]=m1

wA[3*k-1]=m1

wA[3*k]=J1},

{k,2,N1}]

k=N1+1

wA[3*k-2]=m12

wA[3*k-1]=m12

wA[3*k]=J12

Do[

{wA[3*k-2]=m2

wA[3*k-1]=m2

wA[3*k]=J2},

{k,N1+2,N1+N2}]

k=N1+N2+1

wA[3*k-2]=m23

wA[3*k-1]=m23

wA[3*k]=J23

Do[

{wA[3*k-2]=m3

wA[3*k-1]=m3

wA[3*k]=J3},

{k,N1+N2+2,N1+N2+N3}]

k=N1+N2+N3+1

wA[3*k-2]=m34

wA[3*k-1]=m34

wA[3*k]=J34

Do[

{wA[3*k-2]=m4

wA[3*k-1]=m4

wA[3*k]=J4},

{k,N1+N2+N3+2,Ne-1}]

k=Ne

wA[3*k-2]=m4k

wA[3*k-1]=m4k

wA[3*k]=J4k

1. Polecenie rozwiązania liniowego równania ugięcia statycznego ramy – użycie macierzy sztywności MC i wektora sił uogólnionych WP z podstawionymi danymi:

WQs = LinearSolve[numMC, numWP];

Ugięcie statyczne poszczególnych SES:

{0.00000148403324748036, −0.000060122628366477284, −0.0008565217098124666,

0.000008904199484882162, −0.0007443514499637101, −0.004505465450492198,

0.00001632436572228396, −0.0023308660882572965, −0.006995319150580114,

0.000023744531959685762, −0.004461874836825017, −0.008470402244913574,

0.00010246853809133448, −0.006432510702039443, −0.009075034168358694,

−0.006878597162743456, −0.00102189490536133, −0.008768253062643413,

−0.006879034624108425, −0.0019355249074032763, −0.007673212257044223,

−0.006879472085473393, −0.002685346137335691, −0.005789911751561016,

−0.006879909546838362, −0.003184650028171333, −0.0031183515461938134,

−0.006880347008203331, −0.0033494109786012583,  0.00034146835905766537,

−0.006880659480606876, −0.003216498918423674,  0.0027812884880855896,

−0.006880971953010422, −0.0029124305145222253,  0.004766468962022089,

−0.006881284425413968, −0.0024702444884043538,  0.006297009780867255,

−0.006881596897817513, −0.0019256625272558075,  0.007372910944621075,

−0.006881909370221059, −0.0013144063182623342,  0.007994172453283585,

−0.006882221842624605, −0.0006721975486096804,  0.008160794306854765,

0.00007511957929970266, −0.006493359215955717,  0.007134403595983036,

0.00003587090421465957, −0.005096517399065743,  0.005651101508422104,

0.000027330212734978726, −0.003640252414504495, 0.004768240592646999,

0.000018789521255297876, −0.0023677438546065603, 0.004333485450473267,

0.000010248829775617025, −0.0011748490585093572, 0.004194500683716596,

0.000001708138295936171, −0.00029400583943541085,  0.004198950894192523}

1. Wykres ugięcia ramy pod wpływem obciążenia statycznego w odniesieniu to belki nieugiętej:

1. Wyznaczenie postaci i kwadratów częstości drgań własnych ramy przy pomocy funkcji Eigensystem:

rozw = Eigensystem[numZMC]

Wartości drgań własnych w odwrotnej kolejności:

{128459.18293305453,  119202.61081023372,  105079.10981338337,  99778.12645000397,  99391.40282294097,

97328.04731386401,  96702.03051949179,  93374.0458405322,  91722.37877143902,  87310.95069779916,

83722.07151046033,  80663.7029072566,  79042.17596664668,  76587.73220810358,  74128.14261508365,

69845.07312737602,  68648.97529924236,  68424.91742872453,  65688.1663141324,  65266.39677376112,

64515.16271236334,  62836.81532963897,  60296.3683910468,  58303.85869033581,  53790.754782343254,

48636.387195237556,  46539.593385551474,  41578.5648959315,  38364.133095476325,  36997.09802593814,

35953.04341842696,  35577.0154728696,  34946.7459525416,  34466.010521689524,  33370.20944117213,

33301.51355593985,  31914.39719985126,  31164.420104634926,  30344.48796354008,  29329.0673828962,

27698.22373810031,  26784.885777212756,  25813.882003101302,  24740.881335167433,  24114.747840657063,

22841.329016462565,  20878.23263143867,  18448.95854001643,  16485.478943752114,  15202.689346301562,

12874.78662017514,  11416.718461893699,  10527.789610436288,  8946.71672139489,  8110.1543836407,

7542.731114798285,  6482.989725796055,  5365.170070958722,  4848.472919877072,  4780.274628203378,

3591.131772099266,  2378.4501119457623,  1765.7976314594168,  1125.564257793941,  667.8768032124542,

163.40958188126916}

1. Wykresy 3 pierwszych postaci drgań własnych ramy w odniesieniu do nieruchomej belki:
   - pierwsza postać drgań własnych
   - druga postać drgań własnych
   - trzecia postać drgań własnych

Analiza wpływu podziału (dyskretyzacji) ramy na wyniki obliczeń

Wykres ugięcia ramy wyznaczony metodą analityczną:

Wykres ugięcia ramy wyznaczony metodą sztywnych elementów skończonych:

Na podstawie powyższych wykresów można stwierdzić, że metoda sztywnych elementów skończonych (MSES) jest metodą stosunkowo dokładną - otrzymane wykresy są zbliżone do siebie. Ponieważ jest to metoda dyskretna i wyniki przemieszczeń otrzymujemy tylko w miejscach zamocowania elementów sprężystych, jej dokładność w dużym stopniu zależy od przyjętego podziału ramy.

Dla zwiększonego podziału ramy otrzymujemy następujące postaci drgań dla 3 pierwszych częstości:

N1 = 8, N2 = 11, N3 = 13, N4 = 9

- pierwsza postać drgań własnych
- druga postać drgań własnych
- trzecia postać drgań własnych

Na pierwszy rzut oka można dostrzec przeciwne wartości drgań własnych ramy drugiej i trzeciej postaci.

Kolejne zwiększenie stopnia dyskretyzacji ramy skutkuje poniższymi postaciami drgań:

N1 = 15, N2 = 21, N3 = 17, N4 = 14

- pierwsza postać drgań własnych
- druga postać drgań własnych
- trzecia postać drgań własnych

Tutaj z kolei, w porównaniu z najmniejszą liczbą elementów podziału, przeciwny znak ma druga postać drgań.

Z powyższej analizy wynika wniosek, że postać drgań bardzo mocno zależy od dyskretyzacji modelowanego układu i w zasadzie niezmiennikiem są wartości przemieszczeń dla poszczególnych wartości częstości drgań własnych, znak zachowuje się zależnie od liczby (a nawet jej parzystości) elementów, na jakie zostały podzielone elementy układu.