Алгоритм нахождения точки бифуркации

Чтобы найти точку бифуркации, используется идея метода деидеализации [9]. Для данной задачи метод выглядит следующим образом (рис. 6):

­

а

б

в

г

Рис. 6

когда амортизатор находится в ненагруженном состоянии, к нему прикладывается некая нагрузка, обеспечивающая несимметричное состояние, в данном случае прикладываемое усилие или фиксированное смещение по оси (рис. 6а). Далее, амортизатор сжимается нагрузкой до какого-то заданного положения (рис. 6б). Потом, боковая нагрузка убирается (рис. 6в), при этом амортизатор остается в несимметричном положение. После чего, уменьшается нагрузка и амортизатор возвращается в исходное положение (рис. 6г).

Решается задача для систем (1.4) и (1.5) с граничными условиями (1.6). В качестве основного параметра выбран параметр дискретного перемещения верхней пластины по оси , обозначим его .

Расчет проводился для следующей параметров амортизатора: длина верхней пластины – , толщина боковых пластин – , длина боковых пластин – .

Так как в данной работе ищутся точки, в которых амортизатор при вертикальной нагрузке без смещений по другим осям будет склонен принять несимметричную форму, логично будет рассматривать углы, близкие к .

Бифуркация

Рис. 8. Амортизатор в плоскости , несимметричное сжатие

Рис. 7. Амортизатор в плоскости , симметричное сжатие

Рассмотрим сжатие амортизатора с начальным углом наклона боковых пластин к оси равному (рис. 7–8). На всех приведенных ниже рисунках в качестве боковых пластин проиллюстрированы их срединные линии.

На рисунке 8 несимметричная деформация реализована следующим путем: по оси к амортизатору приложено усилие .

P

Далее, амортизатор сжимается дискретно по параметру до момента, когда угол между касательной в точке скрепления пластин и верхней пластиной не становится равным нулю. Потом, используя метод деидеализации, изложенный выше, приравнивается к нулю и сжимающая нагрузка убирается до момента, пока амортизатор не вернется в состояние покоя. На всех диаграммах, приведенных ниже, проиллюстрирован процесс возврата амортизатора в исходное положение, то есть, движение по параметру осуществляется справа налево.

Рис. 9. Диаграмма нагрузка-перемещение

Из рисука 9 видно, что в точке графики несимметричной и симметричной форм сливаются в один. Стоит подметить, что несимметричный вариант является выгоднее симметричного, потому что для его реализации нужна меньшая вертикальная нагрузка. Рассмотрим окрестность точки ветвления ближе (рис. 10):

Рис. 10. Диаграмма нагрузка-перемещение

P

На промежутке (рис. 10) наблюдается численный переско

Рис. 10. Диаграмма нагрузка-перемещение

к. Если в окрестности этой точки шаг параметра уменьшить в 100 раз, то наблюдается то же самое, что и на рисунке 10. При дальнейшем уменьшении шага наблюдается ухудшение сходимости. Следовательно, требуется пройти по другому параметру.

Рис. 17. Диаграмма нагрузка-перемещение

При движении по параметру – параметру вертикального нагружения (рис. 11) несимметричная ветвь снова перескакивает на симметричную.

Рис. 11. Диаграмма нагрузка-перемещение

P

P

Рассмотрим точку ветвления ближе (рис. 12):

Рис. 12. Диаграмма нагрузка-перемещение

Рис. 13. Диаграмма нагрузка-перемещение

Из рисунков 11–12 видно, что графики совпадают, учитывая погрешность вычисления. Чтобы получить больше данных для полноценного анализа, требуется пройти по еще одному параметру.

Рассмотрим движение по параметру – параметру угла поворота (рис. 13). Также, на рисунке 13 симметричное решение продолжено по до конца высоты амортизатора для дальнейших исследований.

Рассмотрим окрестность точки соединения графиков ближе (рис. 14):

Рис. 14. Диаграмма нагрузка-перемещение

Рис. 15. Диаграмма зависимости от трех параметров продолжения решения

P

P

Каким бы мелким ни был шаг, при движении по наблюдается остановка в окрестности этой точки и «разворот» назад по той же кривой в этой плоскости (рис. 14). То есть, следует рассмотреть трехмерный график зависимости всех от трех параметров (рис. 15).

