Метод конечных разностей

АНАЛИТИЧЕСКОЕ РЕШЕНИЕ

Компьютерные программы на языке МАТLAB

% Балка на упругом основании. 11.08.14

% Аналитическое решение.

% Балка из двутавра 18 при равномерной нагрузке q.

% Y1, Y2, Y3, Y4 функции Крылова – вписаны.

% Определяет прогибы – v, углы поворота – f,

% изгибающие моменты – M, поперечные силы – Q,

% реактивный отпор – qr.

% ИСХОДНЫЕ ДАННЫЕ:

% E – модуль упругости стали, b – ширина полки двутавра

% J – осевой момент инерции, k – коэффициент постели,

% v0, fi0 – начальные параметры, F – сосредоточенная сила,

% М – моментная нагрузка, q – равномерно распределённая нагрузка,

% l – длина балки,

% z – продольная координата, ksi – безразмерная координата,

clear; disp('___________________');

E=210e+9; J=1290e-8; l=5,2; c=0.9; n=601; k0=60e+6; b=E*J;

k=k0*c; h=l/(n-1); z=0:h:l; la=sqrt(sqrt(k/4/b)); ksi=la*z;

Y1=zeros(1,n);

F=–30000; M=0; q=–20000;

% Определение значений функций Крылова

Y1=cosh(ksi).*cos(ksi);

Y2=1/2*(cosh(ksi).*sin(ksi)+sinh(ksi).*cos(ksi));

Y3=1/2*sin(ksi).*sinh(ksi);

Y4=1/4*(cosh(ksi).*sin(ksi)-sinh(ksi).*cos(ksi));

% Bычисление начальных параметров

% y1=Y1(n); y2=Y2(n); y3=Y3(n); y4=Y4(n);

A(1,1)=Y1(n); A(1,2)=Y2(n)/la; A(2,1)=-4*la*Y4(n); A(2,2)=Y1(n);

B(1,1)=-1/b*(M/la^2*Y3(n)+F/la^3*Y4(n)+q/4/la^4*(1-Y1(n)));

B(2,1)= -1/b*(M/la*Y2(n)+F/la^2*Y3(n)+q/la^3*Y4(n));

vf0=A\B;

v0=vf0(1); fi0=vf0(2);

for i=1:n;

y1=Y1(i); y2=Y2(i); y3=Y3(i); y4=Y4(i);

v(i)=v0*y1+fi0/la*y2+1/b*(M/la^2*y3+F/la^3*y4+q/4/la^4*(1-y1));

fi(i)=-4*la*v0*y4+fi0*y1+1/b*(M/la*y2+F/la^2*y3+q/la^3*y4);

Mi(i)=-4*b*la^2*v0*y3-4*b*la*fi0*y4+M*y1+F/la*y2+q/la^2*y3;

Q(i)=-4*b*la^3*v0*y2-4*b*la^2*fi0*y3-4*la*M*y4+F*y1+q/la*y2;

qr(i)=-k*v(i);

end;

disp(trapz(z,qr)) % Вычисление интеграла методом трапеций

% В Ы В О Д Г Р А Ф И К О В Н А Э К Р А Н

plot(z, v*100,'k','LineWidth',1); grid; xlabel('z,м'); ylabel('v,см');

figure;

plot(z, fi,'k','LineWidth',1); grid; xlabel('z, м'); ylabel('fi,рад');

figure;

plot(z, Mi/1000,'k','LineWidth',1); grid; xlabel('z,м'); ylabel('M,кНм');

figure;

plot(z, Q/1000,'k','LineWidth',1); grid; xlabel('z,м'); ylabel('Q,кН');

figure;

plot(z, qr/1000,'k','LineWidth',1); grid; xlabel('z,м'); ylabel('qr,кН/м');

disp('___________________');

% Балка на упругом основании. 11.08.14

% Метод конечных разностей.

% Балка из двутавра 18 при равномерной нагрузке q.

% Используется уравнение 4-го порядка: v''''+m*v=q(z), m=k/EJ

% Граничные условия в данной программе:

% v''(0)=M, v'''(0)=F. v(l)=0, v'(l)=0;

