[ «MATLAB В ИНЖЕНЕРНЫХ И НАУЧНЫХ РАСЧЕТАХ Одесса Астропринт 2003 ББК Д УДК 539.3:681.3 Монография посвящена иллюстрации возможностей одной из самых эффективных систем компьютерной математики MATLAB в решении ряда научных и ...»
a(11,12)=1; a(11,13)= 1/2; a(11,14)= 1/6; a(11,16)=1; a(12,12)=1; a(12,13)=1;
a(12,14)= 1/2; a(13,1) = 1; a(13,13)=1; a(13,14)=1; a(14,3) = 1; ; a(14,14)=1;
a(15,5)= 1; a(15,15)= 1; b(11,1) =5/4; b(12,1) =15/2; b(13,1)=30; b(14,1)= 60;
a(16,16)=1; a(16,17)=1; a(16,18)=1/2; a(16,19)=1/6; a(17,17)=1; a(17,18)= 1;
a(17,19)=1/2; a(18,6)= 1; a(18,18)=1; a(18,19)=1; a(19,7)= 1; a(19,19)=1;
a(20,11)= 1; a(20,20)=1; ; b(16,1) =45/8; b(17,1)= 15; ; b(18,1)= 20;
X=a\b Значения граничных параметров рамы сведены в таблицу 3.6.
3.3.2. Построение эпюр напряженно-деформированного состояния Для построения эпюр программируем выражения (3.12), (3.13) х = 0 : 0.001 : 1.0; q = 10.0;
EIv = (X(2,1)*x X(4,1)*x.^ 3 / 6 + q*x.^4 / 24);
EIfi = (X(2,1) X(4,1)*x.^ 2 / 2 + q*x.^3 / 6);
Q = X(4,1) q*x; M = X(4,1)*x q*x.^ 2 / 2;
subplot (2,2,1), plot (x, EIv); axis ([0 1 0.4 0.4]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 1 2 2]); grid on subplot (2,2,3), plot (x, Q); axis ([0 1 10 10]); grid on subplot (2,2,4), plot (x, M); axis ([0 1 5 5 ]); grid on х = 0 : 0.001 : 1.0; q = 10.0;
EIv = (X(17,1)*x X(8,1)*x.^ 2 / 2 X(9,1)*x.^ 3 / 6 + q*x.^4 / 24);
EIfi = (X(17,1) X(8,1)*x X(9,1)*x.^ 2 / 2 + q*x.^3 / 6);
Q = X(9,1) q*x; M = X(8,1) + X(9,1)*x q*x.^ 2 / 2;
subplot (2,2,1), plot (x, EIv); axis ([0 1 0.4 0.4]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 1 2 2]); grid on subplot (2,2,3), plot (x, Q); axis ([0 1 22 22]); grid on subplot (2,2,4), plot (x, M); axis ([0 1 10 10]); grid on х = 0 : 0.001 : 0.75;
EIv = (X(16,1) + X(17,1)*x X(18,1)*x.^ 2 / 2 X(19,1)*x.^ 3 / 6);
EIfi = (X(17,1) X(18,1)*x X(19,1)*x.^ 2 / 2);
Q = X(19,1); M = X(18,1) + X(19,1)*x;
subplot (2,2,1), plot (x, EIv); axis ([0 0.25 3 3]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 0.25 6 6]); grid on subplot (2,2,3), plot (x, Q); axis ([0 0.25 2 2]); grid on subplot (2,2,4), plot (x, M); axis ([0 0.25 14 14]); grid on х = 0.25 : 0.001 : 1.0; m = 20.0;
EIv=(X(16,1)+X(17,1)*xX(18,1)*x.^ 2/2X(19,1)*x.^ 3/6+m*(x0.25).^ 2/2);
EIfi = (X(17,1) X(18,1)*x X(19,1)*x.^ 2 / 2+m*(x0.25));
Q = X(19,1); M = X(18,1) + X(19,1)*x m;
subplot (2,2,1), plot (x, EIv); axis ([0.25 1.0 3 3]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0.25 1.0 6 6]); grid on subplot (2,2,3), plot (x, Q); axis ([0.25 1.0 2 2]); grid on subplot (2,2,4), plot (x, M); axis ([0.25 1.0 15 15]); grid on х = 0 : 0.001 : 0.5;
EIv = (X(16,1) + X(12,1)*x X(13,1)*x.^ 2 / 2 X(14,1)*x.^ 3 / 6);
EIfi = (X(12,1) X(13,1)*x X(14,1)*x.^ 2 / 2);
Q = X(14,1); M = X(13,1) + X(14,1)*x;
subplot (2,2,1), plot (x, EIv); axis ([0 0.5 3 3]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 0.5 6 6]); grid on subplot (2,2,3), plot (x, Q); axis ([0 0.5 2 2]); grid on subplot (2,2,4), plot (x, M); axis ([0 0.5 15 15]); grid on х = 0.5 : 0.001 : 1.0; f = 60.0;
EIv=(X(16,1)+X(12,1)*xX(13,1)*x.^ 2/2X(14,1)*x.^ 3/6+ f *(x0.5).^ 3/6);
EIfi = (X(12,1) X(13,1)*x X(14,1)*x.^ 2 / 2+ f *(x0.5).^ 2/2);
Q = X(14,1) f; M = X(13,1) + X(14,1)*x f *(x0.5);
subplot (2,2,1), plot (x, EIv); axis ([0.5 1.0 3 3]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0.5 1.0 6 6]); grid on subplot (2,2,3), plot (x, Q); axis ([0.5 1.0 60 60]); grid on subplot (2,2,4), plot (x, M); axis ([0.5 1.0 25 25]); grid on Эпюры прогибов EI(x), углов поворота EI(x), поперечных сил Q(x) и изгибающих моментов М(x) представлены на рис. 3. Отметим, что эпюры прогибов EI(x) позволяют уточнить принятую деформированную схему рамы. Если сравнивать рис. 3.15, а и рис. 3.14, d, то видно, что деформированная схема практически совпадает с действительным состоянием рамы. Последнее обстоятельство в немалой степени обеспечивает достоверность ее расчета.
3.3.3. Определение частот собственных колебаний Узлы рамы по рис. 3.16 при свободных колебаниях будут иметь линейные и угловые перемещения, т. е. данная конструкция относится к классу свободных систем. При ее движении возникают силы инерции линейно подвижных стержней 0-1 и 1-2. Учет этих сил инерции можно выполнить по формуле [2] где (х) – прогиб несвободного стержня, связанного с линейно подвижным стержнем; М – сосредоточенная масса линейно подвижного стержня; а – координата сосредоточенной массы.
В узле 1 прикладываем сосредоточенную массу стержня 0-1 и половину массы стержня 1-2, в итоге получится M 1 = 1,5ml. В узле 2 прикладываем половину массы 1-2, т. е. M 2 = 0,5ml (рис. 3.16).
Для стержней 1-3 и 2-4 принимаем, что прогиб приближенно описывается функцией (х ) = cos (x / 2l ), тогда по формуле (3.29) будем иметь В дальнейшем, при динамическом анализе данной рамы, необходимо использовать увеличенные массы стержней 1-3 и 2-4, чтобы учесть возникающие силы инерции.
Частоты собственных колебаний рамы определяются из частотного уравнения краевой задачи А() = 0, где матрица А() берется из уравнения (3.28), где достаточно заменить фундаментальные функции статического изгиба на фундаментальные функции поперечных колебаний (3.14).
А*= При поиске частот принимается EI = m = l = 1, так что аргументы фундаментальных функций поперечных колебаний запишутся так:
стержни 0-1; 1- где частота собственных колебаний. Программа поиска частот принимает вид (обозначения переменных такие же как у неразрезной балки) п = 20; п1 = 300; am = 0.01; dam = 0.01; X = zeros (n1,1); Y = zeros (n1,1);
for m = 1 : n a = zeros (n,n); la1 = sqrt (am); la2 = sqrt (am); la3 = (2*am^2)^0.25;
la4 = (4*am^2)^0.25;
a(1,2) = (sinh (la1)+ sin(la1))/(2*la1); a(1,4) =(sinh (la1) sin(la1))/(2*la1^3);
a(2,2)=(cosh (la1)+cos(la1))/2; a(2,4)=(cosh(la1)cos(la1))/(2*la1^2); a(2,17)=1;
a(3,2) = la1^4*a(2,4); a(3,4) = a(1,2); a(3,8) = 1; a(3,18) = 1;
a(4,2) = la1^4*a(2,4); a(4,4) = a(2,2); a(4,9) = 1; a(4,20) = 1; a(5,10) = 1;
a(5,19) = 1; a(6,8) = (cosh(la2) cos(la2))/(2*la2^2);
a(6,9)=(sinh(la2)sin(la2))/(2*la2^3); a(6,17)=(sinh(la2)+sin(la2))/(2*la2);
a(7,8) = a(6,17); a(7,9) = a(6,8); a(7,12) = 1; a(7,17)=(cosh(la2)+cos(la2))/2;
a(8,8) = a(7,17); a(8,9) = a(6,17); a(8,13) = 1; a(8,17) = la2^4*a(6,9);
a(9,8) = a(8,17); a(9,9) = a(7,17); a(9,15) = 1; a(9,17) = la2^4*a(6,8);
a(10,10) = 1; a(10,14) = 1; a(11,12)=(sinh(la3)+sin(la3))/(2*la3);
a(11,13)=(cosh(la3)cos(la3))/(2*la3^2); a(11,14)=(sinh(la3)sin(la3))/(2*la3^3);
a(11,16)=(cosh(la3)+cos(la3))/2; a(12,12) = a(11,16); a(12,13) = a(11,12);
a(12,14)=a(11,13); a(12,16)= la3^4*a(11,14); a(13,1)=1; a(13,12)= a(12,16);
a(13,13)=a(12,12); a(13,14)=a(11,12); a(13,16)= la3^4*a(11,13); a(13,1)=1;
a(14,12)=a(13,16); a(14,13)=a(13,12); a(14,14)=a(12,12); a(14,15)=la3^4*a(12,13);
a(15,5) = 1; a(15,15) = 1; a(16,16)=(cosh(la4)+cos(la4))/2;
a(16,17)=(sinh(la4)+sin(la4))/(2*la4); a(16,18)=(cosh(la4)cos(la4))/(2*la4^2);
a(16,19)=(sinh(la4)sin(la4))/(2*la4^3); a(17,16)= la4^4*a(16,19);
a(17,17)=a(16,16); a(17,18)=a(16,17); a(17,19)=a(16,18); a(18,6) = 1;
a(18,16)= la4^4*a(16,18); ); a(18,17)=a(17,16);
a(18,18)=a(16,16); a(18,19)=a(16,17); a(19,7) = 1; a(19,16)= la4^4*a(16,17);
a(19,17)=a(18,16); a(19,18)=a(17,16); a(19,19)=a(16,16); a(20,11) = 1;
a(20,20) = 1; d = det(a); X(m,1) = am; Y(m,1) = d; am = am + dam; end;
plot (X,Y); grid оn format long e [X Y] (3.32) Результаты поиска частот собственных колебаний рамы
EI EI EI
3.3.4. Построение форм собственных колебаний Для построения форм собственных колебаний необходимо определить граничные параметры каждого стержня рамы при собственных колебаниях.Воспользуемся кинематическим способом возбуждения колебаний рамы и b(16,1) = EI 1-3( l ) = 1.
Разрешающее уравнение МГЭ рамы примет вид (3.32), где использован кинематический способ возбуждения собственных колебаний, т. е. задается одинаковое единичное смещение неподвижных опор рамы.
Программа определения относительных граничных параметров (они нормируются относительно начального перемещения EI 1-3(о) = EI 2-4(о)) в отдельном М-файле имеет следующий текст a = zeros (20,20); b = zeros (20,1); X = zeros (20,1); x = 2.946045;
la1 = sqrt (x); la2 = sqrt (x); la3 = (2*x^2)^0.25; la4 = (4*x^2)^0.25;
a(1,2) = (sinh (la1)+ sin(la1))/(2*la1); a(1,4) =(sinh (la1) sin(la1))/(2*la1^3);
a(2,2)=(cosh (la1)+cos(la1))/2; a(2,4)=(cosh(la1)cos(la1))/(2*la1^2); a(2,17)=1;
a(3,8)=1; a(3,2)=la1^4*a(1,4); a(3,4)=a(1,2); a(3,18)=1; a(4,2)=la1^4*a(2,4);
a(4,9)=1; a(4,4)=a(2,2); a(4,20)=1; a(5,10)=1; ; a(5,19)=1;
a(6,17)=(sinh(la2)+sin(la2))/(2*la2); a(6,8) = (cosh(la2) cos(la2))/(2*la2^2);
a(6,9)=(sinh(la2)sin(la2))/(2*la2^3); a(7,17)=(cosh(la2)+cos(la2))/2;
a(7,8) = a(6,17); a(7,9) = a(6,8); a(7,12) = 1; a(8,13) = 1;
a(8,17) = la2^4*a(6,9); a(8,8) = a(7,17); a(8,9) = a(6,17); a(9,15) = 1;
a(9,17) = la2^4*a(7,9); a(9,8) = a(8,17); a(9,9) = a(7,17); a(10,14) = 1;
a(10,10) = 1; a(11,12)=(sinh(la3)+sin(la3))/(2*la3);
a(11,13)=(cosh(la3)cos(la3))/(2*la3^2); a(11,14)=(sinh(la3)sin(la3))/(2*la3^3);
a(11,16)=(cosh(la3)+cos(la3))/2; a(12,12) = a(11,16); a(12,13) = a(11,12);
a(12,14)=a(11,13); a(12,16)= la3^4*a(11,14); a(13,1)=1;
a(13,12)= a(12,16); a(13,13)=a(11,16); a(13,14)=a(11,12);
a(13,16)= la3^4*a(11,13); a(14,12)=a(13,16); a(14,13)=a(12,16);
a(14,14)=a(11,16); a(14,3)=1; a(14,16)=la3^4*a(12,13);
a(15,5) = 1; a(15,15) = 1; a(16,16)=(cosh(la4)+cos(la4))/2;
a(16,17)=(sinh(la4)+sin(la4))/(2*la4); a(16,18)=(cosh(la4)cos(la4))/(2*la4^2);
a(16,19)=(sinh(la4)sin(la4))/(2*la4^3); a(17,16)= la4^4*a(16,19);
a(17,17)=a(16,16); a(17,18)=a(16,17); a(17,19)=a(16,18);
a(18,6) = 1; a(18,16)= la4^4*a(16,18); a(18,17)=a(17,16); a(18,18)=a(16,16);
a(18,19)=a(16,17); a(19,7) = 1; a(19,16)= la4^4*a(17,18); a(19,17)=a(18,16);
a(19,18)=a(17,16); a(19,19)=a(16,16); a(20,11) = 1; a(20,20) = 1;
b(11,1)= 1; b(16,1)= 1; X = a\ b; X = X / X(16,1) Значения относительных граничных параметров сведены в таблицу 3.7.
EI 01 (о ) =Х(2,1)= -0. Q 24 (l ) =Х(3,1)= 0.8623107 -10.5288 -155.6993 531.9886 110.7758 449. Q13 (l ) =Х(7,1)= 1.2201107 -14.9016 220.3370 -320.1960 234.3273 72. Q12 (о) =Х(9,1)= -0.4241107 5.1781 -60.1557 98.3542 103.1128 -443. 10 N 12 (о ) =Х(10,1)= -0.2047107 2.4991 -120.4025 420.3761 8.3444 241. 11 N 13 (l ) =Х(11,1)= 0.3733107 -4.5580 -10.8130 -923.8486 -90.5222 604. 14 Q 24 (о ) =Х(14,1)= 0.2047107 -2.4991 120.4025 -420.3761 -8.3444 -241. 15 N 24 (о ) =Х(15,1)= -0.4059107 4.9558 77.7808 81.0184 -69.5101 447. 18 М 13 (о ) =Х(18,1)= -0.2514107 3.0696 10.2235 -112.4317 29.4391 -27. 19 Q13 (о ) =Х(19,1)= -0.2047107 2.4991 -120.4025 420.3761 8.3444 241. N 13 (о ) =Х(20,1)= 0.3733107 -4.5580 -10.8130 -923.8486 -90.5222 604. Из таблицы следует, что единичные перемещения опор приводят к весьма значительным усилиям и перемещениям в элементах рамы.
Формы колебаний стержней рамы, построенные по соотношениям метода начальных параметров (3.17) представлены на рис. 3.17. Протокол их построения запишется следующим образом.
Протокол построения форм колебаний принимает вид х = 0 : 0.001 : 1.0; у = 2.946045; la1 = sqrt (у); la2 = sqrt (у);
la3 = (2*у^2)^0.25; la4 = (4*у^2)^0.25;
EIv1 = (X(2,1)*(sinh(la1*x)+sin(la1*x))/(2* la1) X(4,1)*(sinh(la1*x) … sin(la1*x))/(2* la1^3));
EIv2 = (X(17,1)*(sinh(la2*x)+sin(la2*x))/(2* la2) X(8,1)*(cosh(la2*x) … cos(la2*x))/(2* la2^2) X(9,1)*(sinh(la2*x)sin(la2*x))/(2* la2^3));
EIv3 = (X(16,1)*(cosh(la3*x)+cos(la2*x))/2+ X(12,1)*(sinh(la3*x)+ … sin(la3*x))/(2* la3) X(13,1)* (cosh(la3*x)cos(la*x))/(2* la3^2) … X(14,1)*(sinh(la3*x)sin(la3*x))/(2* la3^3));
Eiv4 = (X(16,1)*(cosh(la4*x)+cos(la4*x))/2 + X(17,1)*(sinh(la4*x) + … sin(la4*x))/(2*la4) X(18,1)*(cosh(la4*x)cos(la4*x))/(2* la4^2) … X(19,1)*(sinh(la4*x)sin(la4*x))/(2* la4^3));
subplot (2,2,1), plot (x, EIv1); axis ([0 1 0.1 0.1]); grid on subplot (2,2,2), plot (x, EIv2); axis ([0 1 0.1 0.1]); grid on subplot (2,2,3), plot (x, EIv3); axis ([0 1 1 1]); grid on subplot (2,2,4), plot (x, EIv4); axis ([0 1 1 1]); grid on Как и при расчете неразрезной балки, перед выполнением данного протокола в окно команд помещается вектор относительных граничных параметров Х, также необходимо корректировать масштабы форм колебаний.
Отметим, что данный расчет спектра частот и форм собственных колебаний кинематическим способом является основой расчета сооружений на сейсмические воздействия. Можно заключить, что представленная методика по МГЭ является наиболее простой и эффективной по сравнению с имеющимися [3].
3.3.5. Расчет на вынужденные колебания Для определения напряженно-деформированного состояния рамы при вынужденных колебаниях с частотой = 0.5 1 принимаем ориентированный граф по рис. 3.18 с тем, чтобы можно было использовать уже запрограммированную матрицу коэффициентов А, добавив к ней элементы матрицы нагрузки В по формулам (3.18 – 3.20). Коэффициенты фундаментальных функций с учетом сил инерции линейно подвижных стержней примут вид (3.31), где вместо нужно подставить частоту вынужденных колебаний. Согласно алгоритму МГЭ разбиваем раму на 4 стержня, нумеруем узлы и стрелками обозначаем начало и конец каждого стержня. Расчет рамы на вынужденные колебания сводится к определению динамических граничных параметров из матричного уравнения А Х = В. Программа решения данного уравнения записывается в отдельном М-файле. После решения матричного уравнения граничные параметры используются для построения эпюр напряженнодеформированного состояния стержней рамы. Для ухода от резонанса частоту вынужденных колебаний принимаем равной Разрешающее уравнение краевой задачи рамы для орграфа по рис. 3. принимает вид (3.33) Программа вычисления граничных параметров рамы имеет следующий текст q = 10.0; m = 20.0; f = 60.0;
a = zeros (20,20); b = zeros (20,1); X = zeros (20,1); t =0.5* 2.946045;
la1 = sqrt (t); la2 = sqrt (t); la3 = (2*t^2)^0.25; la4 = (4*t^2)^0.25;
a(1,2) = (sinh (la1)+ sin(la1))/(2*la1); a(1,4) =(sinh (la1) sin(la1))/(2*la1^3);
a(2,2)=(cosh (la1)+cos(la1))/2; a(2,4)=(cosh(la1)cos(la1))/(2*la1^2);
a(2,17)=1; a(3,8)=1; a(3,2)=la1^4*a(1,4); a(3,4)=a(1,2); a(3,18)=1;
a(4,2)=la1^4*a(2,4); a(4,4)=a(2,2); a(4,9)=1; a(4,20)=1; a(5,10)=1; a(5,19)=1;
b(1,1)= q*(cosh (la1)+cos(la1) 2)/(2*la1^4);
b(2,1)= q*(sinh (la1) sin(la1))/(2*la1^3);
b(3,1)= q*(cosh (la1) cos(la1))/(2*la1^2);
b(4,1)= q*(sinh (la1)+ sin(la1))/(2*la1);
a(6,17)=(sinh(la2)+sin(la2))/(2*la2); a(6,8) = (cosh(la2) cos(la2))/(2*la2^2);
a(6,9)=(sinh(la2)sin(la2))/(2*la2^3); a(7,17)=(cosh(la2)+cos(la2))/2;
a(7,8) = a(6,17); a(7,9) = a(6,8); a(7,12) = 1; a(8,13) = 1;
a(8,17) = la2^4*a(6,9); a(8,8) = a(7,17); a(8,9) = a(6,17); a(9,15) = 1;
a(9,17) = la2^4*a(6,8); a(9,8) = a(8,17); a(9,9) = a(7,17); a(10,14) = 1;
a(10,10) = 1; b(6,1)= q*(cosh (la2)+cos(la2) 2)/(2*la2^4);
b(7,1)= q*(sinh (la2) sin(la2))/(2*la2^3);
b(8,1)= q*(cosh (la2) cos(la2))/(2*la2^2);
b(9,1)= q*(sinh (la2)+ sin(la2))/(2*la2);
a(11,12)=(sinh(la3)+sin(la3))/(2*la3); a(11,13)=(cosh(la3)cos(la3))/(2*la3^2);
a(11,14)=(sinh(la3)sin(la3))/(2*la3^3); a(11,16)=(cosh(la3)+cos(la3))/2;
a(12,12)=a(11,16); a(12,13)=a(11,12); a(12,14)=a(11,13); a(12,16)=la3^4*a(11,14);
a(13,1)=1; a(13,12)= a(12,16); a(13,13)=a(11,16); a(13,14)=a(11,12);
a(13,16)= la3^4*a(11,13); a(14,12)=a(13,16); a(14,13)=a(12,16);
a(14,14)=a(11,16); a(14,3)=1; a(14,16)=la3^4*a(12,13);
a(15,5) =1; a(15,15)=1; b(11,1)=f*(sinh (la3/2) sin(la3/2))/(2*la3^3);
b(12,1)= f*(cosh (la3/2) cos(la3/2) 2)/(2*la3^2);
b(13,1)= f*(sinh (la3/2)+ sin(la3/2))/(2*la3);
b(14,1)= f*(cosh (la3/2) + cos(la3/2))/2;
a(16,16)=(cosh(la4)+cos(la4))/2; a(16,17)=(sinh(la4)+sin(la4))/(2*la4);
a(16,18)=(cosh(la4)cos(la4))/(2*la4^2); a(16,19)=(sinh(la4)sin(la4))/(2*la4^3);
a(17,16)= la4^4*a(16,19); a(17,17)=a(16,16); a(17,18)=a(16,17); a(17,19)=a(16,18);
a(18,6)=1; a(18,16)=la4^4*a(16,18); a(18,17)=a(16,17); a(18,18)=a(16,16);
a(18,19)=a(16,17); a(19,7) = 1; a(19,16)=la4^4*a(17,18); a(19,17)=a(18,16);
a(19,18)=a(17,16); a(19,19)=a(16,16); a(20,11) = 1; a(20,20) = 1;
b(16,1)= m*(cosh (3*la4/4) cos(3*la4/4))/(2*la4^2);
b(17,1)=m*(sinh (3*la4/4) + sin(3*la4/4))/(2*la4);
b(18,1)= m*(cosh (3*la4/4) + cos(3*la4/4))/2;
b(19,1)= m*(la4sinh (3*la4/4) sin(3*la4/4))/2;
X = a\ b Результаты расчета динамических граничных параметров рамы сведены в таблицу 3.8.
3.3.6. Построение эпюр напряженно-деформированного состояния Эпюры прогибов EI(x), углов поворота EI(x), поперечных сил Q(x), изгибающих моментов М(x) строятся по соотношениям метода начальных параметров (3.21), (3.22). Протоколы построения эпюр EI(x), EI(x), Q(x) и М(x) примут вид х = 0 : 0.001 : 1.0; t = 0.5*2.946045; la1 = sqrt (t); q = 10.0;
EIv = (X(2,1)*(sinh(la1*x)+sin(la1*x))/(2*la1) X(4,1)*(sinh … (la1*x)sin(la1*x))/(2*la1^3) +q*(cosh(la1*x)+cos(la1*x) 2) / (2*la1^4);
EIfi = (X(2,1)*(cosh(la1*x)+cos(la1*x)) / 2 X(4,1)* (cosh(la1*x) … cos(la1*x))/(2*la1^2)+q*(sinh(la1*x)sin(la1*x))/(2*la1^3));
Q = X(2,1)*la1^4*(cosh(la1*x)cos(la1*x)) / (2* la1^2) + … X(4,1)*(cosh(la1*x)+cos(la1*x)) / 2 q*(sinh(la1*x)+sin(la1*x))/(2*la1);
M = X(2,1)*la1^4*(sinh(la1*x)sin(la1*x))/(2* la1^3) + … X(4,1)*(sinh(la1*x)+sin(la1*x))/(2*la1)q*(cosh(la1*x)cos(la1*x))/(2*la1^2);
subplot (2,2,1), plot (x, EIv); axis ([0 1 0.6 0.6]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 1 2 2]); grid on subplot (2,2,3), plot (x, Q); axis ([0 1 12 12]); grid on subplot (2,2,4), plot (x, M); axis ([0 1 6 6 ]); grid on х = 0 : 0.001 : 1.0; t = 0.5*2.946045; la2 = sqrt (t); q = 10.0;
EIv = (X(17,1)*(sinh(la2*x)+sin(la2*x))/(2*la2) X(8,1)*(cosh (… la2*x)cos(la2*x))/(2*la2^3)X(9,1)*(sinh(la2*x)sin(la2*x)) / … (2*la2^3)+q*(cosh(la2*x)+cos(la2*x) 2) / (2*la2^4));
EIfi = (X(17,1)*(cosh(la2*x)+cos(la2*x)) / 2 X(8,1)*(sinh(la2*x)+ … sin(la2*x))/(2*la2) X(9,1)*(cosh(la2*x) cos(la2*x))/(2*la2^2)+ … q*(sinh(la2*x)sin(la2*x))/(2*la2^3));
Q =X(17,1)*la2^4*(cosh(la2*x)cos(la2*x))/(2*la2^2)+X(8.1)*la2^4* … (sinh(la2*x)sin(la2*x))/(2*la2^3)X(9,1)*(cosh(la2*x)+cos(la2*x)) / 2 … q*(sinh(la2*x)+sin(la2*x))/(2*la2);
M=X(17,1)*la2^4*(sinh(la2*x)sin(la2*x))/(2*la2^3)+X(8,1)*(cosh(la2*x)+… cos(la2*x)) / 2+X(9,1)*(sinh(la2*x)+sin(la2*x))/(2*la2)q*(cosh(la1*x) … cos(la2*x))/(2*la2^2);
subplot (2,2,1), plot (x, EIv); axis ([0 1 0.4 0.4]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 1 2 2]); grid on subplot (2,2,3), plot (x, Q); axis ([0 1 28 28]); grid on subplot (2,2,4), plot (x, M); axis ([0 1 12 12]); grid on х = 0 : 0.001 : 0.5; t = 0.5*2.946045; la3 = (2*t^2)^0.25;
EIv = (X(16,1)*(cosh(la3*x)+cos(la3*x))/2+X(12,1)*(sinh(la3*x)+ … sin(la3*x))/(2*la3) X(13,1)*(cosh(la3*x)cos(la3*x))/(2*la3^2) … X(14,1)*(sinh(la3*x)sin(la3*x))/(2*la3^3));
EIfi = (X(16,1)*la3^4*(sinh(la3*x) sin(la3*x))/(2*la3^3)+X(12,1)*( … cosh(la3*x)+cos(la3*x))/2X(13,1)*(sinh(la3*x)+sin(la3*x))/(2*la3) … X(14,1)*(cosh(la3*x)cos(la3*x))/(2*la3^3));
Q = X(16,1)*la3^4*(sinh(la3*x)+sin(la3*x))/(2*la3)X(12,1)*( … cosh(la3*x)cos(la3*x))*la3^4/(2*la3^2) + X(13,1)*la3^4*( … sinh(la3*x) sin(la3*x))/(2*la3^3)+X(14,1)*(cosh(la3*x)+cos(la3*x))/2;
M = X(16,1)*la3^4*(cosh(la3*x)cos(la3*x))/(2*la3^2)X(12,1)*la3^4* … (sinh(la3*x)sin(la3*x))/(2* la3^3)+X(13,1)*(cosh(la3*x)+cos(la3*x))/2+… X(14,1)*(sinh(la3*x)+sin(la3*x))/(2*la3);
subplot (2,2,1), plot (x, EIv); axis ([0 0.5 4 4]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 0.5 6 6]); grid on subplot (2,2,3), plot (x, Q); axis ([0 0.5 10 10]); grid on subplot (2,2,4), plot (x, M); axis ([0 0.5 10 10]); grid on б) При 0.5l х l м;
х = 0.5 : 0.001 : 1.0; t = 0.5*2.946045; la3 = (2*t^2)^0.25; f = 60.0;
EIv = (X(16,1)*(cosh(la3*x)+cos(la3*x))/2+X(12,1)*(sinh(la3*x)+ … sin(la3*x))/(2*la3) X(13,1)*(cosh(la3*x)cos(la3*x))/(2*la3^2) … X(14,1)*(sinh(la3*x)sin(la3*x))/(2*la3^3)+f*(sinh(la3*(x0.5)) … sinh(la3*(x0.5)))/(2*la3^3));
EIfi = (X(16,1)*la3^4*(sinh(la3*x) sin(la3*x))/(2*la3^3)+ … X(12,1)*(cosh(la3*x)+cos(la3*x))/2 X(13,1)*(sinh(la3*x)+sin(la3*x))/ … (2*la3) X(14,1)*(cosh(la3*x)cos(la3*x))/(2*la3^2)+ … f *(cosh(la3*(x0.5))cos(la3*(x0.5)))/(2*la3^2));
Q = X(16,1)*la3^4*(sinh(la3*x)+sin(la3*x))/(2*la3)X(12,1)*( … cosh(la3*x)cos(la3*x))*la3^4/(2*la3^2) + X(13,1)*la3^4*(sinh( … la3*x)sin(la3*x))/(2*la3^3)+X(14,1)*(cosh(la3*x)+cos(la3*x))/2 … f *(cosh(la3*(x0.5))+cos(la3*(x0.5)))/2;
M = X(16,1)*la3^4*(cosh(la3*x)cos(la3*x))/(2*la3^2)X(12,1)*la3^4* ( … sinh(la3*x)sin(la3*x))/(2*la3^3)+X(13,1)*(cosh(la3*x)+cos(la3*x))/2+… X(14,1)*(sinh(la3*x)+sin(la3*x))/(2*la3) f *(sinh(la3*(x0.5))+ … sinh(la3*(x0.5)))/(2*la3);
subplot (2,2,1), plot (x, EIv); axis ([0.5 1 4 4]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0.5 1 6 6]); grid on subplot (2,2,3), plot (x, Q); axis ([0.5 1 70 70]); grid on subplot (2,2,4), plot (x, M); axis ([0.5 1 30 30]); grid on х = 0 : 0.001 : 0.25; t = 0.5*2.946045; la4 = (4*t^2)^0.25;
EIv = (X(16,1)*(cosh(la4*x)+cos(la4*x))/2+X(17,1)*(sinh(la4*x)+ … sin(la4*x))/(2*la4) X(18,1)*(cosh(la4*x)cos(la4*x))/(2*la4^2) … X(19,1)*(sinh(la4*x)sin(la4*x))/(2*la4^3));
EIfi = (X(16,1)*la4^4*(sinh(la4*x)sin(la4*x))/(2*la4^3)+X(17,1)*( … cosh(la4*x)+cos(la4*x)/2X(18,1)*(sinh(la4*x)+sin(la4*x))/(2*la4) … X(19,1)*(cosh(la4*x)cos(la4*x))/(2*la4^2));
Q = X(16,1)*la4^4*(sinh(la4*x)+sin(la4*x))/(2*la4)X(17,1)*la4^4*( … cosh(la4*x)cos(la4*x))/(2*la4^2)+X(18,1)*la4^4*(sin(la4*x) … sinh(la4*x))/(2*la4^3)+X(19,1)*(cosh(la4*x)+cos(la4*x))/2;
M = X(16,1)*la4^4*(cosh(la4*x)cos(la4*x))/(2*la4^2)X(17,1)*la4^4*( … sinh(la4*x)sin(la4*x))/(2*la4^3)+X(18,1)*(cosh(la4*x)+cos(la4*x))/2+… X(19,1)*(sinh(la4*x)+sin(la4*x))/(2*la4);
subplot (2,2,1), plot (x, EIv); axis ([0 0.25 4 4]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 0.25 7 7]); grid on subplot (2,2,3), plot (x, Q); axis ([0 0.25 7 7]); grid on subplot (2,2,4), plot (x, M); axis ([0 0.25 18 18]); grid on х = 0.25 : 0.001 : 1.0; m = 20.0; t = 0.5*2.946045; la4 = (4*t^2)^0.25;
EIv =(X(16,1)*(cosh(la4*x)+cos(la4*x))/2+X(17,1)*(sinh(la4*x)+sin(la4*x)) … /(2*la4) X(18,1)*(cosh(la4*x)cos(la4*x))/(2*la4^2)X(19,1)*(sinh(la4*x) … sin(la4*x))/(2*la4^3)+m*(cosh(la4*(x0.25))cosh(la4*(x0.25)))/(2*la4^2));
EIfi = (X(16,1)*la4^4*(sinh(la4*x)sin(la4*x))/(2*la4^3)+X(17,1)*( … cosh(la4*x)+cos(la4*x))/2X(18,1)*(sinh(la4*x)+sin(la4*x))/(2*la4) … X(19,1)*(cosh(la4*x)cos(la4*x))/(2*la4^2)+m*(sinh(la4*(x0.25))+ … sin(la4*(x0.25)))/(2*la4));
Q = X(16,1)*la4^4*(sinh(la4*x)+sin(la4*x))/(2*la4)X(17,1)* la4^4*( … cosh(la4*x)cos(la4*x))/(2*la4^2)+X(18,1)*la4^4*(sinh(la4*x) … sin(la4*x))/(2*la4^3)+X(19,1)*(cosh(la4*x)+cos(la4*x))/2 … m*la4*(sinh(la4*(x0.25))sin(la4*(x0.25)))/2;
M = X(16,1)*la4^4*(cosh(la4*x)cos(la4*x))/(2*la4^2)X(17,1)*la4^4*( … sinh(la4*x)sin(la4*x))/(2*la4^3)+X(18,1)*(cosh(la4*x)+cos(la4*x))/2+… X(19,1)*(sinh(la4*x)+sin(la4*x))/(2*la4) m *(cosh(la4*(x0.25))+ … cos(la4*(x0.25)))/2;
subplot (2,2,1), plot (x, EIv); axis ([0.25 1 3 3]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0.25 1 7 7]); grid on subplot (2,2,3), plot (x, Q); axis ([0.25 1 13 13]); grid on subplot (2,2,4), plot (x, M); axis ([0.25 1 13 13]); grid on Эпюры EI(x), EI(x), Q(x), M(x) представлены на рис. 3. Если сравнить эпюры задачи статики и задачи динамики (рис. 3.15 и рис. 3.19), то видно, что приближение к первой частоте ( = 0.5 1 ) заметно влияет на напряженно-деформированное состояние стержней рамы. Особенно значительное различие наблюдается у стержней, к которым добавлены сосредоточенные массы для учета сил инерции линейно подвижных стержней.
Если же частота вынужденных колебаний далеко отстоит от первой частоты собственных колебаний и к стержням не добавляются сосредоточенные массы, то статический расчет практически не отличается от динамического расчета (см. расчет неразрезной балки).
3.3.7. Определение критических сил рамы Как и для неразрезной балки, определение критических сил потери устойчивости рамы сводится к поиску корней определителя матрицы устойчивости А*(F) = 0, где фундаментальные функции сжатых стержней определяются набором (3.23). Для рамы по рис. 3.20 при потере устойчивости стержень 1-3 испытывает продольно-поперечный изгиб, а другие стержни – статический изгиб. Принимая ориентированный граф расчета, аналогичный задачам статики и динамики (см. рис. 3.14, 3.16), можно существенно облегчить формирование матрицы устойчивости А*(F). Набор компенсирующих элементов остается неизменным, меняются лишь фундаментальные функции стержня 1-3 в сравнении с матрицей задачи статики (3.28). Далее организуется цикл вычисления определителя матрицы устойчивости А*(F), значения которого и сжимающую силу F выводят в окно команд. При просмотре этих массивов данных легко определить точки, где изменяет знак определитель d = А*(F). В этих точках и вычисляются критические силы потери устойчивости данной рамы. При поиске критических сил рекомендуется начинать вычисление определителя матрицы устойчивости с начального значения F0=0.01 EI с шагом F=0.01 EI. Число вычислений определителя п1=200 – позволяет надежно и достаточно точно определить первую и последующие критические силы. При решении задачи устойчивости рамы примем, что EI=1, l = 1 м.
Матрица устойчивости рамы примет вид (см. матрицу А* задачи статики (3.28)) Обозначения переменных в программе соответствуют задаче устойчивости неразрезной балки. Программа запишется следующим образом.
n = 20; n1 = 300; f = 0.01; df = 0.1; X = zeros(n1,1); Y= zeros(n1,1);
for m = 1 : n1 n2 = sqrt(f); a = zeros(n,n);
a(1,2)=1; a(1,4)=-1/6; a(2,2)=1; a(2,4)=-1/2; a(2,17)=-1;
a(3,4)=1; a(3,8)=-1; a(3,18)=-1; a(4,4)=1; a(4,9)=-1; a(4,20)=-1;
a(5,10)=-1; a(5,19)=1; a(6,8)=-1/2; a(6,9)=-1/6; a(6,17)=1; a(7,8)=-1;
a(7,9)=-1/2; a(7,12)=-1; a(7,17)=1; a(8,8)=1; a(8,9)=1; a(8,13)=-1;
a(9,9)=1; a(9,15)=-1; a(10,10)=1; a(10,14)=1; a(11,12)=1;
a(11,13)=-1/2; a(11,14)=-1/6; a(11,16)=1; a(12,12)=1; a(12,13)=-1;
a(12,14)=-1/2; a(13,1)=-1; a(13,13)=1; a(13,14)=1; a(14,3)=-1; a(14;14)=1;
a(15,5)=-1; a(15;15)=1; a(16,16)=1; a(16,17)=sin(n2)/n2;
a(16,18)=-(1-cos(n2))/(n2^2); a(16,19)=-(n2-sin(n2))/(n2^3);
a(17,17)=cos(n2); a(17,18)=-a(16,17); a(17,19)=a(16,18);
a(18,6)=-1; a(18,17)=n2*sin(n2); a(18,18)=a(17,17); a(18,19)=a(16,17);
a(19,7)=-1; a(19,19)=1; a(20,11)=-1; a(20,20)=1; d=det(a);
X(m,1)=f; Y(m,1)=d; f=f+df; end;
plot(X,Y); grid on format long e [X Y] Результаты поиска критических сил потери устойчивости рамы
EI EI EI
F4 = 132.291925 2 ; F5 = 212.607925 2 кН.3.3.8. Построение форм потери устойчивости рамы Для построения форм потери устойчивости рамы необходимо определить граничные параметры рамы при критических силах Fi. Для этого формируется система линейных уравнений краевой задачи по МГЭ, где в правой части в отличие от задачи определения форм собственных колебаний рамы необходимо положить только EI 1-3 ( l )=-1, т.е. b(16,1)=-1. Матричное уравнение данной задачи можно заимствовать из задачи динамики, где необходимо заменить фундаментальные функции поперечных колебаний на фундаментальные функции статического и продольно-поперечного изгибов при граничных значениях координаты х всех стержней (см. (3.35)). Уравнения форм потери устойчивости стержней рамы примут вид (3.12) и (3.26). Программа определения граничных параметров рамы запишется следующим образом:
a=zeros(20,20); b= zeros(20,1); X=zeros(20,1); f=15.098365;
n2=sqrt(f); a(1,2)=1; a(1,4)=-1/6; a(2,2)=1; a(2,4)=-1/2; a(2,17)=-1;
a(3,4)=1; a(3,8)=-1; a(3,18)=-1; a(4,4)=1; a(4,9)=-1; a(4,20)=-1;
a(5,10)=-1; a(5,19)=1; a(6,8)=-1/2; a(6,9)=-1/6; a(6,17)=1; a(7,8)=-1;
(3.35) a(7,9)=-1/2; a(7,12)=-1; a(7,17)=1; a(8,8)=1; a(8,9)=1; a(8,13)=-1; a(9,9)=1;
a(9,15)=-1; a(10,10)=1; a(10,14)=1; a(11,12)=1; a(11,13)=-1/2; a(11,14)=-1/6;
a(11,16)=1; a(12,12)=1; a(12,13)=-1; a(12,14)=-1/2; a(13,1)=-1; a(13,13)=1;
a(13,14)=1;
a(14,3)=-1; a(14,14)=1; a(15,5)=-1; a(15,15)=1; a(16,16)=1;
a(16,17)=sin(n2)/n2; a(16,18)=-(1-cos(n2))/(n2^2);
a(16,19)=-(n2-sin(n2))/(n2^3); a(17,17)=cos(n2); a(17,18)=-a(16,17);
a(17,19)=a(16,18); a(18,6)=-1; a(18,17)=n2*sin(n2); a(18,18)=(17,17);
a(18,19)=a(16,17); a(19,7)=-1; a(19,19)=1; a(20,11)=-1; a(20;20)=1;
b(16,1)=-1; X=a\b; X=X/X(16,1) Значения относительных граничных параметров сведены в таблицу 3.9.
п/п Формы потери устойчивости стержней рамы строятся по соотношениям метода начальных параметров. Для стержня 1-3 это выражение (3.26) продольно-поперечного изгиба, для остальных стержней – выражение (3.12) статического изгиба. Протокол построения форм потери устойчивости примет вид x=0 : 0.001 : 1.0; f=15.098365; n2=sqrt(f);
EIV1=-(X(2,1)*x-X(4,1)*x.^3/6);
EIV2=-(X(17,1)*х-X(8,1)*x.^2/2-X(9,1)*x.^3/6);
EIV3=-(X(16,1)+X(12,1)*x-X(13,1)*x.^2/2-X(14,1)*x.*3/6);
EIV4=-(X(16,1)+X(17,1)*sin(n2*x)/n2-X(18,1)*(1-cos(n2*x))/n2^2- … X(19,1)*(n2*x-sin(n2*x))/n2^3);
subplot (2,2,1), plot (x, EIV1); axis ([0 1 0.1 0.1]); grid on subplot (2,2,2), plot (x, EIV2); axis ([0 1 0.1 0.1]); grid on subplot (2,2,3), plot (x, EIV3); axis ([0 1 1 1]); grid on subplot (2,2,4), plot (x, EIV4); axis ([0 1 1 1]); grid on Перед выполнением данного протокола в окно команд необходимо поместить вектор относительных граничных параметров Х при соответствующей критической силе Fi и не забывать корректировать масштабы форм потери устойчивости. Сами формы потери устойчивости представлены на рис.
3.21.
В заключение отметим еще одну замечательную особенность матрицы А* задач динамики и устойчивости. Эта матрица описывает связь между собой всех начальных и конечных параметров элементов упругой системы. Поэтому при поиске и построении форм собственных колебаний и потери устойчивости можно назначать единичное значение произвольному элементу матрицы нагрузки В. Последующая нормировка компонентов вектора Х* относительно кинематических или статических параметров упругой системы приводит благодаря структуре матрицы А* к одинаковым соотношениям между начальными и конечными граничными параметрами. Другими словами, значения относительных параметров таблиц 3.9, 3.7 и т. д. остаются неизменными независимо от положения единичного элемента матрицы В. Небольшие ограничения касаются тех элементов матрицы В, которые приводят к делению на ноль или. Например, для уравнения (3.35) это b(19.1)=-1.
b(20,1)=-1 и некоторые другие.
3.3.9. Построение эпюр напряженно-деформированного состояния Примеры расчета неразрезной балки и рамы, приведенные в п. 3.1, 3.2, показывают, что необходимо учитывать разрывы 1-го рода при построении эпюр напряженно-деформированного состояния. Данная задача требует применения единичной функции Хевисайда Н(х-а) со сдвигом в точку а (рис.
3.22). В этом случае сокращаются тексты программ, а соответствующие эпюры можно получить сразу для всей длины стержня.
В MATLAB нет встроенной единичной функции Хевисайда, поэтому ее необходимо программировать с применением операторов условного перехода if … end и логического оператора аnd. В качестве примера представим программу, позволяющую построить эпюры прогиба EI(x), угла поворота EI(x), поперечных сил Q(x) и изгибающих моментов М(x) неразрезной балки по рис. 3.1, где выигрыш применения единичных функций наиболее очевиден. Для построения эпюр достаточно взять 500 – 600 точек. С учетом компонентов вектора Х формулы для параметров напряженнодеформированного состояния балки примут вид:
Как видно из этого выражения необходимо применить 5 разных единичных функций. В формуле (3.36) разности в квадратных скобках представляют собой реакции опор 1, 2, 3.
R2 = Q 1 2 (l ) Q 23 (o ) = X (3,1) X (12,1);
R3 = Q 23 (l ) Q 3 4 (o ) = X (5,1) X (16,1).
Изгибающий момент М(x), 0 х 14.0 м.
Программа, реализующая формулы (3.36) – (3.39), организует цикл, в котором применен блок формирования единичных функций. Далее следуют операторы построения эпюр и вывод таблицы значений параметров напряженно-деформированного состояния всей балки. Здесь необходимо подставлять численные значения элементов вектора граничных параметров Х.
Текст программы n=700; EIv=zeros(n,1); EIfi=zeros(n,1); Q=zeros(n,1);
M=zeros(n,1); x=0.0; dx=14.0/n; X=zeros(n,1);
for m = 1 : n if and(x > = 0.0, x < = 2.0) h1=0; h2=0; h3=0; h4=0; h5=0; end;
if and(x > 2.0, x < = 4.0) h1=1; h2=0; h3=0; h4=0; h5=0; end;
if and(x > 4.0, x < = 10.0) h1=1; h2=1; h3=0; h4=0; h5=0; end;
if and(x > 10.0, x < = 12.0) h1=1; h2=1; h3=1; h4=0; h5=0; end;
if and(x > 12.0, x < = 13.0) h1=1; h2=1; h3=1; h4=1; h5=0; end;
if and(x > 13.0) h1=1; h2=1; h3=1; h4=1; h5=1; end;
EIv(m,1)= – (–12.79835*x – (–1.04938)*x.^3/6 + 10*(x – 2).^2*h1 + … (–34.08093)*(x – 4).^3*h2/6 + 10*(x – 4).^4*h2/24 – … 10*(x – 10).^4*h3/24 + 4.36214*(x – 10).^3*h3/6 – … 40*(x – 12).^3*h4/6 + (–51.33059)*(x – 13).^3*h5/6);
EIfi(m,1)= – (–12.79835 – (–1.04938)*x.^2/2 + 20*(x – 2)*h1 + … (–34.08093)*(x – 4).^2*h2/2 + 10*(x – 4).^3*h2/6 – … 10*(x – 10).^3*h3/6 + 4.36214*(x – 10).^2*h3/2 – … 20*(x – 12).^2*h4 + (–51.33059)*(x – 13).^2*h5/2);
Q(m,1)= (–1.04938) – (–34.08093)*h2 10*(x – 4)*h2 + … 10*(x – 10)*h3 4.36214*h3 + 40*h4 – … ( 51.33059)*h5;
M(m,1)= (–1.04938)*x 20*h1-(34.08093)*(x – 4)*h2 … 5*(x – 4).^2*h2 + 5*(x – 10).^2*h3 (4.36214)*(x – 10)*h3 + … 40*(x – 12)*h4 ( 51.33059)*(x – 13)*h5; X(m,1)=x;
x = x + dx; end;
plot(X, EIv); grid on % plot(X, Q); grid on % plot(X, EIfi); grid on % plot(X, M); grid on [X EIv EIfi Q M] На экране компьютера появится окно с первой эпюрой, совпадающей с рис. 3.1, затем, последовательно снимая символ комментария %, строятся другие эпюры, а в окно команд выводится таблица численных значений кинематических и статических параметров балки. Добавим, что совершенно аналогично можно строить эпюры с разрывами 1-го рода и для любых стержней рам. Если вместо оператора plot(X, EIv) использовать оператор stem(X, EIv,filled), то получится закрашенная эпюра, мало отличающаяся от обычного представления такого рода графиков. Можно сократить число используемых единичных функций до одной. Однако в этом случае единичную функцию необходимо создать в отдельном М-файле, например, такого содержания function f = H(t) if t < = 0.0 f = 0.0; else f = 1.0; end;
Этот файл должен быть записан в рабочей папке MATLAB (обычно это папка work) под именем Н. Тогда в головной программе построения эпюр можно использовать одну единичную функцию, но от разных аргументов, т.е.
Н(х – 2), Н(х – 4), Н(х – 10) и т.д. Например, выражение для прогиба (3.36) примет вид EIv(m,1) = ( 12.79835*x ( 1.04938)*x.^3/6 + 10*(x – 2).^2*H(x – 2) + … ( 34.08093)*(x – 4).^3*H(x – 4)/6 + 10*(x – 4).^4*H(x – 4)/24 … 10*(x – 10).^4*H(x – 10)/24 + 4.36214*(x – 10).^3*H(x – 10)/6 … 40*(x – 12).^3*H(x – 12)/6 51.33059*(x – 13).^3*H(x – 13)/6).
При выполнении программы на экране компьютера появляется эпюра по рис. 3.1. Аналогично программируются другие параметры.
3.3.10. Варианты заданий Продолжение рис. 3. Окончание рис. 3. та задания Глава IV. Применение пакетов расширений МATLAB в 4.1. Пакет расширений Symbolic Math Система MATLAB является самой крупной системой компьютерной математики, ориентированной на матричные и численные вычисления. Однако MATLAB имеет также и средства аналитических вычислений. Пакет Symbolic Math Toolbox добавил системе MATLAB качественно новые возможности, связанные с выполнением символьных вычислений и преобразований, которые были доступны только в системе принципиально иного класса, относящихся к компьютерной алгебре. Теперь MATLAB, с учетом новых средств, становится в полной мере универсальной системой. Последняя реализация системы символьной математики Maple 6 в своем ядре и в расширениях имеет около 3000 функций. Система MATLAB с пакетом Symbolic, включающим в себя чуть больше сотни символьных команд и функций, намного уступает Maple по количеству таких команд и функций. Однако в данный пакет включены лишь наиболее важные и широко распространенные функции. Кроме того, есть специальная команда, которая дает доступ к ядру Maple, что заметно расширяет круг используемых функций.
Помимо типовых аналитических вычислений (таких как символьное дифференцирование и интегрирование, упрощение математических выражений, подстановка и т. д.) пакет Symbolic позволяет реализовать арифметические операции с произвольной точностью.
С помощью команды » help symbolic можно получить перечень входящих в пакет команд и функций. Для получения справки по любой команде или функции можно использовать команду » help sym / name.m, где name — это имя соответствующей команды или функции, а name.m — имя m-файла, задающего данную команду или функцию.
С демонстрационными примерами пакета Symbolic можно ознакомиться с помощью директории Symbolic Toolbox.
Создание символьных переменных Поскольку переменные системы MATLAB по умолчанию не определены и традиционно задаются как векторные, матричные, числовые и т. д., то есть не имеющие отношения к символьной математике, для реализации символьных вычислений нужно прежде всего позаботиться о создании специальных символьных переменных. В простейшем случае их можно определить как строковые переменные, заключив имена в апострофы. Например, » sin(x)^2 + соs(х)^ ??? Undefined function or variable ' x '.
» sin('x' )^2 + cos(' x ' ) ^ ans = В первом случае система MATLAB «возмутилась» нашей небрежностью и сообщила, что функция или переменная х не определена и ни о каких вычислениях синуса и косинуса речи быть не может. Вместе с тем она подсказала, как надо поступить — заключить имя переменной в апострофы, ибо таким образом система получает информацию о необходимости включить символьный режим вычислений. Поэтому во второй раз получен вполне осмысленный результат — сумма квадратов синуса и косинуса переменной 'х' выдана равной 1.
Функция создания символьных переменных sym Для работы с командами ядра Meple в MATLAB определён новый тип переменной sym – символьный объект. Фактически это строковые переменные. Для проведения аналитических (символьных) операций нужно, чтобы соответствующие переменные были предварительно объявлены.
• S = sym(A) — возвращает символьный объект S класса 'sym' для входного параметра А. Если А — строка, то будет получена символьная строка или символьная переменная, а если А — это число (скаляр) или матрица, то будут получены их символьные представления.
• х = sym ('x') — возвращает символьную переменную с именем 'х' и записывает результат в х.
Функция создания группы символьных объектов syms Для создания группы символьных объектов служит функция syms:
• syms argl arg2... — создает группу символьных объектов.
Функция создания списка символьных переменных findsym В математических выражениях могут использоваться как обычные, так и символьные переменные. Функция findsym позволяет выделить символьные переменные в составе выражения S:
• findsym(S) — возвращает в алфавитном порядке список всех символьных переменных выражения S. При отсутствии таковых возвращается пустая строка.
Примеры:
» findsym (a*x^2 + b*y + z) Функция вывода символьных выражений pretty MATLAB в отличие от современных систем MathCAD, Maple или Mathematica, пока не способна выводить выражения и результаты их преобразований в естественной математической форме с использованием общепринятых спецзнаков для отображения интегралов, сумм, произведений и т. д. Тем не менее некоторые ограниченные текстовым форматом возможности близкого к математическому виду вывода обеспечивает функция pretty:
• pretty(S) — дает вывод выражения S в формате, приближенном к математическому;
« x = sym(‘x’);
« pretty ( x^2 ) Арифметика произвольной точности Функция digits служит для установки числа цифр в числах арифметики произвольной точности. Она используется в следующем виде:
digits — возвращает число значащих цифр в числах арифметики произвольной точности (по умолчанию 32);
3. Для проведения вычислений в арифметике произвольной точности служит функция vpa:
R = vpa(S) — возвращает результат вычислений каждого элемента символьного массива S, используя арифметику произвольной точности с текущим числом цифр D, установленным функцией digits. Результат R имеет тип sym.
vpa(S,D) — возвращает результат вычислений каждого элемента массива S, используя арифметику произвольной точности с количеством знаков чисел D.
vpa(exp(1),50) 2. Символьные операции с выражениями 1. Функция упрощения выражений – simplify Функция simplify(S) поэлементно упрощает символьные выражения массива S. Если упрощение невозможно, то возвращается исходное выражение.
» V = [sin(х)^2 + соs(х)^2 log(a*b)] » simplify(V) » simplify((a^2 - 2*a*b + b^2) / (a - b)) Дополнительные возможности упрощения обеспечивает функция simple.
2. Функция расширения (раскрытия скобок) выражений - expand Функция expand(S) расширяет выражения, входящие в массив S. Рациональные выражения она раскладывает на простые дроби, полиномы — на полиномиальные разложения и т. д. Функция работает со многими алгебраическими и тригонометрическими функциями.
Примеры:
» S=[(x + 2)*(x + 3)*(x + 4) sin(2*x)];
» expand(S) [ x^3 + 9*x^2 + 26*x + 24, 2*sin(x)*cos(x)] » expand(sin(a + b)) sin(a)*cos(b) + cos(a)*sin(b) » expand((a + b)^:3) 3. Разложение выражения на простые множители – factor Функция factor(S) поэлементно разлагает выражения вектора S на простые множители, а целые числа — на произведение простых чисел. Следующие примеры иллюстрируют применение функции:
» x=sym('х');
» factor(x^7-l) » factor(sym('123456789')) (3 )^2*(3803)*(3607) 4. Комплектование по степеням – collect Функция co11ect(S,v) обеспечивает комплектование выражений в составе вектора или матрицы S по степеням переменной v.
5. Упрощение выражений – simple Функция simple(S) выполняет различные упрощения для элементов массива S и выводит как промежуточные результаты, так и самый короткий конечный результат. В другой форме — [R, HOW] = simple(S) – промежуточные результаты не выводятся.
6. Приведение к рациональной форме – numden Функция [N,D] = numden(A) преобразует каждый элемент массива А в рациональную форму в виде отношения двух неприводимых полиномов с целочисленными коэффициентами. При этом N и D — числители и знаменатели каждого преобразованного элемента массива.
Примеры:
» [n,d] = numden(sym(8/10)) 7. Обеспечение подстановок – subs Одной из самых эффектных и часто используемых операций символьной математики является операция подстановки. Она реализуется функцией subs, имеющей ряд форм записи:
• subs(S) — заменяет в символьном выражении S все переменные их символьными значениями, которые берутся из вычисляемой функции или рабочей области системы MATLAB.
• subs(S,NEW) — заменяет все свободные символьные переменные в S из списка NEW.
• subs (S, OLD, NEW) — заменяет OLD на NEW в символьном выражении S. При одинаковых размеров массивов OLD и NEW замена идет поэлементно. Если S OLD — скаляры, a NEW — числовой массив или массив ячеек, то скаляры расширяются до массива результатов.
» subs(sin(x) + cos(y), [x,y], [a,b]) 8. Обращение функции – finverse Часто возникает необходимость в задании функции, обратной по отношению к заданной функции f. Для этого в Symbolic имеется функция обращения inverse, которая задается в двух формах:
• g = finverse(f) — возвращает функцию, обратную к f. Считается, что f — функция одной переменной, например 'х'. Тогда g(f(x)) = х.
• g = finverse(f,v) — возвращает функцию, обратную к f, относительно заданной переменной v, так что g(f(v)) = v. Эта форма используется, если f — функция нескольких переменных.
Примеры:
» finverse(sinh(x)) » finverse(expCx)) 9. Суперпозиция функций – compose • compose(f, g) — возвращает f(g(y)), где f = f(x) и g = g(y). Независимые переменные х и у находятся с помощью функции findsym.
Символьные операции математического анализа 1. Функция вычисления производных – diff Для вычисления в символьном виде производных от выражения S служит функция diff, записываемая в формате d i f f ( S, ' v ' ) или diff(S, sym('v')) Она возвращает символьное значение первой производной от символьного выражения или массива символьных выражений S по переменной v. Эта функция возdS вращает • diff(S, n) — возвращает п-ю (п — целое число) производную символьного выражения или массива символьных выражений S по переменной v.
• diff(S, V, n) и diff(S, n, V) – возвращает п-ю производную S по переменной v.
» x =sym( 'x' ); y=sym( 'у' );
» slmplify(ans) » diff(s1n(y*x), x, 3) » diff([x^3 sin(x) exp(x)], x) 2. Функция интегрирования – int Функция int вычисляет неопределенные и определенные интегралы • int(S) — возвращает символьное значение неопределенного интеграла от символьного выражения или массива символьных выражений S по переменной, которая автоматически определяется функцией findsym. Если S — скаляр или матрица, то вычисляется интеграл по переменной 'х'.
• int(S, v) — возвращает неопределенный интеграл от S по переменной v.
• int(S, a, b) — возвращает определенный интеграл от S с пределами интегрирования от а до b, причем пределы интегрирования могут быть как символьными, так и числовыми.
• int(S, v, a, b) — возвращает определенный интеграл от S по переменной v с пределами от а до b.
Примеры:
» int(sin(x)^3, x) - l/3*sin(x)^2*cos(x)-2/3*cos(x) » int(log(2*x), x) » int((x^2-2)/(x*3-l), x, l, 2) » int((x^2-2)/(x*3-l), x, 2, 5) - 2/3*1og(2) + 2/3*1og(31) + 2/3*3^(l/2)*atan(11/3*3^(l/2)) -...
2/3*log(7) - 2/3*3^(1/2)*atan(5/3*3^(l/2)) » int([x^3 sin(x) exp(x)], x) » int(log(sin(x)),x,0,pi/2) С помощью функции int можно вычислять имеющие аналитическое решение сложные интегралы, например с бесконечными пределами (или одним из пределов), а также кратные интегралы.
» int(log(1+exp(-x),x,0,inf) » int(int(int(x^2 + y^2)*z, x, 0, a), y, 0, a), z, 0, a) 3. Функция вычисления пределов – limit Для вычисления пределов аналитически заданной функции F(x) служит функция limit limit(F, x, a) – возвращает предел символьного выражения F в точке х = а.
»limit(sin(x)/x, x, 0) 4. Функция разложения в ряд Тейлора – taylor taylor (f) – возвращает шесть членов ряда Маклорена функции f ;
taylor (f, n, x, a) – возвращает n членов ряда Тейлора в точке х = а ;
taylor (f, n) – возвращает ( n - 1) членов ряда Маклорена функции f ;
taylor (f, a) – возвращает шесть членов ряда Тейлора функции f в точке а.
»taylor(sin(x)) »taylor(int(sin(x)) 5. Функция вычисления суммы рядов – symsum Для аналитического вычисления суммы ряда служит команда symsum:
• symsum(S) — возвращает символьное значение суммы бесконечного ряда по переменной, найденной автоматически с помощью функции findsym;
• symsum(S, v) — возвращает сумму бесконечного ряда по переменной v;
• symsum(S, a, b) и symsum(S, v, a, b) — возвращают конечную сумму ряда в пределах номеров слагаемых от а до b.
» symsum(x^2) 1/3*x^3-1/2*x^2+1/6*x » symsum([x, x^2, x^3], 1, 5) 6. Функция решение алгебраических уравнений – solve Для решения систем алгебраических уравнений и одиночных уравнений служит функция solve:
• solve(expr1, expr2,... exprN, var1, var2,... varN) — возвращает значения переменных var1, при которых соблюдаются равенства, заданные выражениями exprI. Если в выражениях не используются знаки равенства, то полагается ехргI = 0;
• solve(expr1, expr2,..... exprN) — аналогична предшествующей функции, но переменные, по которым ищется решение, определяются функцией findsym.
Примеры решения уравнений:
» so1ve(x^3 -1, x) [ -1/2+1/2*i*3^(1/2)] [ -1/2-1/2*i*3^(1/2)] » solve(a*x^2+b*x+c) [ 1/2/a*(-b+(b^2-4*a*c)^(1/2))] [ 1/2/a*(-b-(b^2-4*a*c)^(1/2))] »S = solve(‘x+y=3’, ‘x*y^2=4’, x, y) » solve(‘sin(x)=0.5’, x) 0. 7. Решение дифференциальных уравнений – dsolve Для решения дифференциальных уравнений в форме Коши MATLAB имеет следующую функцию:
• dsolve( 'eqn1', 'eqn2',...) — возвращает аналитическое решение системы дифференциальных уравнений с начальными условиями. Они задаются равенствами eqnl (вначале задаются уравнения, затем начальные условия).
По умолчанию независимой переменной считается переменная 't', обычно обозначающая время. Можно использовать и другую переменную, добавив ее в конец списка параметров функции dsolve. Символ D обозначает производную по независимой переменной, то есть d / dt, при этом D2 означает d2 / dt2 и т. д. Имя независимой переменной не должно начинаться с буквы D.
Начальные условия задаются в виде равенств 'у(а)=b' или 'Dy(a)=b', где у — независимая переменная, а и b — константы. Если число начальных условий меньше, чем число дифференциальных уравнений, то в решении будут присутствовать произвольные постоянные С1, С2 и т. д.
Примеры применения функции dsolve:
» dsolve('D2x = -2*x') Cl*cos(2^(1/2)*t) + C2*sin(2^(l/2)*t) » dsolve('D2y = -2*x + y', 'у(0) = 1, 'х') (2*х*ехр(х) + (- С2 + 1)*ехр(х)^2 + С2 / ехр(х) Графические возможности Simbolic 1. Графопостроитель – funtool Команда funtool создает интерактивный графический калькулятор, позволяющий быстро построить две функции одной переменной - f(x) и g(x).
Например, одна может задавать собственно функцию, а другая — ее производную. Функции обозначаются как ' f = ' и ' g = ' и после знака равенства можно набрать функции с помощью клавиш калькулятора в его нижней части. С помощью полей 'х = ' и ' а = ' можно задать диапазон изменения переменной х и значение масштабирующего параметра а.
При запуске команды funtool появляются окна для двух функций и окно калькулятора (рис. 4.1). По умолчанию заданы функции f(x) = х и g(x) = 1, предел изменения х от -2 до 2 и а = 1/2.
Верхний ряд кнопок вычислителя относится только к функции f(x) и задает следующие операторы:
• df/dx — символьное дифференцирование f(x);
• int f — символьное интегрирование f(x) при наличии замкнутой формы;
• simple f — упрощение выражения, если таковое возможно;
• num f — выделение числителя рационального выражения;
• den f — выделение знаменателя рационального выражения;
• 1/f — замена f(x) на 1 / f(x);
• finv — замена f(x) инверсной функцией.
Второй ряд клавиш выполняет операции масштабирования и сдвига f(x) с применением параметра 'а'.
Третий ряд клавиш предназначен для осуществления бинарных операций над функциями f(x) и g(x).
Четвертый ряд клавиш служит для работы с памятью калькулятор и иных операций:
• Insert — помещает текущую функцию в список функций.
• Cycle — выполняет текущую функцию из списка.
• Delete — удаляет выделенную функцию из списка.
• Reset — устанавливает f, g, x, а и fxl i st в исходное состояние.
• Help — выводит описание калькулятора.
• Demo — запускает демонстрационный пример.
• Close — завершает работу с калькулятором.
Благодаря описанным средствам вычислитель позволяет задать инересующую вас функцию, выполнить ее преобразования (например, дифференцирование и интегрирование) и, наконец, построить график функции и результатов ее преобразования (рис.4.2) Рис. 4.2. Построения графиков некоторых функций Таким образом, графопостроитель funtool является весьма удобным средством визуализации графиков самых различных функций.
2. Графики поверхностей – ezsurf и ezsurfc Команда ezsurf служит для построения графиков поверхностей, задаваемых функциями двух переменных f(x, у):
• ezsurf(f) — построение поверхности f(x,y) с параметрами х и у, меняющимися по умолчанию от -2 до 2;
• ezsurf(f,domain) — построение поверхностиf(x,y) с пределами изменения х и у, заданными параметром domain;
• ezsurf(x,y,z) — построение поверхности, заданной параметрически зависимостями x(s, t), y(s, t), z(s, t) при s и t, меняющихся в интервале от -2 до 2;
• ezsurf(x,y,z,[smin, smax, tmin, tmax]) — построение поверхности, заданной параметрически зависимостями x(s, t), y(s, t), z(s, t) при s и t меняющихся в заданном интервале.
Следующий пример показывает действие этой команды:
» ezsurf(rea1(asec(x+i*y))) Рис. 4.3. Пример построения графика поверхности командой ezsurf Аналогичная по синтаксису записи группа команд ezsurfc строит еще и контурный график поверхности на плоскости, лежащей под поверхностью.
Применение средств символьных вычислений в теоретической механике 1. Решение задач статики. При решении некоторых прикладных задач оптимизации желательно получать аналитические решения, устанавливающие функциональную связь между заданными параметрами и величинами, подлежащими определению. В этом случае целевая функция оптимизации имеет явное выражение и её максимум (минимум) определяется обычными методами математического анализа. Рассмотрим следующую задачу о равновесии конструкции:
З а д а н и е С 1. Определение наибольшего значения реакции опоры для На плоскую конструкцию, состоящую из двух невесомых элементов, которые в точке С соединены шарниром, наложены следующие связи:
1) Шарнирная опора на катках в точке А ;
2) Неподвижный цилиндрический шарнир в точке D ;
3) Горизонтальный невесомый стержень в точке К.
На конструкцию действуют:
1. Пара сил с моментом М = 5 кН м ( в точке Е).
2. Две сосредоточенные силы F1 = 4 кН ( в точке D ) и F2 = 6 кН ( в точке А ), образующие в точке приложения с контуром элемента соответственно углы и 300.
3. Перпендикулярно к действующему участку ( участок АВ ) равномерно распределённая нагрузка интенсивности q = 10 кН / м.
Длина каждого участка конструкции а = 0,4 м. Предполагая, что конструкция находится в равновесии, найти при каком значении угла наклона силы F1 реакция шарнира D ( R D = X 2 + YD ) имеет наибольшее значение, опредеD лить для этого случая реакции всех связей и усилия в шарнире С.
Решение. Для определения искомых реакций можно поступить двояко, либо мысленно расчленить систему тел (составную конструкцию) на отдельные твёрдые тела и рассмотреть равновесие каждого из тел в отдельности, либо, применяя принцип затвердевания, вначале рассмотреть равновесие всей конструкции, а затем равновесие её отдельных элементов.
I. Прежде рассмотрим равновесие всей конструкции.
1. Объект равновесия – составная конструкция.
2. Связи – шарнирная опора на катках в точке А, неподвижный цилиндрический шарнир в точке D и невесомый стержень в точке К. На конструкцию действуют сосредоточенные силы F1 и F2, равномерно распределённая нагрузка интенсивности q и пара сил с моментом М.
3. Действие связей заменим их реакциями RA, ХD, УD, S. Равномерно распределённую нагрузку заменим равнодействующей Другие активные силы - F1, F2 и пара сил с моментом М.
4. Составим уравнения равновесия (система координат Dху изображена на рисунке 4.5a) Данная система уравнений содержит четыре неизвестных XD, УD, S, RA.
II. Рассмотрим теперь равновесие правой части конструкции ( Рис. 4.5b ).
1. Объект равновесия – правая часть конструкции.
2. Связи – шарнирная опора на катках в точке А и соединительный шарнир С. На объект равновесия действует сосредоточенная сила F2, равнодействующая распределённой нагрузки Q и реакции связей XC, УC, RA.
3. Составим уравнения равновесия ( в системе координат Cху ) Решаем систему шести уравнений (а) – (б) с шестью неизвестными, система (б) является независимой и содержит три неизвестных. Находим, из третьего уравнения из второго уравнения Решаем систему уравнений (а), для этого воспользуемся символьными вычислениями системы MATLAB Q1=solve('XD-6+4*sin(al)-8*0.5+S=0','YD-8*0.866+3.464-4*cos(al)=0',...
'-6*0.5*0.4+4*0.4*cos(al)+4*2*0.4*sin(al)-5-YD*0.4+3.464*0.4*0.866+S*0.4+XD*2*0.4=0'...
1. 5.96418 – 4.*sin(al) 3.46400 + 4.*cos(al) ((5.96418 – 4.*sin(al))^2 + (3.46400 + 4.*cos(al))^2)^(1/2) ans = Таким образом, наибольшее значение полная реакция RD шарнира D имеет при = - 1,04461, а наименьшее – при = 2,09698.
Реакции связей RA, XC, YC, S не зависят от угла величины а наибольшее и наименьшее значения полной реакции шарнира D соответственно равны :
Варианты задания. Конструкция состоит из двух частей. Установить, при каком значении угла для силы F1 реакция опоры А ( R A = X 2 + YA ) имеет наибольшее значение, определить для этого случая реакции всех опор, а также соединения С.
Необходимые для расчёта данные приведены в таблице 4. Схемы конструкций для каждого из вариантов приведены на рисунках 4.6 – 4. 2. Решение задач кинематики точки.
З а д а н и е К1. Определение скорости и ускорения точки по заданным По заданным уравнениям движения точки М x = x ( t ) y = y ( t ) установить вид её траектории и для момента времени t = t1 ( c ) найти положение точки на траектории, её скорость, полное, касательное и нормальное ускорение, а также радиус кривизны траектории.
Необходимые для решения данные приведены в таблице 4.2.
Пример выполнения задания К1. Движение точки на плоскости Оху Рис. 4. Рис. 4. Рис. 4. Задан момент времени t1 = 0.
Необходимо:
1. Определить уравнение траектории точки.
2. Построить траекторию и указать на траектории положение точки в заданный момент времени t1.
3. Для заданного момента времени найти скорость и ускорение точки, её тангенциальное и нормальное ускорения, значение радиуса кривизны траектории.
4. Векторы скорости и ускорения точки показать на траектории.
Решение.
1. Движение точки задано координатным способом.
Для определения уравнения траектории точки, соответствующей заданным уравнениям движения, исключаем время t. Из уравнений движения имеем Применяя основное тригонометрическое тождество, получим Таким образом, траекторией точки является эллипс.
2. Определяем скорость точки При малых значениях t ( 0 < t syms x y t Vx Vy V ax ay a at an ro >> x=4*t;y=16*t^2-1;
>> [Vx,Vy,V,ax,ay,a,at,an,ro]=cinema2d(x,y) 4*(1+64*t^2)^(1/2) 256*t/(1+64*t^2)^(1/2) 32*(1-64*t^2/(1+64*t^2))^(1/2) 1/32*(16+1024*t^2)/(1-64*t^2/(1+64*t^2))^(1/2) Вид траектории определяется командой - ezplot(x2-1,-2,2) В момент времени t1 = 0,5 c кинематические характеристики определяются командой подстановки - subs >> K=[Vx,Vy,V,ax,ay,a,at,an,ro];
Columns 1 through 4.0000 16.0000 16.4924 0 32.0000 32. Columns 7 through 31.0446 7.7611 35. Аналогичным образом решается задача кинематики точки при её движении в пространстве: в этом случае m – функция задаётся в виде function[Vx,Vy,Vz,V,ax,ay,az,a,at,an,ro]=cinema3d(x,y,z);
%кинематика точки %входные параметры - х(t), y(t), z(t) - уравнения движения точки Vx=diff(x,'t');
Vy=diff(y,'t');
Vz=diff(z,'t');
V=sqrt(Vx^2+Vy^2+Vz^2);
ax=diff(Vx,'t');
ay=diff(Vy,'t');
az=diff(Vz,'t');
a=sqrt(ax^2+ay^2);
at=(Vx*ax+Vy*ay+Vz*az)/V;
an=sqrt(a^2-at^2);
ro=V^2/an;
Рассмотрим пример >> x=-2*t^2+3;y=-5*t;z=3*t;
>> [Vx,Vy,Vz,V,ax,ay,az,a,at,an,ro]=cinema3d(x,y,z) 16*t/(16*t^2+34)^(1/2) 4*(1-16*t^2/(16*t^2+34))^(1/2) 1/4*(16*t^2+34)/(1-16*t^2/(16*t^2+34))^(1/2) Траектории движения – ezplot3(x,y,z,[0,8]) В момент времени t1 = 0,5 c >> K=[Vx,Vy,Vz,V,ax,ay,az,a,at,an,ro];
Columns 1 through -2.0000 -5.0000 3.0000 6.1644 -4.0000 Columns 7 through 0 4.0000 1.2978 3.7836 10. варианта 4.2. Пакет оптимизации Optimization Toolbox Назначение и возможности пакета Пакет оптимизации (Optimization Toolbox) — это библиотека функций, расширяющих возможности системы MATLAB по численным вычислениям и предназначенная для решения задач оптимизации и систем нелинейных уравнений. Поддерживает основные методы оптимизации функций ряда переменных [ 6 ]:
• Безусловная оптимизация нелинейных функций.
• Метод наименьших квадратов.
• Решение нелинейных уравнений.
• Линейное программирование.
• Квадратичное программирование.
• Условная минимизация нелинейных функций.
• Методы минимакса.
• Многокритериальная оптимизация.
Рассматриваемый пакет дает возможности решать задачи минимизации функций, нахождения решений уравнений, задачи аппроксимации («подгонки» кривых под экспериментальные данные).
Различные типы таких задач вместе с применяемыми для их решения функциями пакета Optimization Toolbox приведены в таблице 4.3.
Типы задач, решаемых средствами пакета Optimization Toolbox
MATLAB
(одномерная) минимизация Тип задачи Математическая запись ФункцияMATLAB
(без ограничений) программирование А х < b, Aeq x=beq, Квадратичное min программирование Минимизация min f(x) при условиях fmincon при наличии ограничений Достижение цели min при условиях Полубесконечная min f (х) при условиях fseminf минимизацияMATLAB
уравнения Нелинейное уравнение одной переменной уравнения многих п неизвестных переменных Наименьших квадратов (МНК) линейный МНК при наличии Принятые обозначения:• а — скалярный аргумент; х, у — в общем случае векторные аргументы;
• f(а), f(х) — скалярные функции; F(x), c(x), ceq(x), K(x, w) — векторные функции;
• А, Aeq, С, Н — матрицы;
• b, beq, d, f, w, goal, xdata, ydata — векторы;
• xL, xU — соответственно, нижняя и верхняя границы области изменения аргумента.
Применяемые алгоритмы В рамках пакета Optimization Toolbox все задачи оптимизации делятся на две группы: задачи малой и средней размерности и задачи большой размерности. Соответственно, такое же деление принято для алгоритмов решения данных задач. Разумеется, это не означает, что для решения задач средней размерности нельзя применять алгоритмы большой размерности и наоборот. Просто алгоритмы той или иной группы более эффективны для задач своей размерности.
В пакете Optimization Toolbox реализован широкий набор алгоритмов для решения задач оптимизации средней и малой размерности. Основными для задач без ограничений являются симплексный метод Нелдера—Мида и квазиньютоновские методы. Для решения задач с ограничениями, минимакса, достижения цели и полубесконечной оптимизации использованы алгоритмы квадратичного программирования.
Задачи, сводящиеся к нелинейным МНК, решаются с помощью алгоритмов Ньютона—Рафсона и Левенберга—Марквардта.
Вспомогательные процедуры одномерной (скалярной) оптимизации используют алгоритмы квадратичной (параболической) и кубической интерполяции.
Рассмотрим наиболее часто применяемые в задачах оптимизации функции линейного (linprog) и нелинейного (fmincon) программирования.
1. Функция fmincon furincon — функция поиска минимума скалярной функции многих переменных при наличии ограничений вида - задача нелинейного программирования.
Функция записывается в виде х = fmincon( fun, x0, A, b) х = fmincon( fun, x0, A, b, Aeq, beq) х = fmincon( fun, x0, А, b, Aeq, beq, lb, ub) х = fmincon( fun, x0, A, b, Aeq, beq, lb, ub, nonlcon) х = fmincon( fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options) х = fmincon( fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, opt1ons, Pl, P2,...) [x, fval] = fmincon(...) – возвращается не только оптимальное значение векторного аргумента, но и значение целевой функции в точке минимума fval ;
[x, fval, exitflag] = fmincon(…) – то же, что и предыдущая функция, но возвращается ещё информация о характере завершения вычисления exitflag ;
[x, fval, exitflag, output] = fmincon(…) – то же, что и предыдущая функция, но возвращается ещё информация о результатах оптимизации (выходная структура) output ;
[x, fval, exitflag, output, lambda] = fmincon(...) – то же, что и предыдущая функция, но возвращаются ещё множители Лагранжа lambda ;
[x, fval, exitflag, output, lambda, grad] = fmincon(…) – то же, что и предыдущая функция, но возвращается ещё величина градиента функции в точке минимума grad ;
[x, fval, exitflag, output, lambda, grad, hessian] = fmincon(…) – то же, что и предыдущая функция, но возвращается ещё величина гессиана H функции в точке минимума.
Аргументы функции:
fun — векторная функция векторного аргумента. Должна быть задана либо с помощью функции inline, например:
» fun = inline('sin(x.*x)'):
либо как m -файл, например:
function F = myfun(x) nonlcon — функция, возвращающая значения функций-ограничений, а при необходимости (если задано options = optimset('Grad-Constr', 'on')) и их градиентов; должна быть оформлена в виде m-файла, например:
function [c, ceq] = mycon(x) с =... % Вычисление левых частей нелинейных неравенств ceq =... % Вычисление левых частей нелинейных равенств function [с, ceq, GC, GCeq] = mycon(x) с =... % Вычисление левых частей нелинейных неравенств ceq =... % Вычисление левых частей нелинейных равенств GC =... % Градиенты неравенств GCeq =... % Градиенты равенств Options – опции (их можно изменять, используя функцию optimset):
о DerivativeCheck — задает проверку соответствия производных, определенных пользователем, их вычисленным оценкам в виде первых разностей;
o Diagnostics — вывод диагностической информации о минимизируемой функции;
o DiffMaxChange — максимальные значения изменений переменных при определении первых разностей;
о DiffMinChange — минимальные значения изменений переменных при определении первых разностей;
o Display — уровень отображения: 'off' — вывод информации отсутствует, 'iter' — вывод информации о поиске решения на каждой итерации, 'final' — вывод только итоговой информации;
o Goal ExactAchieve — определяет количество целей, которые должны быть достигнуты «точно»;
о GradConstr — использование градиентов для ограничений (опция имеет смысл в случае применения аргумента nonlcon ), возможные значения — 'off и 'on';
o GradObj — использование градиента для целевой функции, определяемого пользователем (возможные значения — 'off и 'on');
о MaxFunEvals — максимальное число вычислений функции;
о Maxlter — максимальное допустимое число итераций;
o MeritFunction — устанавливает вид функции оценки качества достижения цели (возможные значения 'multiobj' или 'singleobj');
o TolCon — допуск останова вычислений при нарушении ограничений;
o TolFun — допуск останова вычислений по величине изменений функции;
o ТоlХ — допуск останова вычислений по величине изменений х;
• exitflag — информация о характере завершения вычислений: если эта величина положительна, то вычисления завершились нахождением решения х, если она равна нулю, то останов произошел в результате выполнения предельного числа итераций, если данная величина отрицательна, то решение не найдено;
• lambda — множители Лагранжа, соответственно, для различных типов ограничений:
o 1ambda. lower — для нижней границы lb;
o 1ambda. upper — для верхней границы ub;
o 1ambda. ineqlin — для линейных неравенств;
o 1ambda. eqlin — для линейных равенств;
o 1ambda. ineqnonlin — для нелинейных неравенств;
o 1ambda. eqnonlin — для нелинейных равенств;
• output — информация о результатах оптимизации:
o output. iterations — число выполненных итераций;
o output. funcCount — число вычислений функции;
o output. algorithm — используемый алгоритм;
o output. cgiterations — число PCG-итераций (только при использовании алгоритма большой размерности);
o output. stepsize — величина конечного шага поиска (только при использовании алгоритма средней размерности);
o output. firstorderopt — мера оптимальности первого порядка (норма вектора градиента в точке минимума) — только при использовании алгоритма большой размерности) o LargeScale — может принимать значения 'off' (по умолчанию) и 'on'. В первом слу;чае используется алгоритм средней размерности, во втором — алгоритм большой размерности.
Следующие опции используются только при работе с алгоритмом средней размерности (описание см. выше):
o DerivativeCheck;
o DiffMaxChange;
o DiffMinChange;
o LineSearchType — задание вида алгоритма одномерной оптимизации.
Опции, используемые только в алгоритме большой размерности:
o Hessian — гессиан (в случае матрицы Гессе, задаваемой пользователем, — см. выше);
o HessPattern — задание гессиана как разреженной матрицы (это может привести к существенному ускорению поиска минимума);
o MaxPCGIter — максимальное число итераций PCG-алгоритма (preconditioned conjugate gradient, см. выше);
o PrecondBandWidth — верхняя величина начальных условий для PCG-алгоритма;
o TolPCG — допуск на завершение PCG-итераций;
o TypicalX — типовые величины х;
5) возвращаемая величина output в данном случае имеет дополнительные компоненты:
Пример. Пусть требуется найти минимум функции f(х) = - x1 x2 x3 при начальном значении х = [10; 10; 10] и при наличии ограничений Решение. Вначале создадим m-файл, определяющий целевую функцию:
function f = myfun(x) f =- x(l)*x(2)*x(3) ;
Затем запишем ограничения в виде неравенств:
или в матричной форме:
Теперь нахождение решения представляется следующим образом:
» А = [-1 -2 -2; 1 2 2];
» b = [0; 72];
» x0 = [10; 10; 10]; % Стартовое значение » [x, fval] = fmincon('myfun', x0, A, b) 24. 12. 12. fval = -3.4560e+ 2. Функция linprog Функция linprog обеспечивает решение задачи линейного программирования.
Функция записывается в виде х = linprog(f,A,b,Aeq,beq) х = linprog(f,A,b,Aeq,beq,lb,ub) х = linprog(f,A,b,Aeq,beq,lb,ub,x0) х = linprog(f,A,b,Aeq,beq,lb,ub,x0,options) [x,fval] = linprog(...) [x,fval,exitflag] = linprog(...) [x,fval,exitflag,output] = linprog(...) [x,fval,exitflag,output,lambda] = linprog(...) Аргументы и возвращаемые величины здесь аналогичны рассмотренным ранее для функции fmincon одним исключением: здесь введен дополнительный аргумент f — вектор коэффициентов линейной целевой функции. Функция может использовать алгоритм большой размерности lipsol или алгоритм средней размерности (метод проекций).
Рассмотрим примеры задач оптимизации.
Пример 1. Требуется найти решение задачи, описывающейся соотношениями (задача линейного программирования) >> f = [ -5, -4, -6]; % Вектор коэффициентов линейной целевой функции >> % Матрица коэффициентов ограничений – неравенств >> A = [1 -1 1; 3 2 4; 3 2 0];
>> b = [20; 42; 30]; % Вектор ограничений – неравенств >> % Задание нижних границ переменных (нулей) >> lb = zeros(3,1);
>> % Поиск решения >> [x, fval, exitflag, output, lambda] = linprog(f, A, b, [ ], [ ], lb) Optimization terminated successfully 0. 15. fval = -78. exitflag = output = iterations: cgiterations: algorithm: 'lipsol' lambda = ineqlin: [3x1 double] eqlin [0xl double] upper [3x1 double] lower [3x1 double] Из возвращенной информации, в частности, следует:
• что оптимальное решение х = [0.0000 15.0000 3.0000];
• минимальное значение целевой функции равно -78.0000;
• вычисления завершились нахождением решения (exitflag>0);
• всего было выполнено 6 итераций;
• был использован алгоритм lipsol Пример 2. Цех малого предприятия должен изготовить 100 изделий трёх типов. Каждого изделия нужно сделать не менее 20 штук. На изделия уходит соответственно 4, 3,4 и 2 кг металла при его общем запасе 340 кг, а также по 4,75, 11 и 2 кг пластмассы при её общем запасе 700 кг. Сколько изделий каждого типа х1, х2 и х3 надо выпустить для получения максимального выпуска в денежном выражении, если цена изделия составляет по калькуляции 40, 30 и 20 грн ?
Данная задача сводится к вычислению максимума функции f(x1,x2,x3) = 40 x1 + 30 x2 + 20 x при наличии ограничений Имеем задачу линейного программирования.
f=[40;30;20];% Вектор коэффициентов линейной целевой функции % Матрица коэффициентов ограничений – неравенств A=[4 3.4 2; 4.75 11 2];
b=[340;700]; %вектор ограничений – неравенств Aeq=[1 1 1]; %матрица коэффициентов ограничений – равенств beq=[100]; %вектор ограничений – равенств % Задание нижних границ переменных (нулей) lb=[20;20;20];
% Поиск решений, определение максимум целевой функции [x,fval,exitflag,output]=linprog(-f,A,b,Aeq,beq,lb) Optimization terminated successfully 56. 20. 24. fval = -3.3200e+ exitflag = output = iterations: cgiterations: algorithm: 'lipsol' Замечание. Максимизация с помощью функций fminbnd, fminunc, fminsearch, linprog, fmincon, fgoalattain, fminimax, lsqcurvefit и lsqnonlin достигается путём замены целевой функции f(x) на - f(x) ( т. е. изменить знак). Аналогично в задачах квадратичного прогаммирования: там необходимо заменить матрицу Н и вектор f соответственно на - Н и -f.
Пример 3. Оптимальный наклон плоского солнечного коллектора [ 9 ].
Расчёт интенсивности солнечной радиации проводится по следующей формуле [Даффи, Бекма] где Hb, Hd – прямая и рассеянная составляющие солнечной радиации на () - приведенная поглощательная способность d - диффузная отражательная способность. Для системы покрытий из одного, двух, трёх и четырёх листов стекла d приблизительно равна 0,16, 0,24, 0,29, 0, - направленная поглощательная способность поглощающей пластины - пропускательная способность r - пропускательная способность без учёта поглощения n1, n2 - показатели преломления;
1, 2 - углы падения и преломления (рис.1).
Если среда 1 воздух, то n1 = 1, а показатель преломления стекла n2 = 1,526.
a - пропускательная способность, учитывающая только поглощение;
L - толщина стекла, см;
К - коэффициент ослабления стекла, см -1.
Значение К изменяется от 0,04 см –1 для высокопрозрачного стекла (с бесцветной кромкой) до 0,32 см –1 для плохого (с зеленоватым оттенком кромки) стекла. Для обычного оконного стекла L = 0,32 см, К = 0,161 см –1.
Подставляя соответствующие значения = 0,94, d = 0,16 (для одиночного покрытия), a = 0,9498 в ( 3 ) и ( 2 ), для приведенной поглощательной способности получим Величина Rb в выражении ( 1 ) определяется следующим соотношением здесь 1, z - углы падения соответственно для наклонной и горизонтальной Рис. 2. Ход лучей для горизонтальной и наклонной поверхностей.
- широта местности ( для Одессы = 46,50 );
- угол между рассматриваемой плоскостью и горизонтальной поверхностью ( наклон коллектора);
- склонение (угловое положение Солнца);
n - порядковый номер дня года;
- время суток в часах.
Таким образом, интенсивность солнечной радиации, определяемая формулой (1) с учётом зависимостей (5), (6), (8) – (11) представляет собою функцию времени года n, времени суток и угла наклона солнечного коллектора При расчёте оптимального значения угла наклона рассматривается усреднённое значение по времени интенсивности солнечной радиации за период использования системы горячего водоснабжения (в течение года апрель – октябрь и в течение суток с 6 часов утра до 19 часов вечера).
здесь i – месяцы года (апрель – октябрь); k – время суток (от 6 до 19 часов).
Для условий города Одессы ( = 46,5 градусов северной широты) прямая Hb и рассеянная Hd составляющие солнечной радиации на горизонтальной поверхности определялись по метеорологическим данным. Результаты расчётов средних значений солнечной инсоляции (12) в зависимости от наклона коллектора представлены на рис. Рис. 3. Зависимость средней интенсивности солнечной радиации от угла наклона приёмной пластины коллектора S ср ( ) Функция Sср = S () представляется следующей зависимостью По данной зависимости находим значение угла наклона солнечного коллектора, при котором средняя солнечная инсоляция имеет наибольшее значение, используя функцию fminbnd.
Вначале создаётся m – файл:
%Вычисление максимального угла наклона солнечного коллектора %Суммарная инсоляция function S=maxins(x) a=5.6667347;b=0.0018194204;c=0.017329411;
d=-0.00014610493;e=-0.0009964058;f=-1.3104374e-07;
m=a+c*x+e*x.^2;
n=1+b*x+d*x.^2+f*x.^3;
S=exp(n/m);
Теперь определяются значения угла:
>> x=fminbnd('maxins',0,48) 27. Для данного угла наклона * = 27,50 интенсивность солнечной радиации для плоскости приёмной пластины коллектора за весь период эксплуатации системы горячего водоснабжения (апрель – октябрь с 6 утра до 19 часов вечера) имеет наибольшее значение.
Пример 4. Расчёт оптимальной конструкции противоточной Выбор рациональной конструкции градирни и режимам её работы как тепломассообменного аппарата должен опираться на достаточно большом числе вариантов принятия решений.
Эта задача не однозначна и её решение зависит от целевой оптимизации. Однако в большинстве известных задач оптимизации тепломассообменных аппаратов предпочтение отдаётся экономическим критериям, как наиболее объективным.
В общем случае целевая функция оптимальности П тепломассообменного аппарата имеет вид здесь площадь тепломассообменной поверхности F(x,y,…,k), мощности NН (x,y,…,k), NB (x,y,…,k) насоса и вентилятора, GH(x,y,…,k), GB(x,y,…,k) – массовые потери теплоносителей зависят от ряда независимых величин x,y,…,k.
Для определения эксплуатационных и конструктивных параметров градирни экономическая эффективность может быть записана в виде где Э – годовые эксплуатационные расходы, грн / год; – капиталовложения, приходяно щиеся на один год нормативного срока окупаемости но, грн / год;
Минимизацию целевой функции П следует проводить при ограничении, устанавливающем зависимость температуры жидкости на выходе из аппарата от независимых переменных задачи В работе [ ] температура жидкости на выходе из аппарата представляется следующим образом:
Данная задача оптимизации П при ограничении ( 2 ), как задача нелинейного программирования, решается с помощью функции fmincon:
формирование целевой функции:
function f=zfun1(x) %задание целевой функции %стоимость градирни за год эксплуатации % х(1)- nL (произведение количества листов и длины насадки),м % х(2) - Н (высота насадки),м, х(3)- Gmg (массовый расход воздуха),кг/ч Se=0.15; %стоимость электроэнергии, грн/(кВт ч) Sv=1.0; %стоимость воды, грн/м tg=7000; %годовое число часов работы оборудования tno=5; %нормативный срок окупаемости, 5 лет del=0.0002; %толщина алюминиевого листа насадки - 0,2 мм ro=2710; %плотность алюминия, кг/м Kcto=23800; %цена единицы массы аппарата, грн/т Gl=200; %объемный расход жидкости, м3/ч rol=998.2; %плотность жидкости, кг/м Gml=Gl*rol; %массовый расход жидкости, кг/ч qul=4.5; %плотность орошения, кг/(м2 с) rog=1.165; %плотность газа, кг/м nug=16e-6; %кинематическая вязкость газа, м2/с Sv=Gml/3600/qul; %площадь горизонтального сечения насадки, м Vg=x(3)/3600/rog/Sv; %скорость газа, м/с b=Sv/x(1); %расстояние между листами, м def=2*b; %эквивалентный диаметр, м Re=x(3)/1800/rog/nug/x(1); %число Рейнольдса lam=0.3164/Re^0.25; %коэффициент сопротивления dp=lam*rog*Vg^2*x(2)/2/def; %гидравлические потери, Па P1=0.35*Kcto*ro*del*x(1)*x(2)/1000; %стоимость аппарата av=599.426; %коэффициенты, определяющие удельные капиталовложения bv=0.3058; % вентиляторы в зависимости от его мощности ett=0.7; %к.п.д. вентилятора Nv=x(3)*dp/3600/rog/ett/1000; %мощность вентилятора P3=Nv*(0.35*av*Nv^bv+Se*tg); %стоимость вентилятора, монтажа и годовой эксплуатации формирование ограничения:
function [c, ceq]=confuneq2(x) %нелинейное ограничение %задана температура жидкости на выходе из аппарата % х(1)- nL (произведение количества листов и длины насадки) % х(2) - Н (высота насадки), х(3)- Gmg (массовый расход воздуха) %-------------------------------------------------------------------начальные условия t0=35; tH=30; %температура жидкости на входе и выходе из аппарата, 0С tgH=28; %температура воздуха на входе, 0С pH=2645.3; %парциальное давление пара на входе в аппарат, Па mp=-4062.9; np=276.1; %линейная аппроксимация парциального давления насыщенного пара %другие данные Gl=200; %объемный расход жидкости, м3/ч rol=998.2; %плотность жидкости, кг/м Gml=Gl*rol; %массовый расход жидкости, кг/ч qql=Gml/7200/x(1); %удельный расход жидкости, приходящийся на ед. ширины листа насадки qqg=x(3)/7200/x(1); %удельный расход газа - полурасход между листами насадки qul=4.5; %плотность орошения, кг/(м2 с) rog=1.165; %плотность газа, кг/м nug=16e-6; %кинематическая вязкость газа, м2/с Sv=Gml/3600/qul; %площадь горизонтального сечения насадки, м Vg=x(3)/3600/rog/Sv; %скорость газа, м/с b=Sv/x(1); %расстояние между листами, м Re=4*qqg/rog/nug; %число Рейнольдса %коэффициенты уравнений ТМО Nu=0.018*Re^0.8; %число Нуссельта a1=0.611e-5*Nu/def/qql;
b1=0.1013e-6*Nu/def/qql;
a2=0.2476e-4*Nu/def/qqg;
b2=0.2794e-4*Nu/def/qqg;
n1=(a1+b1*np-a2-b2)/2; k=sqrt(a2*b2-a1*b2-a2*b1*np);
lamb1=-n1+sqrt(n1^2-k^2);
lamb2=-n1-sqrt(n1^2-k^2);
d1=tgH+t0;d2=pH-mp-np*t0;
r1=a2*exp(lamb1*x(2))/(a2-lamb1)+1;r2=a2*exp(lamb2*x(2))/(a2-lamb2)+1;
s1=b2*np*exp(lamb1*x(2))/(b2-lamb1)-np;s2=b2*np*exp(lamb2*x(2))/(b2-lamb2)-np;
del=r1*s2-r2*s1;
del1=d1*s2-d2*r2;
del2=d2*r1-d1*s1;
C1=del1/del; C2=del2/del;
%нелинейное ограничение при заданной температуре воды на выходе из аппарата х=Н c=[];
ceq=tH-t0-C1*(exp(lamb1*x(2))-1)-C2*(exp(lamb2*x(2))-1);
определение оптимальных параметров:
x0=[80,0.1,30000];
options=optimset('LargeScale','off');
[x,fval,exitflag,output]=fmincon('zfun1',x0,[],[],[],[],[180,0.4,200000],[],'confuneq2') Получены следующие значения эксплуатационной ( Gг ) и конструктивных ( n L, H ) переменных (таблица 1), соответствующих минимуму капиталовложений в аппарат и вентиляторы, отнесённых к году окупаемости. При оптимизации функции П кроме ( 2 ) накладывалось также ограничение H 4 м..
5. Задачи линейного и квадратичного программирования [ 11 ]:
4.3. Simulink – система визуального моделирования Назначение пакета Simulink [ 12 ] Пакет Simulink служит для блочного моделирования различных систем и устройств. Для построения функциональной блок-схемы моделируемых устройств Simulink имеет обширную библиотеку блочных компонентов и удобный редактор блок-схем. Последний основан на использовании возможностей графического интерфейса пользователя и по существу является типичным средством визуального программирования. Используя палитры компонентов (наборы) блок-схем, пользователь с помощью мышки переносит нужные компоненты с палитр на рабочий стол пакета Simulink и соединяет линиями входы и выходы блоков. Таким образом, готовится блок-схема моделирования нужной системы или устройства.
Но этим уникальность пакета не исчерпывается, так как Simulink автоматизирует следующий наиболее трудоемкий этап моделирования: он автоматически составляет и решает сложные системы алгебраических и дифференциальных уравнений, описывающих заданную функциональную схему (модель), обеспечивая удобный и наглядный визуальный контроль за поведением созданного пользователем виртуального устройства. Вам достаточно уточнить (если нужно) вид анализа и пустить Simulink в режим симуляции (откуда и название пакета – Simulink) созданной модели системы или устройства. Изучение работы устройств и систем по моделям принято называть моделированием. Средства визуализации результатов моделирования в пакете Simulink настолько наглядны, что порою создаётся ощущение, что созданная в виде блок-схемы модель и впрямь работает как "живая". Более того, Simulink почти мгновенно меняет математическое описание модели по мере ввода ее новых блоков даже в том случае, когда это изменяет порядок системы уравнений и ведет к существенному качественному изменению поведения системы. Впрочем, именно это и является одной из важнейших целей применения пакета Simulink.
Ценность Simulink заключается и в обширной открытой для изучения и модификации библиотеке компонентов. Она включает в себя источники сигналов с практически любыми временными зависимостями (от скачка до произвольной функции, задаваемой пользователем), масштабирующие, линейные и нелинейные преобразователи с разнообразными формами передаточных характеристик, квантующее устройство, интегрирующие и дифференцирующие блоки и так далее.
В библиотеке имеется и приличный набор регистрирующих устройств – от простых измерителей типа вольтметра или амперметра до универсальных осциллоскопов, позволяющих с удобством просматривать временные зависимости выходных параметров моделируемых систем, например, токов и напряжений, перемещений, давления и тому подобное. Есть даже графопостроитель для построения фигур в полярной системе координат, например, фигур Лиссажу или фазовых портретов колебаний. Simulink имеет обширные средства анимации и даже звукового сопровождения. А в дополнительных библиотеках можно отыскать и такие "дорогие приборы", как анализаторы спектра сложных сигналов, многоканальные самописцы и средства анимации графиков.
Интеграция одной из самых быстрых математических систем MATLAB с удобным и мощным пакетом Simulink открывает новые возможности использования самых современных математических методов для решения задач динамического и ситуационного моделирования сложных систем и устройств. А средства анимации Simulink позволяют строить на его основе виртуальные физические лаборатории с очень наглядным представлением результатов моделирования. Возможности Simulink позволяют выполнять математическое моделирование сложных динамических систем в физике, электро- и радиотехнике, в биологии и химии, словом, во всех областях науки и техники. Это объясняет популярность данного пакета как в крупных университетах и институтах, так и научных лабораториях всего мира.
И, наконец, важным достоинством пакета является возможность задания в блоках произвольных математических выражений, что позволяет подчас решать типовые задачи, пользуясь примерами пакета Simulink и просто задавая новые выражения, описывающие работу моделируемых пользователем систем и устройств. Полезным свойством пакета является и возможность задания системных S-функций и полноценное их включение в состав библиотек Simulink.
К этому надо добавить и возможность моделирования устройств и систем в реальном масштабе времени.
Как программное средство Simulink – типичный представитель визуально-ориентированного языка программирования. На всех этапах работы, особенно при подготовке схем моделей, пользователь вообще не имеет дела с программированием как таковым. Нужная программа автоматически генерируется по мере ввода выбранных блоков компонентов, их соединений друг с другом и по мере задания параметров компонентов.
Важное достоинство Simulink – интеграция как с самой системой MATLAB, так и рядом других пакетов расширения, что обеспечивает по существу неограниченные возможности в применении Simulink для решения практически любых задач имитационного моделирования, причем как простых, так и самых сложных. Совместно с Simulink используется ряд пакетов. Ниже приведены краткие сведения о наиболее важных из них:
Communications – пакет моделирования коммуникационных устройств (более 100 функций и 150 компонентов);
Signal Processing – пакет для проектирования сложных устройств обработки сигналов;
Real Time Workshop – пакет для моделирования в реальном масштабе времени, обеспечивающий создание кода на языке программирования С;
Stateflow – пакет для моделирования событийно-управляемых систем на базе теории конечных автоматов;
DSP BIockset – библиотека моделирования цифровых сигнальных процессоров DSP (Digital Signal Processor) - цифровых фильтров.
Интеграция пакета Simulink с системой MATLAB После инсталляции Simulink (отдельно от MATLAB или в составе этой системы) он интегрируется с ней и не требует никаких дополнительных операций для работы. Внешне интеграция проявляется появлением кнопки New Simulink Model в конце панели инструментов (перед кнопкой справочной системы ?). Нажатие этой кнопки открывает окна редактора функциональных блок-схем и оглавления библиотеки компонентов - рисунке 4.11.
Рис. 4.11. Интеграция средств пакета Simulink с системой MATLAB Для решения автоматически составленной системы нелинейных дифференциальных уравнений первого порядка (ODE) Simulink использует решатель дифференциальных уравнений, построенный в виде программного цифрового интегратора. Решатель работает в двух основных режимах:
Variable-step solvers – решение с переменным шагом. Fixed-step solvers - решение с фиксированным шагом. Как правило, лучшие результаты дает решение с переменным шагом (обычно по времени, но не обязательно). В этом случае шаг автоматически уменьшается, если скорость изменения результатов в процессе решения возрастает. И, напротив, если результаты меняются слабо, шаг решения автоматически увеличивается. Это исключает (опять-таки, как правило) накапливание ошибки, которое нередко случается при фиксированном шаге.
Параметры решателя устанавливаются с помощью окна, которое появляется при исполнении команды Parameters в позиции Simulation главного меню пакета Simulink. В нем же можно установить конкретный метод решения дифференциальных уравнений: ode45, ode23, rk45 (метод Рунге-Кутта), odel13 (метод Адамса), ode15s и odel (метод Эйлера). Своими возможностями пакет Simulink во многом обязан специальному аппарату создания и применения так называемых системных S-функций (System Functions). Эти функции позволяют в ходе решения осуществлять сложные функциональные преобразования по различным математическим алгоритмам, например алгоритмам решения систем дифференциальных уравнений. Для подготовки S-функций Simulink имеет специальный редактор. Создав свою S-функцию, пользователь фактически создает блок библиотеки, который может использоваться по всем правилам применения блоков. Это значит, что блок можно переносить в окно редактирования мышкой, менять заданную S-функцию, использовать необходимые связи между блоками и так далее.
Особенности интерфейса Simulink Интерфейс новой версии Simulink 4 полностью соответствует стилю интерфейса типовых приложений под Windows 95/98, в том числе интерфейсу системы MATLAB. В тоже время он концептуально строг, дабы не досаждать пользователю многочисленными "излишками" стандартного интерфейса Windows 95/98. Главное меню системы имеет следующие позиции:
File – работа с файлами моделей и библиотек (их создание, сохранение, считывание и печать).
Edit – операции редактирования, работа с буфером промежуточного хранения и создание субсистем.
View – вывод или удаление панелей инструментов и статуса.
Simulation – управление процессом моделирования (старт, ввод паузы и вывод окна настройки параметров симуляции).
Format – операции форматирования модели (смена фонтов, редактирование надписей, повороты блоков, использование тени от блоков операции с цветами линий блоков, их фоном и общим фоном).
Tools – управление видом анализа (в линейной области и в режиме реального времени RTW).
Первые три позиции главного меню содержат общепринятые для Windows-приложений команды и операции, так что мы не будем их обсуждать.
Остальные позиции будут рассмотрены по мере знакомства с системой Simulink.
Как уже отмечалось, вместе с рабочим окном Simulink выводится окно с перечнем разделов основной библиотеки компонентов. Это окно – важная специфическая часть интерфейса Simulink. Оно открывает доступ к множеству других подобных по оформлению окон, дающих доступ к новым обширным пакетам компонентов (Blocksets&Toolboxes) и примеров их применения (Demos).
Это позволяет пользователю постепенно знакомиться все с новыми и новыми областями применения Simulink.
Демонстрация возможностей Simulink Даже не знакомый с Simulink пользователь может быстро оценить возможности этого пакета, воспользовавшись множеством интересных и поучительных примеров его применения. Они находятся в довольно далекой директории MATLAB/TOOLBOX/SIMULINK/SIMDEMO10.
Для загрузки примеров используется (как обычно) команда Open в позиции File главного меню основного окна пакета Simulink. В уже знакомом вам окне, правда, с названием Open Model (Открыть Модель), надо войти в указанную директорию с примерами и выбрать подходящий пример.
Выбрав нужный пример, нужно нажать кнопку Открыть окна. Появится графическая модель моделируемого устройства. На рисунке 2.3 показана выбранная нами модель, позволяющая проанализировать поведение массивного кубика, который через пружину привязан к палочке (файл и пример onecart).
«