Из рисунков 9, 11 и 13 видно, что при движении по любому параметру не осуществляется плавного перехода на симметричное решение, точка является предельной для всех параметров. Также, если в окрестности этой точки пройти по симметричному варианту с мелким шагом, то наблюдается сильно ухудшение сходимости вплоть до остановки алгоритма. Следовательно, эта точка – точка бифуркации, то есть точка разветвления решения. В ней при нагружении будет происходить переход симметричного варианта на несимметричный.

Рис. 16. Бифуркационная диаграмма нагрузка-поворот

P

Рассмотрим диаграмму (рис. 15) в плоскости (рис. 16).

Зафиксируем и в точках пересечения с диаграммой (рис. 16) построим формы амортизатора (рис. 17–18).

Точка С 1)

Точка С 2)

Рис. 17. Формы симметричного варианта

Точка Н 1)

Точка Н 2)

Рис. 18. Формы несимметричного варианта

Также, было выяснено, что при данных расчетных параметрах амортизатора бифуркация имеет место быть при . То есть, при будет реализовываться симметричный вариант.

Изолированное решение

Помимо бифуркации, методом деидеализациии было найдено изолированное решение.

Рис. 19. Амортизатор в плоскости , симметричное сжатие

Рис. 20. Амортизатор в плоскости , несимметричное сжатие

Рассмотрим сжатие амортизатора с начальным углом наклона боковых пластин к оси , равному (рис. 19–20).

Рис. 21. Диаграмма нагрузка-перемещение

P

На рисунке 20 продемонстрирована несимметричная деформация, полученная немного иным путем, чем в главе 3.2: верхней пластине дается фиксированное смещение по оси , обозначим его . Далее все происходит аналогично главе 3.2.

Из рисунка 21 видно, что в точке графики нагрузки-перемещения для симметричной и для несимметричной форм сливаются в один. Рассмотрим окрестность этой точки ближе (рис. 22).

Рис. 22. Диаграмма нагрузка-перемещение

P

На рисунке 22 шаг параметра равен . Видно, что нагрузка убывает, потом снова возрастает и на промежутке осуществляется численный перескок с несимметричной ветви на симметричную. В этом промежутке сходимость ухудшается, поэтому требуется поменять параметр продолжения решения.

Если в окрестности точки перескока пройти по параметру – параметру угла поворота, то получается изолированное решение, изображенное на рисунке 23. Реализовываться этот вариант, скорее всего, не будет, тем более, симметричный вариант в этом случае выгоднее несимметричного по нагрузке.

Рис. 23. Диаграмма нагрузка-перемещение

P

Выводы

Составлен алгоритм решения нелинейной краевой задачи методом стрельбы в сочетании с методом продолжения по параметру. Написана программа, реализующая этот алгоритм. Отладка программы проводилась на примерах, известных из литературы.

Была рассмотрена задача деформации мостичного амортизатора при условии шарнирного опирания у верхнего основания и заделки у нижнего основания. Также, предполагалось, что основания при сжатии остаются параллельными, но верхнее основание может перемещаться по горизонтали. Задача в такой постановке другими авторами не рассматривалась.

Найдены бифуркационные ветви перехода от симметричного деформирования к несимметричному. Показано, при каком значении начального угла будет возможна бифуркация при расчетных параметрах.

В развитие работы можно рассмотреть другие толщины резиновых пластин; поставить условие возможности поворота верхнего основания при шарнирном опирании; рассмотреть криволинейные формы резиновых пластин.

Заключение

В ходе выполнения была написана программа, реализующая решение поставленной задачи. Найдены ветви потери устойчивости амортизатора и диапазон угла, при котором она реализуется для заданных расчетных параметров.

Список литературы

1. Бочкарева Н.Л., Колпак Е.П. Об устойчивости арочного амортизатора // Вестник СПбГУ, 1993. Сер. 1, вып. 4. С. 49–53.

2. Колпак Е.П. Устойчивость мостичного амортизатора из резиноподбного материала // Теория и методы расчета нелинейных пластин и оболочек, 1981. С.: Изд-во Саратовского ун-та, 1981. С. 43–45.

