Пример текста М-функции

⇐ ПредыдущаяСтр 5 из 7Следующая ⇒

Пример текста М-функции

function beam_FEM_cub

% ввод,задание и формирование исходных данных

L=input('ввести длину балки L=');

hb=input('ввести высоту поперечного сечения hb=');

bb=input('ввести ширину поперечного сечения bb=');

Eb=input('ввести модуль упругости Eb=');

kb=input('ввести коэффициент отпора грунта kb=');

bt=kb*bb;

J=bb*hb^3/12;

EJ=Eb*J;

n=input('ввести число точек n=');

ne=n-1; %количество элементов

h=L/ne; % длина элемента

mp=2; %число неизвестных в узле: y(xi),y'(xi)

me=2*mp;%число неизвестных на элементе: y(x_i),y'(x_i),y(x_i+1),y'(x_i+1)

ns=mp*n; %общее количество неизвестных

Pf=input('ввести параметр нагрузки Pf=');

xi=(0:h:L)'; % координаты узлов

bnd_str=sprintf('%s\n','1 - y(0)=0,y"(0)=0;y(L)=0,y"(L)=0',...

'2 - y(0)=0,y''(0)=0;y(L)=0,y''(L)=0',...

'3 - y(0)=0,y"(0)=0;y"(L)=0,y"''(L)=0',...

'4 - y(0)=0,y''(0)=0;y"(L)=0,y"''(L)=0');

fprintf('\nввести номер граничных условий bnd_cond:\n%s\n',bnd_str)

bnd_cond=input('bnd_cond=');

% формирование вспомогательных матриц и векторов

Ag=[1 0 0 0;0 1 0 0;1 1 1 1;0 1 2 3]

A=inv(Ag)

S=diag([1 2 3]',-1)

T0=[ 1 1/2 1/3 1/4

1/2 1/3 1/4 1/5

1/3 1/4 1/5 1/6

1/4 1/5 1/6 1/7]

TT=420*T0

Di=diag([1 h 1 h]')

t0=[1 1/2 1/3 1/4]'

% формирование локальных матриц и векторов

K0=h*Di*A'*T0*A*Di

K2=1/h^3*Di*A'*S^2*T0*(S^2)'*A*Di

Ke=EJ*K2+bt*K0 %локальная матрица жесткости

R0=h*Di*A'*t0 %локальный вектор нагрузки

% формирование глобальной матрицы жесткости

Kg=zeros(ns);

for i=1:mp:ns-mp-1

Kg(i:i+me-1,i:i+me-1)=Kg(i:i+me-1,i:i+me-1)+Ke(:,:);

end

% формирование глобального вектора нагрузки

Rg=zeros(ns,1);

Rg(ns/2)=Pf;

% учет граничных условий

if bnd_cond==1

Kg(1,:)=0;Kg(:,1)=0;Kg(ns-1,:)=0;Kg(:,ns-1)=0;

Kg(1,1)=1;Kg(ns-1,ns-1)=1;

elseif bnd_cond==2

Kg(1:2,:)=0;Kg(:,1:2)=0;Kg(ns-1:ns,:)=0;Kg(:,ns-1:ns)=0;

Kg(1,1)=1;Kg(2,2)=1;Kg(ns-1,ns-1)=1;Kg(ns,ns)=1;

elseif bnd_cond==3

Kg(1,:)=0;Kg(:,1)=0;

Kg(1,1)=1;

else bnd_cond==4

Kg(1:2,:)=0;Kg(:,1:2)=0;

Kg(1,1)=1;Kg(2,2)=1;

end

%Kg

%Rg

U=Kg\Rg;

yi=U(1:mp:ns);

dyi=U(2:mp:ns);

for i=1:ne

xs(i)=h*(i-0.5);

d2y(i)=df(i,2,0.5);

d3y(i)=df(i,3,0.5);

end

figure()

plot(xi,yi),grid on,title('progib Y')

figure()

plot(xi,dyi),grid on,title('ugol dY')

figure()

plot(xs,-EJ*d2y),grid on,title('moment -EJ*d2Y')

figure()

plot(xs,-EJ*d3y),grid on,title('poperechnaja sila -EJ*d3Y')

function dpyi=df(i,p,t)

% вычисление dpy в точке t на i-м элементе

Dp=(S^p)'*A*Di;

ti=[1;t;t^2;t^3];

in=1+mp*(i-1);ik=in+me-1; z=U(in:ik);

dpyi=dot(Dp*z,ti)/h^p;

end

1. Краевые условия – балка, шарнирно опертая с двух сторон

⇐ Предыдущая 1 2 3 456 7 Следующая ⇒