% Задача приводится к алгебраической системе уравнений

% A v = p, (1)

% p=qh^4

% Система (1) решается с помощью подпрограммы

% вычислительного комплекса MATLAB

clear; disp('___________________');

% ИСХОДНЫЕ ДАННЫЕ:

l=6; E=210e+09; J=1290e-08; n=1001; с=0.9; k0=60e+06;

F=-30000; M=0; q=-20000;

b=E*J; q0=q/b;

k=с*k0; n1=n-1; h=l/n1; z=[0: h: l];

m=k/b; mu=6+m*h^4;

v=zeros(1,n); p=zeros(1,n); fi=zeros(1,n); A=zeros(n,n);

MO=zeros(1,n); Q=zeros(1,n);

% Метод конечных разностей

% Формирование матрицы A и вектора р

% Граничные условия

A(1,1)=2; A(1,2)=-5; A(1,3)=4; A(1,4)=-1;

A(2,1)=-5; A(2,2)=18; A(2,3)=-24; A(2,4)=14; A(2,5)=-3;

A(n1,n-2)=1; A(n1,n-1)=-4; A(n1,n)=3;

A(n,n)=1; p(1)=h^2*M/b; p(2)=2*h^3*F/b;

for i=3:n-2; p(i)=q0*h^4;

A(i,i-2)=1; A(i,i-1)=-4; A(i,i)=mu; A(i,i+1)=-4; A(i,i+2)=1;

end;

% Изогнутая ось и реактивный отпор основания

v=A\p'; qr=-k*v;

% Вычисления для эпюр

% По левому концу

fi(1)=(-3*v(1)+4*v(2)-v(3))/2/h; fi(2)=(-v(1)+v(3))/2/h;

MO(1)=(2*v(1)-5*v(2)+4*v(3)-v(4))/h^2;

MO(2)=(v(1)-2*v(2)+v(3))/h^2;

Q(1)=(-5*v(1)+18*v(2)-24*v(3)+14*v(4)-3*v(5))/2/h^3;

Q(2)=(-3*v(1)+10*v(2)-12*v(3)+6*v(4)-v(5))/2/h^3;

% По правому концу

fi(n)=(v(n-2)-4*v(n-1)+3*v(n))/2/h;

fi(n-1)=(-v(n-2)+v(n))/2/h;

MO(n)=(-v(n-3)+4*v(n-2)-5*v(n-1)+2*v(n))/h^2;

MO(n-1)=(v(n-2)-2*v(n-1)+v(n))/h^2;

Q(n)=(3*v(n-4)-14*v(n-3)+24*v(n-2)-18*v(n-1)+5*v(n))/2/h^3;

Q(n-1)=(v(n-4)-6*v(n-3)+12*v(n-2)-10*v(n-1)+3*v(n))/2/h^3;

% В регулярных точках

for i=3:n-2;

fi(i)=(v(i+1)-v(i-1))/2/h;

MO(i)=(v(i-1)-2*v(i)+v(i+1))/h^2;

Q(i)=(-v(i-2)+2*v(i-1)-2*v(i+1)+v(i+2))/2/h^3;

end;

% Изгибающий момент и поперечная сила в сечениях

M1=MO*b; Q1=Q*b;

% В Ы В О Д Г Р А Ф И К О В Н А Э К Р А Н

plot(z, v*100,'k','LineWidth',1); grid; xlabel('z,м'); ylabel('v,см');

figure;

plot(z, fi,'k','LineWidth',1); grid; xlabel('z, м'); ylabel('fi,рад');

figure;

plot(z, M1/1000,'k','LineWidth',1); grid; xlabel('z,м'); ylabel('M,кНм');

figure;

plot(z, Q1/1000,'k','LineWidth',1); grid; xlabel('z,м'); ylabel('Q,кН');

figure;

plot(z, qr/1000,'k','LineWidth',1); grid; xlabel('z,м'); ylabel('qr, кН/м');

disp('___________________');

Дата добавления: 2017-01-13