3. Черных К.Ф. Нелинейная теория упругости в машиностроительных расчетах. Л.: Машиностроение, 1986. 336 с.

4. Кабриц С.А., Михайловский Е.И., Товстик П.Е., Черных К.Ф., Шамина В.А. Общая нелинейная теория упругих оболочек. СПб.: СПбГУ, 2002. 386 с.

5. Kabrits S.A., Kolpak E.P. Finding bifurcation branches in nonlinear problems of statics of shells numerically, 2015 International conference «Stability and Control processes» in Memory of Zubov V.I., 2015. P. 389–391.

6. Kabrits S.A., Kolpak E.P. Numerical Study of Convergence of Nonlinear Models of the Theory of Shells with Thickness Decrease, AIP Conference Proceedings, 2016. Vol. 1738, №160006.

7. Бахвалов Н.С. Численные методы (анализ, алгебра, обыкновенные дифференциальные уравнения). М.: Наука, 1975. 631 с.

8. Кабриц С.А., Терентьев В.Ф. О численном построении диаграмм нагрузка-перемещение в одномерных нелинейных задачах теории стержней и оболочек // Актуальные проблемы нелинейной механики сплошных сред, 1977. Вып. 1. С. 155–171.

9. Алфутов Н.А. Основы расчета на устойчивость упругих систем. М.: Машиностроение, 1978. 312 с.

10. Пановко Я.Г., Губанова И.И. Устойчивость и колебания упругих систем. М.: Наука, 1987. 352 с.

11. Белкин А.Е., Семенов В.В., Семенов В.К. Численный анализ больших плоский деформаций арочного амортизатора // Вестник МГТУ им. Н.Э. Баумана, 2011. №2. С. 55–64.

12. Белкин А.Е., Хоминич Д.С. Расчет больших деформаций арочного амортизатора с учетом объемной сжимаемости резины // Вестник МГТУ им. Н.Э. Баумана, 2012. №2. С. 3–11.

Приложение

Краткое описание всех функций:

main –скрипт, запускающий работу программы;

BALKAM – функция, осуществляющая продолжения по параметру (сжатие амортизатора, несимметричный случай);

BALKAM1 – функция, осуществляющая продолжения по параметру (сжатие амортизатора, несимметричный случай);

PSI1 / PSI2 – функция от вектора начального приближения, задающая начальные условия для задачи Коши для правой / левой половинки (1.3);

PSI / PSITh – функция от вектора начального приближения, задающая граничные условия на верхнем основании при параметре продолжения решения вертикального перемещения / угла поворота (1.3);

FRESHE – функция от функции PSI1/PSI2и вектора начального приближения для вычисления невязки и матрицы Якоби в методе Ньютона;

PSI_PRINT1 / PSI_PRINT2 –функция от вектора деформации, осуществляющая построение формы правой / левой половинки амортизатора;

Runge –модуль метода Рунге – Кутты – Мерсона для системы дифференциальных уравнений (1.1);

solve_lambda –функция, вычисляющая методом касательных кратность удлинения срединной линии резиновых пластин амортизатора;

RVO –функция задания начальной геометрии и констант материала;

kgx –функция от вектора начального приближения, вычисляющая правые части системы дифференциальных уравнений (1.1).

function [g] = kgx(s,y)

%вычисление правых частей системы дифференциальных уравнений

% y[1]=v y[2]=u y[3]=Theta y[4]=M y[5]=Tx y[6]=Tz L1=lambda1

global lg

lg = 1;

[ R0,Z0,F0,K0,h0,mu,n,beta] = RVO(s);

fi=y(3)+F0; g=zeros(6,1);

Tn=y(5).*sin(fi)+y(6).*cos(fi);

Ts=y(5).*cos(fi)-y(6).*sin(fi);

lambda10=1;

lambda2=lg;

alfa=mu.*h0./n.*((1+beta)+(1-beta).*lambda2.^n);

L1=solve_lambda(Ts,lambda10,alfa,n,lambda2);

lambda3=1./(L1.*lambda2);

g(1)= L1.*cos(fi)-cos(F0);

g(2)=-L1.*sin(fi)+sin(F0);

b=(L1*lambda2)^n;

a=1/(L1.^2.*lambda2).*((1+beta)./b+(1-beta).*b);

