Автор работы: Пользователь скрыл имя, 19 Декабря 2011 в 10:47, лабораторная работа
Такая схема называется явной. Она удобна для организации вычислительного процесса и минимальна по количеству вычислительных операций. Единственным, но очень существенным недостатком явной схемы является то, что она устойчива лишь при весьма жестких ограничениях шага по времени. Мелкое дробление шага , не связанное с требованием точности, приводит к неоправданно большим вычислительным затратам. По этой причине использование явных схем для решения одномерных и тем более многомерных задач малоэффективно.
Второй
полушаг по времени: расчет изменения
температуры на новом временном
слое по горизонтали (по индексу i)
для каждой фиксированной строчке (индекс
k)
.
(30)
Решение систем уравнений (29) и (30) осуществляется методом прогонки.
Прогоночные коэффициенты для расчета изменения промежуточной температуры по вертикали рассчитываются снизу вверх по формулам
Необходимые для начала счета значения определяются из граничного условия - существование постоянного потока q на нижнем слое (k = Nk). В результате имеем
После расчета прогоночных коэффициентов от нижнего слоя к верхнему производится расчет промежуточной температуры от верхнего слоя к нижнему по формуле
, . (31)
При
расчете промежуточной
Полученное
распределение промежуточной
Расчет начинается с вычисления прогоночных коэффициентов «обратным счетом» от правой границы по формулам
Необходимые начальные значения определяются из граничных условий (отсутствие горизонтальной составляющей потока) на правом краю среды. Тогда .
Далее с помощью этих коэффициентов определяется изменение температуры горизонтальных слоев Расчет проводится слева от известного, приведенного в выражении (24), граничного распределения температуры по вертикали
. (32)
Таким
образом, все значения температуры
на (n+1)-м временном слое известны и
можно перейти к вычислению значений
на очередном (n+2)-м слое. При переходе
к (n+2) временному слою матрица значений
температур (n+1) слоя считается начальной
и алгоритм возвращается на начало продольно
поперечной схемы.
Программа расчета поля температур и плотности потока
#include <conio.h>
#include <stdio.h>
#include <stdlib.h>
int i,k,p;
int Ni=100,Nk=60,L=30000;
double hx=0.01,hz=0.0166666666;
double tau=0.001 ;
double lam[102][62];// описание переменных
double q=0.03;
int Q=0;
FILE *f1,*data;
main()
{
double alf[61],bet[61];
double A[61],B[61],C[61],W[61];
double q1[101];
double T[101][61],Tt[101][61],T1[101]
double alf1[101],bet1[101];
double A1[101],B1[101],C1[101],W1[
f1=fopen("a:\\das.txt","w");
data=fopen("a:\\plotn.txt","w"
for(k=0;k<=61;k++)
{
for(i=0;i<=101;i++)
{
if(i>43 && i<=57)
{//задание матрицы лямбда
if(k>=20 && k<28)
lam[i][k]=2;
else lam[i][k]=1;
}
else lam[i][k]=1;
}
}
for (k=0;k<=60;k++)
{
for(i=0;i<=100;i++)//задание начального распределения температуры
{
T[i][k]=(q*hz*k*L)/(2*lam[0][
T1[i][k]=(q*hz*k*
Tt[i][k]=(q*hz*k*
}
}
for(p=1;p<=1000;p++)//начало цикла по времени
{ //start of circle of p (time)
for(i=1;i<=99;i++)//начало прогонки по вертикали
{
for(k=1;k<=60;k++)
{ //вычисление коэф-тов А,В,С,W
A[k]=((2*lam[i][k+1]*lam[i][k]
B[k]=(1/(0.5*tau))+(((2*lam[i]
C[k]=(2*lam[i][k]*lam[i][k-1])
W[k]=(T[i][k]/(0.5*tau))+Q+(((
}
alf[59]=0.9999; //задали граничные альфа и бета
bet[59]=(hz*q*L)/(2*lam[i][60]
for(k=59;k>=1;k--)
{ //вычисление коэф-тов альфа и бета снизу вверх
alf[k-1]=(C[k]/(B[k]-(A[k]*
bet[k-1]=((W[k]+(A[k]*bet[k]))
}
for(k=1;k<=60;k++)
{ //вычисление промежуточной температуры сверху вниз
Tt[i][k]=(alf[k-1]*Tt[i][k-1])
}
} //end of circle i //конец прогонки по вертикали
for(k=1;k<=59;k++)
{//начало прогонки по горизонтали
alf1[99]=0.999999;
bet1[99]=0.0;
for(i=1;i<=100;i++)
{ //вычисление коэф-тов А1,В1,С1,W1
A1[i]=((2*lam[i+1][k]*lam[i][
B1[i]=(1/(0.5*tau))+(((2*lam[
C1[i]=(2*lam[i][k]*lam[i-1][k]
W1[i]=(Tt[i][k]/(0.5*tau))+Q+(
}
for(i=99;i>=1;i--)
{ //вычисление коэф-тов альфа-1 и бета-1 слева направо
alf1[i-1]=(C1[i]/(B1[i]-A1[i]*
bet1[i-1]=((W1[i]+A1[i]*bet1[
}
for(i=1;i<=100;i++)
{ //вычисление искомой температуры справа налево
T1[i][k]=(alf1[i-1]*T1[i-1][k]
}
} //end of circle k
//конец прогонки по
горизонтали
for(i=1;i<=99;i++)
{
T1[i][60]=Tt[i][
}
for(i=1;i<=100;i++)
{
for(k=1;k<=60;k++)
{//переобозначение для возврата на начало
T[i][k]=T1[i][k];
Tt[i][k]=T1[i][k];
}
}
}//end of circle of p//конец цикла по времени
for(k=0;k<=60;k++)
{
for(i=0;i<=100;i++)
{
//вывод конечной матрицы температур
fprintf(f1,"%f ",T1[i][k]);
}
fprintf(f1,"\n");
}
for(i=0;i<=100;i++)
{//вычисление плотности потока между нулевым и первым уровнями
q1[i]=((T1[i][1]-T1[i][0])*2)/
fprintf(data,"%.10f\n",q1[i]);
}
fclose(f1);fclose(data); return 0;}