S(i,j)=S0;

end

end

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

for j=1:m

u(1,j)=0;v(1,j)=0;C(1,j)=C1e;

u0(1,j)=0;v0(1,j)=0;C0(1,j)=C1e;

end

for j=1:n

u(i,1)=0;v(i,1)=0;C(i,1)=C1e;

u0(i,1)=0;v0(i,1)=0;C0(i,1)=C1e;

end

% Задание данных для нестационарной задачи

dt=0.01; tk=10; t=0;

% Начало временного цикла

while t
%Вычисление коэффициентов уравнений b и c

for i=2:n-1

for j=2:m-1

Fe=r0*u(i,j);

Fw=r0*u(i-1,j);

De=D*dy(j)/hx(i);

Dw=D*dy(j)/hx(i-1);

Pe=Fe/De;

Pw=Fw/Dw;

Ape=max(0,(1-0.1*Pe)^5);

Apw=max(0,(1-0.1*Pw)^5);

b(i,j)=De*Ape+max(-Fe,0);

c(i,j)=Dw*Apw+max(Fw,0);

end

end

%коэффициентов уравнений d и e

for i=2:n-1

for j=2:m-1

Fn=r0*v(i,j);

Fs=r0*v(i,j-1);

Dn=D*dx(i)/hy(j);

Ds=D*dx(i)/hy(j-1);

Pn=Fn/Dn;

Ps=Fs/Ds;

Apn=max(0,(1-0.1*Pn)^5);

Aps=max(0,(1-0.1*Ps)^5);

d(i,j)=Dn*Apn+max(-Fn,0);

e(i,j)=Ds*Aps+max(Fs,0);

end

end

%Вычисление коэф. для прохода по горизонтальным линиям

for i=2:n-1

for j=2:m-1

ap0(i,j)=2*r0*dx(i)*dy(j)/dt;

a(i,j)=b(i,j)+c(i,j)+ap0(i,j);

d1(i,j)=d(i,j)*C0(i,j+1)+e(i,j)*C0(i,j-1)-(d(i,j)+e(i,j)-ap0(i,j))*C0(i,j)+S(i,j)*dx(i)*dy(j);

end

end

%Проход по горизонтальным линиям

%Открываем цикл по j

for j=2:m-1

%При x=L:dc/dx=0 учтем это Гу

%Уравнение на правой границе а(n,j)*C(n,j)=c(n,j)*C(n-j)

%Тогда

a(n,j)=1;

c(n,j)=1;

b(n,j)=0;

d1(n,j)=0;

P(1)=0;

Q(1)=0;

for i=2:n

P(i)=b(i,j)/(a(i,j)-c(i,j)*P(i-1));

Q(i)=(d1(i,j)+c(i,j)*Q(i-1))/(a(i,j)-c(i,j)*P(i-1));

end

C(n,j)=Q(n);

i=n-1;

while i>1

C(i,j)=P(i)*C(i+1,j)+Q(i);

i=i-1;

end

for i=2:m

C0(i,j)=C(i,j);

end

end

for i=2:n-1

for j=2:m-1

a(i,j)=d(i,j)+e(i,j)+ap0(i,j);

d1(i,j)=b(i,j)*C0(i+1,j)+c(i,j)*C0(i-1,j)-(b(i,j)+c(i,j)-ap0(i,j))*C0(i,j)+S(i,j)*dx(i)*dy(j);

end

end

for i=2:n-1

%при y=h: dc/dy=0 учтем это ГУ

%уравнение на верхней границе a(i,m)*C(i,m)=e(i,m)*c(i,m-1);

a(i,m)=1;

e(i,m)=1;

d(i,m)=0;

d1(i,m)=0;

P(1)=0;

Q(1)=0;

%Решаем TDMA систему уравнений:

for j=2:m

P(j)=d(i,j)/(a(i,j)-e(i,j)*P(j-1));

Q(j)=(d1(i,j)+e(i,j)*Q(j-1))/(a(i,j)-e(i,j)*P(j-1));

end

C(i,m)=Q(m);

j=m-1;

while j>1

C(i,j)=P(j)*C(i,j+1)+Q(j);

j=j-1;

end

for j=2:m

C0(i,j)=C(i,j);

end

end

t=t+dt;

end

%Рисование графика на данный момент времени

contour(x,y,C','k')

xlabel('x,м'); ylabel('y,м')

---

Смотрите также файлы

• огневые тренировки и стрельбы.docx

• Методические указания по практической работе для студентов направлений 20. 03. 01 Техносферная безопасность.docx

• Практическая работа к теме 1 Ценности и их роль для человека и общества Темы занятий День Земли.docx

• Создание комплексноцелевой программы управления качеством образования.docx

• Сколько чашек теперь там.docx