kappa_s=6.*y(4)./(mu.*h0.^3)./a;

g(3)=L1.^2.*lambda2.*(kappa_s+K0)-K0;

g(4)=L1*Tn;

g(5)=0;

g(6)=0;

end

function [R0,Z0,F0,K0,h0,mu,n,beta] = RVO(s)

%начальная геометрия и константы материала

global ab

F0=85/180*pi;

R1=0.5; R0=R1+s.*cos(F0);

H=ab*sin(F0); Z0=H-s.*sin(F0);

K0=0;

h0=1/7;

mu=1;

n=2;

beta=1;

end

function [L1] = solve_lambda(T,L10,alfa,n,L2)

%вычисление кратности удлинения срединной линии резиновых пластин методом касательных

global lg

L1=L10; eps=1; b=L2.^(-n);

while eps>1E-10

a=L1.^(n-1); c=b./(a*L1*L1);

f=alfa.*(a-c)-T;

df=alfa.*((n-1)*a+(n+1).*c)/L1;

L11=L1-f/df;

eps=abs(L11-L1);

L1=L11;

lg = L11;

end

end

function [ Y ] = Runge(FT,TT0,TTK,TTH,EET,XX)

%Метод Рунге - Кутты - Мерсона с автоматическим шагом

T=TT0; C=abs(TTK-TT0)/3; NN=size(XX,1);

if C<=10^(-8)*abs(TTH)

return

end;

TTH=abs(TTH).*sign(TTK-TT0);

ib=1;

while (ib==1)

if (T+TTH-TTK)*sign(TTH)>0

TTH=TTK-T; ib=0;

end;

F0=FT(T,XX);

Y1=XX+F0.*TTH./3;

F1=FT(T+TTH./3,Y1); Y1=XX+(F0+F1).*TTH./6;

F2=FT(T+TTH/3,Y1); Y1=XX+(F0+3.*F2).*TTH./8;

F3=FT(T+TTH/2,Y1); Y1=XX+(F0-3.*F2+4.*F3).*TTH./2;

F1=FT(T+TTH,Y1);

Y1=abs(F0-4.5.*F2+4.*F3-0.5.*F1).*C;

II=0; JJ=0;

for I=1:NN

if EET(I)~=0

II=II+1;

if Y1(I) >= EET(I)*5

TTH=TTH./2;

continue;

else

if Y1(I)<= EET(I)/32

JJ=JJ+1;

end

end

end

end

T=T+TTH;

XX=XX+(F0+4.*F3+F1).*TTH./6;

if (II~=0) && (II==JJ)

TTH=2.*TTH;

end

if abs(T-TTK)>1.E-6.*C

continue

end

end

Y=XX;

end

function [Y] = PSI_PRINT1(X)

%построение формы правой пластинки амортизатора

global ab shag dlt ln

Z=zeros(6,1);

[R0,Z0,F0,K0,h0,mu,n,beta] = RVO(0);

Z(1)=0; %начальные условия

Z(2)=0;

Z(3)=0;

Z(4)=X(1);

Z(5)=X(2);

Z(6)=X(3);

Ep=zeros(6,1);

m=50;

YY=zeros(m,6);

x=linspace(ab,0,m);

for i=1:m-1

[R0,Z0,F0,K0,h0,mu,n,beta] = RVO(x(i));

if i==1

YY(i,1)=R0+Z(1);

YY(i,2)=Z0+Z(2);

end

ZZ=Runge(@kgx,x(i),x(i+1),shag,Ep,Z);

Z=ZZ;

[R0,Z0,F0,K0,h0,mu,n,beta] = RVO(x(i+1));

YY(i+1,1)=R0+Z(1);

YY(i+1,2)=Z0+Z(2);

end

ln(1) = R0+Z(1);

ln(2) = Z0+Z(2);

hold on; grid on;

if (dlt<0)

subplot(2,2,1);

plot(YY(:,1),YY(:,2),'LineWidth',4);

else

subplot(2,2,2);

plot(YY(:,1),YY(:,2));

end

hold on; grid on;

end

function [Y] = PSI_PRINT2(X)

%построение формы левой пластинки амортизатора

global ab shag dlt ln

Z=zeros(6,1);

