Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- uses crt;
- const mf=102;
- type
- vector1=array[1..mf] of real;
- vector2=array[1..mf,1..mf] of real;
- var {раздел описания переменных, которые мы будем использовать в
- программе}
- i, j, Nx, Ny : integer;
- T : vector2;
- alfa, beta : vector1;
- ai, bi, ci, fi : real;
- a, lamda, ro, c : real;
- hx, hy, tau, t_end, time : real;
- T0, L, H, Th, Tc : real;
- f, g : text;
- begin
- clrscr;
- {с клавиатуры вводим все необходимые входные параметры}
- Writeln('Введите количество пространственных узлов в пластине по оси х, Nx');
- Readln(Nx);
- Writeln('Введите количество пространственных узлов в пластине по оси y, Ny');
- Readln(Ny);
- Writeln('Введите окончание по времени, t_end');
- Readln(t_end);
- Writeln('Введите длину пластины, L');
- Readln(L);
- Writeln('Введите толщину пластины, H');
- Readln(H);
- Writeln('Введите коэффициент теплопроводности материала пластины, lamda');
- Readln(lamda);
- Writeln('Введите плотность материала пластины, ro');
- Readln(ro);
- Writeln('Введите теплоемкость материала пластины, c');
- Readln(c);
- Writeln('Введите температуру на границе х = 0 области решения, Th');
- Readln(Th);
- Writeln('Введите температуру на границе х = L области решения, Tc');
- Readln(Tc);
- Writeln('Введите начальную температуру, T0');
- Readln(T0);
- {определяем расчетные шаги сетки по пространственным координатам}
- hx:=L/(Nx-1);
- hy:=H/(Ny-1);
- {определяем коэффициент температуропроводности}
- a:=lamda/(ro*c);
- {определяем расчетный шаг сетки по времени}
- tau:=t_end/100.0;
- {определяем поле температуры в начальный момент времени}
- for i:= 1 to Nx do
- for j:= 1 to Ny do
- T[i,j]:=T0;
- {проводим интегрирование нестационарного уравнения
- теплопроводности}
- time:=0;
- while time<t_end do {используем цикл с предусловием}
- begin
- {увеличиваем переменную времени на шаг τ}
- time:=time+tau;
- {решаем СЛАУ в направлении оси Ох для определения поля температуры на промежуточном временном слое}
- for j:=1 to Ny do
- begin
- {определяем начальные прогоночные коэффициенты на основе левого граничного условия}
- alfa[1]:=0.0;
- beta[1]:=Th;
- {цикл с параметром для определения прогоночных коэффициентов по формуле (8)}
- for i:= 2 to Nx-1 do
- begin
- {ai, bi, ci, fi – коэффициенты канонического представления СЛАУ с трехдиагональной матрицей}
- ai:=lamda/sqr(hx);
- bi:=2.0*lamda/sqr(hx)+ro*c/tau;
- ci:=lamda/sqr(hx);
- fi:=-ro*c*T[i,j]/tau;
- {alfa[i], beta[i] – прогоночные коэффициенты}
- alfa[i]:=ai/(bi-ci*alfa[i-1]);
- beta[i]:=(ci*beta[i-1]-fi)/(bi-ci*alfa[i-1]);
- end;
- {определяем значение температуры на правой границе на основе правого граничного условия}
- T[Nx,j]:=Tc;
- {используя соотношение (7) определяем неизвестное поле температуры на промежуточном временном слое}
- for i:= Nx-1 downto 1 do
- T[i,j]:=alfa[i]*T[i+1,j]+beta[i];
- end;
- решаем СЛАУ в направлении оси Оу для определения поля температуры на целом временном слое}
- for i:=2 to Nx-1 do
- begin
- {определяем начальные прогоночные коэффициенты на основе нижнего граничного условия, используя соотношения (20) при условии, что q1 = 0}
- alfa[1]:=2.0*a*tau/(2.0*a*tau+sqr(hy));
- beta[1]:=sqr(hy)*T[i,1]/(2.0*a*tau+sqr(hy));
- {цикл с параметром для определения прогоночных коэффициентов по формуле (8)}
- for j:= 2 to Ny-1 do
- begin
- {ai, bi, ci, fi – коэффициенты канонического представления СЛАУ с трехдиагональной матрицей}
- ai:=lamda/sqr(hy);
- bi:=2.0*lamda/sqr(hy)+ro*c/tau;
- ci:=lamda/sqr(hy);
- fi:=-ro*c*T[i,j]/tau;
- {alfa[j], beta[j] – прогоночные коэффициенты}
- alfa[j]:=ai/(bi-ci*alfa[j-1]);
- beta[j]:=(ci*beta[j-1]-fi)/(bi-ci*alfa[j-1]);
- end;
- {определяем значение температуры на верхней границе, используя соотношение (21) при условии, что q2 = 0}
- T[i,Ny]:=(2.0*a*tau*beta[Ny-1]+sqr(hy)*T[i,Ny])/(2.0*a*tau
- *(1.0-alfa[Ny-1])+sqr(hy));
- {используя соотношение (7) определяем неизвестное поле температуры на промежуточном временном слое}
- for j:= Ny-1 downto 1 do
- T[i,j]:=alfa[j]*T[i,j+1]+beta[j];
- end;
- end; {цикл с предусловием окончен}
- {выводим результат в файл}
- Assign(f,'res.txt');
- Rewrite(f);
- Writeln(f,'Длина пластины L = ',L:6:4);
- Writeln(f,'Толщина пластины H = ',H:6:4);
- Writeln(f,'Число узлов по пространственной координате x в пластине Nx = ',Nx);
- Writeln(f,'Число узлов по пространственной координате y в пластине Ny = ',Ny);
- Writeln(f,'Коэффициент теплопроводности материала пластины lamda =',lamda:6:4);
- Writeln(f,'Плотность материала пластины ro = ',ro:6:4);
- Writeln(f,'Теплоемкость материала пластины с = ',c:6:4);
- Writeln(f,'Начальная температура T0 = ',T0:6:4);
- Writeln(f,'Температура на границе x = 0 области решения Th = ',Th:6:4);
- Writeln(f,'Температура на границе x = L области решения Tc = ',Tc:6:4);
- Writeln(f,'Результат получен с шагом по координате x hx = ',hx:6:4);
- Writeln(f,'Результат получен с шагом по координате y hy = ',hy:6:4);
- Writeln(f,'Результат получен с шагом по времени tau = ',tau:6:4);
- Writeln(f,'Температурное поле в момент времени t = ',t_end:6:4);
- close(f);
- Assign(g,'tempr.txt');
- Rewrite(g);
- for i:=1 to Nx do
- for j:=1 to Ny do
- writeln(g,' ',hx*(i-1):10:8,' ',hy*(j-1):10:8,' ',T[i,j]:8:5);
- close(g);
- end.
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement