Řešení soustav rovnic (numericky), numerická integrace 6.1 Soustava lineárních rovnic – příklad ze statiky 6.1.1 Jednoduchá 2D úloha.
Rovnice statické rovnováhy
− F cos α + FAx − FBn cos β = 0 FAy + FBn sin β − F sin α = 0 − F cos α a + F sin α a + Μ − FB sin β a = 0 Řešení Ax = b F cos α ⎡1 0 − cos β ⎤ ⎡ FAx ⎤ ⎡ ⎤ ⎢0 1 ⎥ ⎢ ⎥ ⎢ ⎥ sin β ⎥ ⎢ FAy ⎥ = ⎢ F sin α ⎢ ⎥ ⎢⎣ 0 0 − sin β a ⎥⎦ ⎢⎣ FBn ⎥⎦ ⎢⎣ F cos α a − F sin α a − Μ ⎥⎦
Matlab % example 1 A = [1 0 -cosd(beta); 0 1 sind(beta); 0 0 -sind(beta)*a];
F = 10; % Newtons beta = 0; % degrees alpha = 45; % degrees M = 10; % Nm a = 0.5; % meters
% x = [Fax,Fay,Fb];
1/19
B = [F*cosd(alpha); F*sind(alpha); F*cosd(alpha)*a-F*sind(alpha)*a-M]; x = A\B
x= 41.7121 -12.9289 40.0000
2/19
Limitní případy β = 90 (a) a 0 (b) stupňů (a)
(b)
x= 7.0711 -12.9289 20.0000
x= NaN NaN Inf
Důvod – singulární matice A Změna výsledků při změně úhlu β % ruzne uhly beta counter = 1; for beta = 0:1:90 A = [1 0 -cosd(beta); 0 1 sind(beta); 0 0 -sind(beta)*a]; x = A\B; resultsFax(counter) = x(1); resultsFay(counter) = x(2); resultsFb(counter) = x(3); counter = counter +1; end plot(0:90,resultsFax,0:90,resultsFay,0:90,resultsFb); legend('Fax','Fay','Fb'); xlabel('beta angle [degrees]'); ylabel('Force [N]');
3/19
6.1.2 Jednoduchá 3D úloha.
Rovnice statické rovnováhy Fx : FAx + FD = − F sin α Fy : FAy − FB = F cos α Fz : FAz − FC = 0 a a − FC = − F cos α a 2 2 a a M y : − FC + FD = − F sin α a 2 2 a a M z : FB − FD = − F cos α a − F sin α a 2 2 M x : FB
Řešení Ax = b
⎡1 ⎢0 ⎢ ⎢0 ⎢ ⎢0 ⎢ ⎢ ⎢0 ⎢ ⎢ ⎢0 ⎣
0 0 0 1 0 −1 0 1 0 0 0
a 2
0 0
0
0 0
a 2
0 0 −1 a 2 a − 2 −
0
⎤ ⎥ − F sin α ⎤ ⎥ ⎡ FAx ⎤ ⎡ ⎢ ⎥ ⎢ ⎥ ⎥ F F cos α ⎥ ⎥ ⎢ Ay ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ F 0 ⎥ Az 0 ⎥ ⎥⎢F ⎥ = ⎢ − F cos α a ⎢ B⎥ ⎢ ⎥ ⎥ a ⎢ ⎥ ⎢ ⎥ ⎥ FC − F sin α a ⎥ 2 ⎥⎢ ⎥ ⎢ ⎢ F ⎥ − F cos α a − F sin α a ⎦ a⎥⎣ D ⎦ ⎣ − ⎥ 2⎦ 1 0 0
4/19
Matlab F = 10; % Newtons alpha = 45; % degrees a = 1; % meters A = [1 0 0 0 0 1; 0 1 0 -1 0 0; 0 0 1 0 -1 0; 0 0 0 a/2 -a/2 0; 0 0 0 0 -a/2 a/2; 0 0 0 a/2 0 -a/2]; % x = [Fax,Fay,Faz,Fb,Fc,Fd]; B = [-F*sind(alpha); F*cosd(alpha); 0; -F*cosd(alpha)*a; -F*sind(alpha)*a; -F*sind(alpha)*a - F*cosd(alpha)*a]; x = A\B % A je singulární matice % Pohněme vazbou B. Vzdálenost od roviny xy je označena c counter = 1; for c = 0:a/100:a A = [1 0 0 0 0 1; 0 1 0 -1 0 0; 0 0 1 0 -1 0; 0 0 0 c -a/2 0; 0 0 0 0 -a/2 a/2; 0 0 0 a/2 0 -a/2]; x = A\B; resultsFax(counter) = x(1); resultsFay(counter) = x(2); resultsFaz(counter) = x(3); resultsFb(counter) = x(4); resultsFc(counter) = x(5); resultsFd(counter) = x(6); counter = counter +1; end xaxisdata = 0:a/100:a; plot(xaxisdata,resultsFax,xaxisdata,resultsFay,xaxisdata,resultsFaz,... xaxisdata,resultsFb,xaxisdata,resultsFc,xaxisdata,resultsFd); legend('Fax','Fay','Faz','Fb','Fc','Fd'); xlabel('distance from xy plane [m]'); ylabel('Force [N]');
5/19
% Jakou roli hraje posunutí vazby B podél os x a y % Vzdálenosti jsou označeny c a d figure; clear resultsFax; clear resultsFay; clear resultsFaz; clear resultsFb; clear resultsFc; clear resultsFd;
counter1 = counter1 +1; end counter1 = 1; counter2 = counter2 + 1; end subplot(3,2,1); surf(resultsFax); title('Fax'); subplot(3,2,2); surf(resultsFay); title('Fay'); subplot(3,2,3); surf(resultsFaz); title('Faz'); subplot(3,2,4); surf(resultsFb); title('Fb'); subplot(3,2,5); surf(resultsFc); title('Fc'); subplot(3,2,6); surf(resultsFd); title('Fd');
counter1 = 1; counter2 = 1; for c = 0:a/15:a for d = 0:a/15:a A = [1 0 0 0 0 1; 0 1 0 -1 0 0; 0 0 1 0 -1 0; 0 0 0 c -a/2 0; 0 0 0 0 -a/2 a/2; 0 0 0 d 0 -a/2]; x = A\B; resultsFax(counter1,counter2) = x(1); resultsFay(counter1,counter2) = x(2); resultsFaz(counter1,counter2) = x(3); resultsFb(counter1,counter2) = x(4); resultsFc(counter1,counter2) = x(5); resultsFd(counter1,counter2) = x(6);
6/19
% Pohybujme vazbou B podel diagonály % Vysledky by mely byt singulární % Použijme velmi jemný krok figure; counter = 1; for c = 0:a/5000:a
7/19
A = [1 0 0 0 0 1; 0 1 0 -1 0 0; 0 0 1 0 -1 0; 0 0 0 c -a/2 0; 0 0 0 0 -a/2 a/2; 0 0 0 c 0 -a/2]; x = A\B; resultsFax(counter) = x(1); resultsFay(counter) = x(2); resultsFaz(counter) = x(3); resultsFb(counter) = x(4); resultsFc(counter) = x(5); resultsFd(counter) = x(6); counter = counter +1; end xaxisdata = 0:a/5000:a; plot(xaxisdata,resultsFax,xaxisdata,resultsFay,xaxisdata,resultsFaz,... xaxisdata,resultsFb,xaxisdata,resultsFc,xaxisdata,resultsFd); legend('Fax','Fay','Faz','Fb','Fc','Fd'); xlabel('distance from xy plane [m]'); ylabel('Force [N]'); % Jak je videt, Matlab obcas resení najde, ale je velmi velmi nepresne....
8/19
6.1.3 Soustava těles
Rovnice statické rovnováhy Fx 2 : FAx − FBx = 0 Fy 2 : FAy − FBy = 0 M z 2 : FBx a − FBy b − M 1 = 0 Fx 3 : FBx − FCx = 0 Fy 3 : FBy − FCy = 0 M z 3 : − FBy d + FBx c − M 2 = 0 Řešení Ax = b
⎡1 ⎢0 ⎢ ⎢0 ⎢ ⎢0 ⎢0 ⎢ ⎣0
0 −1 0 1 0 −1 0
a
−b
0
1
0
0
0
0
c
1 −d
0 ⎤ ⎡ FAx ⎤ ⎡ 0 ⎤ ⎢ ⎥ 0 ⎥⎥ ⎢ FAy ⎥ ⎢⎢ 0 ⎥⎥ 0 0 ⎥ ⎢ FBx ⎥ ⎢ M 1 ⎥ ⎥⎢ ⎥ = ⎢ ⎥ −1 0 ⎥ ⎢ FBy ⎥ ⎢ 0 ⎥ 0 −1⎥ ⎢ FCx ⎥ ⎢ 0 ⎥ ⎥⎢ ⎥ ⎢ ⎥ 0 0 ⎦ ⎣⎢ FCy ⎦⎥ ⎢⎣ M 2 ⎥⎦ 0 0
Matlab M1 = 2; %Nm M2 = 3; %Nm a = 2; % meters b = 3; c = 4; d = 0.5;
9/19
A = [1 0 -1 0 0 0; 0 1 0 -1 0 0; 0 0 a -b 0 0; 0 0 1 0 -1 0; 0 0 0 1 0 -1; 0 0 c -d 0 0] % x = [Fax,Fay,Fbx,Fby,Fcx,Fcy]; B = [0; 0; M1; 0; 0; M2] x = A\B %Řešení je nalezeno bez problémů. Nicméně úloha obsahuje skrytou pastičku. % Pokud a=b=c=d (dva čtverce) potom výsledná silová soustava je soustava sil protínajících % jednu přímku. % Tím nám odpadne jedna podmínka statické rovnováhy a chybí jedna rovnice, což vede na % singularitu matice A a neplatné řešení. b = a; c = a; d = a; A = [1 0 -1 0 0 0; 0 1 0 -1 0 0; 0 0 a -b 0 0; 0 0 1 0 -1 0; 0 0 0 1 0 -1; 0 0 c -d 0 0] x = A\B
10/19
6.2 Soustava nelineárních rovnic – příklad ze statiky
6.2.1 Rovnovážná poloha Zátěžná síla F je konstantní, stále kolmá na těleso. Úlohou je najít rovnovážnou polohu tělesa (příslušný úhel α ).
Fx : − F cos α + FAx = 0 Fy : F sin α + FAy − G = 0 M z : GR sin α − 2 RF = 0 Tato soustava rovnic nejde zapsat ve formě Ax = b , jelikož je nelineární (neznámý úhel α je argumentem goniometrické funkce). V ruce je řešení snadné: GR sin α − 2 RF = 0 2 RF sin α = GR 2F α = arcsin G Matlab obsahuje toolbox pro symbolickou matematiku, který využívá kernelu Maple.
Matlab % rovnovazna poloha – symbolicke reseni F = 1; R = 1; Fg = 4; syms Fax Fay alpha
11/19
% rovnice eqs = '-F*cos(alpha) + Fax, F*sin(alpha)+Fay-Fg, Fg*R*sin(alpha)-F*2*R'; result = solve(eqs,alpha,Fax,Fay) % result obsahuje symbolicke reseni % Fax = F*(-(-Fg^2+4*F^2)/Fg^2)^(1/2) % Fay = -(2*F^2-Fg^2)/Fg % alpha = asin(2*F/Fg) % Pro ziskani numerických hodnot provedeme substituci F,R a Fg Fax = subs(result.Fax) Fay = subs(result.Fay) alpha = subs(result.alpha) % Fax = 0.8660 % Fay = 3.5000 % alpha = 0.5236
Matlab obsahuje také numerický řešič nelineárních rovnic. Aby se dal použít, musí být nejprve definována funkce která určuje řešenou úlohu: ********* function Fun = e4fun(x); % x = [Fax,Fay,Alpha]; F = 1; % Newtons R = 1; % meters G = 4; % Newtons Fun = [-F*cos(x(3))+x(1); F*sin(x(3))+x(2)-G; G*R*sin(x(3))-2*R*F];
********** Řešič je použit následujícím způsobem: % rovnovazna poloha - numericky x0 = [0,0,0]; % odhad řešení – nemame tuseni, nastavime na nulu options = optimset('Display','iter'); % nastaveni options [x,fval] = fsolve(@e4fun,x0,options); % volani resice % vysledek je ulozen v poli x % x = 0.8660 3.5000 0.5236
12/19
6.2.2 – pasivní odpory Určete hodnotu síly F tak aby se těleso Z pohybovalo v určeném směru a smyslu. Uvažujte pasivní odpory ve vazbách A a C, součinitele jsou známé: fC , rA , f A
Rovnice statické rovnováhy
Případ 1 – je uvažována jako pasivní pouze vazba C, případ 2 – je uvažován i odporový moment M A . Fx 2 : FBx + FCT = 0 Fy 2 : FBy + FCn − F = 0 M z 2 : FCn a − F ( a + b ) = 0 Fx 3 : FAx − FCT = 0 Fy 3 : FAy − FCn − FZ | = 0 M z 3 : FCT R − FZ r + M A = 0 Přídavné rovnice FCT = f C FCn MA =
0 (case 1) f AcrC FAx2 + FAy2 (case 2)
Případ 1 – M A = 0 , soustava lineárních rovnic, Případ 2 – nelineární rovnice ! Případ 1: Ax = b
13/19
⎡0 ⎢0 ⎢ ⎢0 ⎢ ⎢1 ⎢0 ⎢ ⎢⎣ 0
0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0
⎤ ⎡ FAx ⎤ ⎡ 0 ⎤ ⎥ ⎢F ⎥ ⎢ 0 ⎥ ⎥ ⎢ Ay ⎥ ⎢ ⎥ a − ( a + b ) ⎥ ⎢ FBx ⎥ ⎢ 0 ⎥ ⎥⎢ ⎥ = ⎢ ⎥ 0 ⎥ ⎢ FBy ⎥ ⎢ 0 ⎥ − fC 0 ⎥ ⎢ FCn ⎥ ⎢ FZ ⎥ −1 ⎥⎢ ⎥ ⎢ ⎥ fC R 0 ⎥⎦ ⎣ F ⎦ ⎢⎣ FZ r ⎥⎦ fC 1
0 −1
Matlab Fz = 10; a = 0.3; b = 0.5; r = 0.4; R = 0.6; fc = 0.2;
% Newtons % meters % meters % meters % meters % Friction coefficient
A = [0 0 1 0 fc 0; 0 0 0 1 1 -1; 0 0 0 0 a -(a+b); 1 0 0 0 -fc 0; 0 1 0 0 -1 0; 0 0 0 0 fc*R 0]; % x = [Fax,Fay,Fbx,Fby,Fcn,F]; B = [0; 0; 0; 0; Fz; Fz*r]; x = A\B; % vysledek % x = 6.6667 43.3333 -6.6667 -20.8333 33.3333 12.5000 % je vazba F / Fz linearni ? for i = 1 : 50 Fz = i; B = [0; 0; 0; 0; Fz; Fz*r]; x = A\B; force(i) = x(6); end % for plot(force); xlabel('Load force'); ylabel('Force needed');
14/19
Případ 2 – nelineární rovnice, numerické řešení Definice funkce function Fun = e5fun(x); % x = [Fax,Fay,Fbx,Fby,Fcn,F]; Fz = 10; % Newtons a = 0.2; % meters b = 0.5; % meters r = 0.4; % meters R = 0.6; % meters fc = 0.2; % Friction coefficient, coupling C fa = 0.3; % Friction coefficient, coupling A c = 1.27; % pin wear-in coefficient (1.57 for new) ra = 0.02; % pin radius Fun = [x(3)+fc*x(5); x(4)+x(5)-x(6); a*x(5)-(a+b)*x(6); x(1)-fc*x(5); x(2)-x(5)-Fz; fc*R*x(5)-Fz*r+fa*c*ra*sqrt(x(1)^2+x(2)^2) ];
Volání řešiče x0 = [0,0,0,0,0,0]; % odhad řešení options = optimset('Display','iter'); [x,fval] = fsolve(@e5fun,x0,options); % result % x = 6.1437 40.7184 -6.1437 -21.9417 30.7184 8.7767
Jak ovlivní výslednou sílu součinitel smykového tření v ložisku ? Definice funkce musí být modifikována, aby mohla získat Fz a fa zvnějšku: function Fun = e5fun_2(x); % x = [Fax,Fay,Fbx,Fby,Fcn,F]; c = 1.27; % pin wear-in coefficient (1.57 for new) ra = 0.02; % pin radius a = 0.3; % meters b = 0.5; % meters r = 0.4; % meters R = 0.6; % meters fc = 0.2; % Friction coefficient global fa; % souc. smyk. treni, vazba A, definovana jako globalni global Fz; % Newtons, definovano jako globalni
15/19
Fun = [x(3)+fc*x(5); x(4)+x(5)-x(6); a*x(5)-(a+b)*x(6); x(1)-fc*x(5); x(2)-x(5)-Fz; fc*R*x(5)-Fz*r+fa*c*ra*sqrt(x(1)^2+x(2)^2) ];
Volání řešiče v jednoduchém cyklu: clear all; global Fz; global fa; fa = 0.00; for j = 1:5 for i = 1:10 Fz = i; x0 = [0,0,0,0,0,0]; % odhad řešení options = optimset('Display','off'); % potlačení výpisu iterací [x,fval] = fsolve(@e5fun_2,x0,options); results(j,i) = x(6); fzlist(i) = Fz; end falist(j) = fa; fa = fa + 0.1; end; plot(results'); legend('fa = 0.0','fa = 0.1','fa = 0.2','fa = 0.3','fa = 0.4'); xlabel('Load force (Fz) [N]'); ylabel('Resulting break force (F) [N]');
14
Resulting break force (F) [N]
12
10
8 fa = fa = fa = fa = fa =
6
4
0.0 0.1 0.2 0.3 0.4
2
0
1
2
3
4
5 6 Load force (Fz) [N]
7
8
16/19
9
10
6.3 Numerická integrace
Pro spojité zatížení najděte odpovídající výslednou sílu (velikost a polohu) 6.3.1 – jednoduchá funkce spojitého zatížení (půlkružnice) f ( x ) = R2 − ( x − R )
2
Nejprve je nutno definovat funkce které se budou integrovat (jak pro sílu tak moment). Určení polohy je z jednoduchého vztahu Fv* x* = M V % ********** Hlavni rutina ************* clear all; global R; R = 5; % vykresleni průběhu sily a momentu x = 0:0.05:2*R; Fx = e6fun_1(x); Mx = e6fun_2(x); subplot(2,1,1); [AX,H1,H2] = plotyy(x, Fx, x, Mx); set(get(AX(1),'Ylabel'),'String','Force'); set(get(AX(2),'Ylabel'),'String','Moment'); xlabel('x'); title('case 1'); % vypocet hodnot F a M numerickou integraci F = quad(@e6fun_1,0,2*R) M = quad(@e6fun_2,0,2*R) % vypocet pusobiste sily x_star = M / F % numericke vysledky % F = 39.2699, M = 196.3495, x_star = 5.0000 % DEFINICE FUNKCI ************************* function FunF = e6fun_1(x); global R; % parametr R je dodan zvenci % musí byt pouzit operator tect, aby byla funkce vektorizovana % (napriklad e6fun_1([0,1]) musí pracovat ok) FunF = sqrt(R^2-(x-R).^2); % *********** function FunM = e6fun_2(x);
17/19
global R; FunM = x.*sqrt(R^2-(x-R).^2);
6.3.2 – složitější funkce (není snadné odhadnout výsledek
( )
⎛1 2 ⎞ ⎜ ( x −3) sin x ⎟ ⎠
f ( x ) = 2 − e⎝ 3
% vykresleni průběhu sily a momentu x = 0:0.01:5; Fx = e6fun_3(x); Mx = e6fun_4(x); subplot(2,1,2); [AX,H1,H2] = plotyy(x, Fx, x, Mx); set(get(AX(1),'Ylabel'),'String','Force'); set(get(AX(2),'Ylabel'),'String','Moment'); xlabel('x'); title('case 2'); % vypocet hodnot F a M numerickou integraci F = quad(@e6fun_3,0,5) M = quad(@e6fun_4,0,5) % vypocet pusobiste sily x_star = M / F
% numericke vysledky % F = 5.3283, M = 12.7457, x_star = 2.3921 % DEFINICE FUNKCI ************************* function FunF = e6fun_3(x); FunF = 2-(exp(1/3.*(x-3) .* sin(x.^2))); % *********** function FunM = e6fun_4(x); FunM = x.*(2-(exp(1/3.*(x-3) .* sin(x.^2))));
18/19
19/19