[ R0,Z0,F0,K0,h0,mu,n,beta] = RVO(0);

Z(1)=0; %начальные условия

Z(2)=0;

Z(3)=0;

Z(4)=X(4);

Z(5)=X(5);

Z(6)=X(6);

Ep=zeros(6,1);

m=50;

YY=zeros(m,6);

DD1=zeros(m,6);

DD2=zeros(m,6);

YY1=zeros(m,6);

x=linspace(ab,0,m);

for i=1:m-1

[R0,Z0,F0,K0,h0,mu,n,beta] = RVO(x(i));

if i==1

YY(i,1)=-R0-Z(1);

YY(i,2)=Z0+Z(2);

end

ZZ=Runge(@kgx,x(i),x(i+1),shag,Ep,Z);

Z=ZZ;

[R0,Z0,F0,K0,h0,mu,n,beta] = RVO(x(i+1));

YY(i+1,1)=-R0-Z(1);

YY(i+1,2)=Z0+Z(2);

end

ln(3) = -R0-Z(1);

ln(4) = Z0+Z(2);

hold on; grid on;

if (dlt<0)

subplot(2,2,1);

plot(YY(:,1),YY(:,2),'LineWidth',4);

else

subplot(2,2,2);

plot(YY(:,1),YY(:,2));

end

hold on; grid on;

end

function [Y, DF] = FRESHE(PSI,X);

%вычисление Y=F(x) (невязка) и DF (матрица якоби dPSI/dx)

nfresh=size(X,1);

Y=PSI(X); %вектор деформации

for k=1:nfresh

A=X(k); P=A*1.E-4;

if abs(P)<1.E-8

P=1.E-8;

end

X(k)=A+P;

Z=PSI(X);

X(k)=A;

for j=1:nfresh

DF(j,k)=(Z(j)-Y(j))./P;

end

end

end

function [ZZ] = PSI1(X)

%начальные условия для решения задачи коши правой половинки

global ab shag

Z=zeros(6,1);

Z(1)=0;

Z(2)=0;

Z(3)=0;

Z(4)=X(1);

Z(5)=X(2);

Z(6)=X(3);

Ep=zeros(6,1);

ZZ=Runge(@kgx,ab,0,shag,Ep,Z); %решение задачи коши

end

function [ZZ] = PSI2(X)

%начальные условия для решения задачи коши левой половинки

global ab shag

Z=zeros(6,1);

Z(1)=0;

Z(2)=0;

Z(3)=0;

Z(4)=X(4);

Z(5)=X(5);

Z(6)=X(6);

Ep=zeros(6,1);

ZZ=Runge(@kgx,ab,0,shag,Ep,Z); %решение задачи коши

end

function [Y] = PSI(X)

% граничные условия на верхнем основании при параметре продолжения вертикального перемещения

global delta ab shag deltax

Z1 = PSI1(X); %правая половинка

Z2 = PSI2(X); % левая половинка

Y(1) = Z1(4); %граничные условия

Y(2) = Z2(4);

Y(3) = Z1(1) + Z2(1);

Y(4) = Z1(2) - Z2(2);

Y(5) = Z2(5) - Z1(5) + deltax;

Y(6) = Z1(2) - delta;

end

function [Y] = PSITh(X)

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

global tzz

Z1 = PSI1(X);

Z2 = PSI2(X);

Y(1) = Z1(4); Y(2) = Z2(4);

Y(3) = Z1(1) + Z2(1);

Y(4) = Z1(2) - Z2(2);

Y(5) = Z2(5) - Z1(5);

Y(6) = Z2(3) - tzz;

end

function [] = BALKAM

%функция, осуществляющая продолжения по параметру (сжатие амортизатора, несимметричный случай)

global q1 obrshag mas lnln p1 tzz distance TzL TzR grdlt delta ab dlt shag deltax ln t ogr

ln = zeros(4);

lnln = [];

[R0,Z0,F0,K0,h0,mu,n,beta] = RVO(0);

nfresh=3; %размерность

ab=1; %длина боковых пластин

shag=0.05; %шаг интегрирования

kol = 6; %размерность вектора системы дифференциальных уравнений

X=zeros(kol,1);

Xp=zeros(kol,1);

X0=zeros(kol,1);

X1=zeros(kol,1);

