PROJEKT Z OPTYMALIZACJI KONSTRUKCJI

Autor:

Prowadzący: dr M. Rodak

Poznań 2012

SPIS TREŚCI

1. Wstęp

2. Model matematyczny konstrukcji

   1. Wykres sił poprzecznych

   2. Wykres momentów zginających

   3. Równanie linii ugięcia belki

   4. Charakterystyka geometryczna konstrukcji

3. Model optymalizacyjny

   1. Zmienne decyzyjne

   2. Ograniczenia

4. Rozwiązanie zadania

   1. Rozwiązanie numeryczne

   2. Rozwiązanie graficzne

5. Analiza rozwiązania

1. Wstęp

Celem projektu była optymalizacja belki podanej przez prowadzącego zajęcia. Wszystkie założenia i ograniczenia narzucone przez prowadzącego przedstawione zostaną w rozdziale 2. W wykonaniu projektu pomocny był program Matlab w którym została rozwiązana funkcja celu, oraz program Microsoft Excel, w którym wykonane zostało graficzne rozwiązanie zadania. Wykresy dotyczące belki zostały wykonane w programie Belka - wersja 1.0

Optymalizacja dotyczyła belki przedstawionej na rysunku 1:

Rysunek Optymalizowana belka

Przekrój poprzeczny belki przedstawiony został na rysunku 2:

Rysunek Przekrój poprzeczny belki

1. Model matematyczny konstrukcji

   1. Wykres sił poprzecznych

By przedstawić wykres sił poprzecznych należało wcześniej policzyć siły reakcyjne w podporach.

Suma sił wzdłuż osi y:

$\sum_{}^{}{F_{y} = \ R_{A} + \ R_{B} = 0}$

Suma momentów względem punktu A:

$$\sum_{}^{}{M_{A} = \ M - 2M - \ \ R_{B} = 0}$$

Z tego dostajemy:

R_A = − R_B oraz $R_{B} = \ - \frac{M}{L}$ = - 16,67 kN

Obrazowo jest to pokazane na rysunku 3

Rysunek Reakcje podporowe [kN]

Po wstawieniu myślowego przekroju o długości x od podpory A dostajemy:

Siła poprzeczna T(x) = R_A

Siła poprzeczna ma wartość stałą na całej długości belki. Jest to schemtycznie na przedstawione na rysunku 4

Rysunek Wykres sił poprzecznych [kN]

2.2 Wykres momentów zginających

Momenty zginające wyznaczone zostały przy pomocy tego samego przekroju myślowego.

M(x) = R_Ax + M

Wartość momentu w początku belki wynosi :

M(x) = 16, 67 * 0 + 20 = 20 kNm

W punkcie B belki x=L = 1,2m, czyli momenty wynosi:

M(x) = 16,67*1,2 + 20 = 40 kNm

Tak więc wartość momentu zmienia się od wartości 20kNm do 40kNm. Widać to na rysunku 5.

Rysunek Wykres momentów zginających [kN]

3. Równanie linii ugięcia belki

Równanie różniczkowe linii ugięcia belki wyraża się wzorem:

$$\text{EI}_{z}v"(x) = - M(x)$$

Moduł Younga dla stali z której wykonana jest belka wynosi 2,1 * 10⁵ MPa

Podstawiając do powyższego wzoru dostajemy:

$\backslash n\text{EI}_{z}v"(x) =$ −R_Ax - Mx⁰

Jest to postać potrzebna do obliczenia linii ugięcia metodą Clebscha

Całkujemy równanie

$\text{EI}_{z}v"(x) =$ −R_Ax -Mx⁰/ ∫dx

EI_zv′(x)= $- \frac{R_{A}x^{2}}{2}$ -Mx + C

Po kolejnym scałkowaniu

EI_zv′(x)= $- \frac{R_{A}x^{2}}{2}$ -Mx + C /∫dx

EI_zv(x) = − $\frac{R_{A}x^{3}}{6}\ $- $\frac{Mx^{2}}{2}$+ Cx + D

W punkcie A z powodu podpory v = 0 i x=0 więc D = 0

W punkcie B v = 0 i x = L więc 0 = − $\frac{R_{A}L^{3}}{6}$- $\frac{ML^{2}}{2}$+ Cl

$\frac{16,67*{1,2}^{2}}{6}$+ $\frac{20*1,2}{2}$+ =C

Z tego otrzymujemy C = 4,008 + 12 = 16,008 kNm²

By policzyć miejsce w którym ugięcie przyjmuje wartość maksymalną v_max należy równanie z v’ przyrównać do 0 i obliczyć gdzie jest ekstremum.