q1 = 1;

mas = []; %массив данных

TzL = [];

TzR = [];

Tz = [];

grdlt = []; %вектор параметров

rr = 0;

er = 0;

ch = 0; %счетчик

while (abs(delta)<(ab*sin(F0)-abs(dlt))) && (ogr-abs(dlt)>abs(delta))&& (q1~=0); %цикл сжатия

iter=0; %счетчик

s=1; %параметр сходимости

Xp = X;

rrr = 0;

if (er~=1)

delta = delta+dlt; %продолжение по параметру

end

if (delta<distance) %проверка, дошел ли процесс до заданного параметра и последующий разворот назад

deltax = 0;

dlt = -dlt/obrshag;

rrr = 1;

rr = 1;

end;

if (rr == 1) && (delta>-p1-0.0000001) && (delta<-p1+0.0000001) %условие смена параметра

er = 1;

end;

if (delta>-p1-0.0000001) && (er == 1) && (delta<-p1+0.0000001) %смена параметра

ll = length(Tz);

thxR = (Tz(ll)-Tz(ll-1))/100;

thxRR = thxR;

thetaxR = Tz(ll);

end;

if (er == 1) && (rr == 1) && (delta>-0.007) %обратная смена параметра при успешном проходе

er = 0;

end;

if (delta>-dlt+dlt/10) && (dlt>0) %условие остановки алгоритма

q1 = 0;

end;

if (er == 1) && (delta<-0.007) %экстраполяция параметра

tzz = thetaxR + thxR;

thxR = thxR +thxRR;

end;

while s>0.0001 %цикл сжатия

iter=iter+1;

if (er == 1) && (delta<-0.007) %определение параметра

if s>0.1 %условие сходимости

[Y,DF]=FRESHE(@PSITh,X); %вычисление правых частей

else [Y]=PSITh(X); %вычисление правых частей

end

else

if s>0.1

[Y,DF]=FRESHE(@PSI,X);

else [Y]=PSI(X);

end

end

WW=DF\Y'; %решение системы нелинейных уравнений

X=X-WW; %вектор деформации

s=0;

for i=1:nfresh %сходимость

s=s+abs(WW(i))./(abs(X(i))+1e-10);

end

end

if (er == 1) && (delta<-0.007) %вычисление параметра вертикального смещения при движении по другому

zzz = PSI2(X);

delta = zzz(2);

end;

if (dlt>0)

zzz = PSI2(X); %вектор левой половинки

z11 = PSI1(X); %вектор правой половинки

z22 = PSI2(X);

uzz = z11(6) + z22(6); %вертикальная суммарная нагрузка