V’(x) = $\frac{1}{EI_{z}}$(-$R_{A}\frac{x^{2}}{2} - \text{Mx} + C\ $) = 0

Z warunku tego wynika, że -$R_{A}\frac{x^{2}}{2} - Mx + C\ $ = 0

-R_Ax² − 2Mx + 2C = 0 

-16,67 x² − 40x + 32 = 0

Δ = 1600 + 2132,48 = 3732,48

$\sqrt{\Delta}$ = 61,09

x = $\frac{40 \pm 61,09}{- 33,34}$

x = 0,63 v x= - 3, 03

Rozwiązanie -3,03 jest sprzeczne z założeniami bo x ≥ 0. W miejscu 0,63 m znajduje się więc wartość v_max

3. Charakterystyka geometryczna konstrukcji

Do charakterystyki geometrycznej konstrukcji zalicza się pole powierzchni przekroju poprzecznego A, oraz jego moment bezwładności I_z

Pole powierzchni przekroju poprzecznego jest sumą pól prostokątów składających się na dwuteownik.

A₁₌0, 8x₁ * 0, 2x₂ = 0, 16x₁x₂

A₂₌0, 2x₁ * 0, 6x₂ = 0, 12x₁x₂

A₃₌x₁ * 0, 2x₂ = 0, 20x₁x₂

A = A₁ + A₂ + A₃₌0, 48x₁x₂

Odległości środków ciężkości C₁,   C₂,   C₃ od osi x to y₁,   y₂,   y₃

y₁ = 0, 9x₂

y₂ = 0, 5x₂

y₃ = 0, 1x₂

Moment statyczny S_x =  S_x1 +  S_x2 +  S_x3   =  A₁y₁ +  A₂y₂ +  A₃y₃ =  0,16x₁x₂ * 0, 9x₂ + 0,12x₁x₂ * 0, 5x₂ +  0, 20x₁x₂ * 0, 1x₂ = 0, 144x₁x₂² + 0, 06x₁x₂² +  0, 02x₁x₂² =  0, 244x₁x₂²

Współrzędna środka ciężkości $y_{c} = \ \frac{S_{x}}{A}$ = $\frac{0,224x_{1}{x_{2}}^{2}}{0,48x_{1}x_{2}}$ = 0,46 x₂

Oś y jest osią symetrii

Dla elementów prostokątnych $I_{x0} = \ \frac{bh^{3}}{12}$

Korzystając z tego wzoru i z twierdzenia Steinera obliczony został moment bezwładności

I_z = 0, 062x₁x₂³ 

$W_{z} = \ \frac{I_{z}}{y_{\max}}$ = $\frac{0,062x_{1}{x_{2}}^{3}\ }{0,54x_{2}}$ = 0,115x₁x₂²

1. Model optymalizacyjny

   1. Zmienne decyzyjne

Zmiennymi decyzyjnymi w zadaniu są wartości x₁ i x₂. Całą konstrukcja optymalizowana jest ze względu na pole przekroju poprzecznego A. Wartość pola musi być jak najmniejsza.

1. Ograniczenia

Ograniczeniami przy optymalizowaniu belki są:

a) Dopuszczalne naprężenia σ_dop = 130 MPa

b) Dopuszczalne ugięcia belki f_max = 4mm

Wartości uzyskane dla belki muszą być mniejsze lub równe od dopuszczalnych.

c) $1 \leq \ \frac{x_{1}}{x_{2}} \leq 4$

d) x₁,  x₂ > 0

1. Rozwiązanie zadania

4.1 Rozwiązanie numeryczne

Zadanie rozwiązane zostało w programie Matlab. Poniżej znajdują się kody źródłowe z m-skryptów z programy wraz z komentarzem.

Skrypto

M=20*10^6 ; %Nmm

L=1.2*10^3; %mm

E=2.05*10^5; %MPa

sdop=130; %MPa

ydop=4; % mm

param=[M,L,E,sdop,ydop];

x0=[100,100]; % punkt początkowy poszukiwań

A=[1,-1;-4,1];

b=[0;0];

Aor=[];

bor=[];

lb=[0,0]; % x1>0 i x2>0

nb=[];

options=optimset('LargeScale','off');

[xm,fm]=fmincon(@fopt,x0,A,b,Aor,bor,lb,nb,@gopt,options,param)

gopt

function[onr,or]=gopt(x,param)

M=param(1);

L=param(2);

E=param(3);

sdop=param(4);

ydop=param(5);

Mmax=40*10^6; %Nmm

EIymax=5.42*10^9;

Wz=(0.115)*x(1)*x(2)^2;