mas = [mas;[-delta,iter,uzz,z11',z22',X']]; %вектор данных

TzR=[TzR,X(3)]; %вектор нагрузки правой половинки

TzL=[TzL,X(6)]; %вектор нагрузки левой половинки

Tz = [Tz,zzz(3)]; %вектор параметра угла поворота

grdlt = [grdlt,-delta]; %вектор параметра вертикального смещения

end

t = 1;

PSI_PRINT1(X); %построение формы правой половинки

PSI_PRINT2(X); %построение формы левой половинки

hold on; grid on;

if (dlt<0) %постоение верхней пластины

subplot(2,2,1)

plot([ln(1),ln(3)],[ln(2),ln(4)],'r');

else

lnln = [lnln;[ln(1),ln(3),ln(2),ln(4)]];

end

hold on; grid on;

[delta,iter]

X1 = X;

X = (X-X0).*(-2*dlt)./(-dlt+(1e-20)) + X0; %экстраполяция начального приближения методом касательных

X0 = X1;

ch = ch + 1;

end

hold on; grid on;

subplot(2,2,2);

for i=1:length(lnln)

plot([lnln(i,1),lnln(i,2)],[lnln(i,3),lnln(i,4)],'r'); %постоение верхней пластины

end;

hold on; grid on;

end

function [] = BALKAM1

%функция, осуществляющая продолжения по параметру (сжатие амортизатора, симметричный случай)

global q1 ii mas1 razb ThetaS distance TzL1 TzR1 grdlt1 delta ab dlt shag ln t ogr

ln = zeros(4);

[R0,Z0,F0,K0,h0,mu,n,beta] = RVO(0);

nfresh=3;

ab=1;

shag=0.05;

kol = 6;

X=zeros(kol,1);

Xp=zeros(kol,1);

X0=zeros(kol,1);

X1=zeros(kol,1);

q1 = 1;

mas1 = [];

TzL1 = [];

TzR1 = [];

ThetaS = [];

grdlt1 = [];

ii = 1;

while (abs(delta)<(ab*sin(F0)-abs(dlt))) && (ogr-abs(dlt)>abs(delta)) && (q1~=0) && (distance<delta)

iter=0; s=1;

Xp = X;

rrr = 0;

delta = delta+dlt;

if (delta>-0.0080001) && (delta<-0.0079999)

dlt = dlt/razb;

end;

if (delta>-0.01200001) && (delta<-0.01199999)

dlt = dlt*razb;

end;

while s>0.0001

iter=iter+1;

if s>0.1

[Y,DF]=FRESHE(@PSI,X);

else [Y]=PSI(X);

end

WW=DF\Y';

X=X-WW;

s=0;

for i=1:nfresh

s=s+abs(WW(i))./(abs(X(i))+1e-10);

end

end

TzR1=[TzR1,X(3)];

TzL1=[TzL1,X(6)];

grdlt1 = [grdlt1,-delta];

ThetaS = [ThetaS,X(1)];

t = 1;

if (X(3)+X(6) > 0.0099) && (X(3)+X(6) < 0.011)

PSI_PRINT1(X);

PSI_PRINT2(X);

end;

[delta,iter]

z11 = PSI1(X);

z22 = PSI2(X);

uzz = z11(6) + z22(6);

mas1 = [mas1;[-delta,iter,uzz,z22',X']];

X1 = X;

X = (X-X0).*(-2*dlt)./(-dlt+(1e-20)) + X0;

X0 = X1;

end

%main

%скрипт, запускающий работу программы

global delta distance1 razb p1 ln mas mas1 lnln ii arr1 arr ThetaS obrshag distance grdlt TzL TzR grdlt1 TzL1 TzR1 ab ThetaL ThetaR dlt deltax ogr

hold off;

[R0,Z0,F0,K0,h0,mu,n,beta] = RVO(0);

dlt = -0.001; %задание шага параметра продолжения

deltax = 0.001; %задание неидеальности

p1 = 0.011; %точка при достижении которой осуществляется смена параметра

obrshag = 1; %разбиение шага при движении обратно

distance = -0.4; %точка предела сжатия

razb = 10; %разбиение шага при движении по другому параметру

delta = -dlt; %начальное значение параметра (обнулится при заходе в цикл)

ogr = ab*sin(F0)+3*dlt; %ограничение для того, чтобы верхняя пластина не перерезала нижнюю

BALKAM; %основная функция несимметричного сжатия

hold on; grid on;

Tt = TzL + TzR; %вектор суммарной вертикальной нагрузки

% yy = [];

% yy=[yy;[grdlt;Tt]];

% hold on; grid on;

subplot(2,2,4); %диаграмма нагрузка-перемещение для правой и левой половинок по отдельности

hold on; grid on;

plot(grdlt, TzL,'r');

plot(grdlt, TzR,'b');

legend('левая','правая',0);

plot(grdlt, TzL,'r*',grdlt, TzR,'*');

ylabel('Tz');

xlabel('delta');

hold on; grid on;

deltax = 0; %обнуление неидеальности

distance = -ogr;

dlt = -0.001;

delta = -dlt;

BALKAM1; %основная функция симметричного сжатия

hold on; grid on; %построение диаграммы нагрузка-перемещения для несиммеричного варианта

subplot(2,2,3);

Tt = TzL + TzR;

plot(grdlt, Tt,'r','LineWidth',2);

%plot(grdlt, Tt,'r*');

ylabel('Tz');

xlabel('Uz');

%yy = [];

%yy=[yy;[grdlt;Tt]];

hold on; grid on;

Tt = TzL1 + TzR1;

plot(grdlt1, Tt,'--b','LineWidth',2); %построение диаграммы нагрузка-перемещение для симметричного варианта

legend('несимметричная','симметричная',0);

%plot(grdlt, Tt,'b*');

hold on; grid on;