Iz=(0.062)*x(1)*x(2)^3;

if Wz>=0

onr(1)=1-sdop*Wz/Mmax;

else

onr(1)=sdop*Wz/Mmax-1;

end

if Iz>=0

onr(2)=1-ydop*E*Iz/EIymax;

else

onr(2)=ydop*E*Iz/EIymax-1;

end

or=[];

fopt

function Aval=fopt(x,param)

Aval=0.48*x(1)*x(2);

Plik Skrypto zawiera dane charakteryzujące belkę oraz funkcję fmincon, która łączy ze sobą wszystkie ograniczenia i pozostałe dwie funkcje gopt i fopt. Jest to funkcji optymalizująca. Skrypt gopt zawiera warunki ograniczające natomiast plik fopt zawiera funkcję calu z uwagi na którą dokonywana jest optymalizacja.

Program matlab po wywołaniu pliku Skrypto wygenerował rozwiązanie. Spełnia ono podane ograniczenia. Poniżej znajduje się rozwiązanie z programu:

Warning: Options LargeScale = 'off' and Algorithm = 'trust-region-reflective' conflict.

Ignoring Algorithm and running active-set algorithm. To run trust-region-reflective, set

LargeScale = 'on'. To run active-set without this warning, use Algorithm = 'active-set'.

> In fmincon at 445

In skrypto at 15

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in

feasible directions, to within the default value of the function tolerance,

and constraints are satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

Active inequalities (to within options.TolCon = 1e-006):

lower upper ineqlin ineqnonlin

2 1

xm =

55.0934 220.3736

fm =

5.8277e+003

Wszystkie dane w m-skryptach zostały podane w takich jednostkach by wynik wyszedł w milimetrach.

x₁ = 55, 0934mm natomiast x₂ = 220, 3736mm

Pole powierzchni wynosi 5827,7mm²

4.2 Rozwiązanie graficzne

Rozwiązanie graficzne zostało wykonane w programie Microsoft Excel.

Linie proste zostały wyznaczone z warunków $1 \leq \ \frac{x_{1}}{x_{2}} \leq 4$

Hiperbola została wyznaczona z warunku wytrzymałościowego:

$$\frac{M_{\max}}{W_{z}} \leq \ \sigma_{\text{dop}}$$

Po dokonaniu przekształceń otrzymujemy:

$x_{2} \geq \ \sqrt{\frac{M_{\max}}{0,116x_{1}\sigma_{\text{dop}}}}$

Po wstawieniu odpowiednich danych i przeliczeniu jednostek otrzymujemy warunek :

$x_{2} \geq \ \sqrt{\frac{2,56*\ 10^{6}}{x_{1}}}$

Uzyskany wynik podany jest w [mm].

Kolejnym krokiem było wyznaczenie warstwic dla funkcji celu A = 0, 48x₁x₂

Po przekształceniach dostajemy:

$x_{2} = \ \frac{A}{0,48x_{1}}$

Z tego warunku podstawiając różne wartości pola powierzchni A dla zadanego przedziału wartości x₁ możemy uzyskać żądane warstwice.

Na kolejnych stronach przedstawione są wykresy uzyskane w programie Microsoft Excel. Pierwszy z nich przedstawia ogólny wygląd ograniczeń, obszaru dopuszczalnego i warstwic funkcji celu. Jednak ciężko z tego wykresu odczytać optymalną wartość, dlatego kolejny rysunek prezentuje okolice brzegów obszaru dopuszczalnego – miejsce w którym warstwice celu przecinają funkcji ograniczeń. W przypadku drugiego wykresu skupimy się na zakresie 35-250 dla x₂ i zakresie 40-72 dla x₁.

1. Analiza rozwiązania.

Na wykresie 2 widać, że Warstwica dla pola powierzchni 5500 mm² nie mieści się w zbierze rozwiązań dopuszczalnych. Warstwica 6000 spełnia warunki, ale daje gorszy wynik niż Warstwica 5800. Funkcji ta przecina przecina ograniczenie nr. 2 blisko miejsca w którym przecina je hiperbola dla warunku wytrzymałościowego. Jest to wynik zbliżony do tego uzyskanego w Matlabie. Dalsza optymalizacja mogłaby polegać na wyrysowaniu warstwic dla wartości bliskiej 5800 i sprawdzaniu czy istnieją lepsze rozwiązania niż odczytane z Wykresu 2 wartości x₁ = 54, 8 mm oraz  x₂ =  220, 48 mm. Metoda graficzna jest metodą przybliżoną, dlatego można zadowolić się uzyskanymi wynikami.