Уравнение теплопроводности с подвижной границей. Ур-е теплопроводности с подвижной границей. #Одномерное уравнение теплопроводности с подвижной границей (промерзание влажного грунта)
 Скачать 16.16 Kb.
|
|
#`Одномерное уравнение теплопроводности с подвижной границей (промерзание влажного грунта)` restart; with(plots); N:=10; #`количество пространственных узлов` tk := 300; #`конечное время ` L := 2; #`нижняя граница (глубина грунта) x (в cм)` To := 275.15; # начальная температура грунта (в Кельвинах) (2^(0)C) Tc := 271.15; # установившаяся температура на границе в начальный момент времени, которая ниже температуры замерзания (-2^(0)C) Tz := 271.65; # температура замерзания грунта (-1.5^(0)C) Q := 400000; #`теплота фазового перехода (Дж|кг) #для воды 333000` c1 := 2.5; #`удельная теплоемкость промерзшей зоны грунта ` c2 := .8; #`удельная теплоемкость талой зоны грунта ` p1 := 450; #`плотность промерзшей зоны грунта (кг|м^(3))` p2 := 300; #`плотность талой зоны грунта (кг|м^(3))` lyamb1 := 1.0; #`коэфициент теплопроводности промерзшей зоны грунта` lyamb2 := 0.85; #`коэфициент теплопроводности талой зоны грунта` a[1] := 9.5; #`Коэффициент температуропроводности промерзшей зоны грунта` a[2] := 7; #`Коэффициент температуропроводности талой зоны грунта` w := 0.3; #`влажность грунта ` eps := 0.001; #`Определяем поле температуры в начальный момент времени :` tn := 0; for i to N do T[i] := To: od: #` Шаг по пространственной координате:` h := evalf(L/(N-1)); k := 1: #`Положение границы фазового перехода` #`Проводим интегрирование нестационарного уравнения теплопроводности:` #`Используем цикл с предусловием while: ` while(tn #`1) запоминаем поле температуры на предыдущем временном слое:` for i to N do Tn[i] := T[i] : od: #`Граница фазового перехода на каждом соответствующем временном шаге смещается на пространственный шаг вправо` k := k+1: #`Цикл с постусловием, позволяющий итерационно вычислять поле температуры,вследствие наличия нелинейностив граничном условии:` for i to N do Ts[i] := T[i]: od: #` Определяем соответствующий шаг по временной координате:` tau :=  :#`Определяем начальные прогоночные коэффициенты на основе левого граничного условия:` alfa[1] := 0: betta[1] := Tc: #`Цикл с параметром для определения прогоночных коэффициентов по формуле (8):` for i from 2 to (k-1) do #`a[i], b[i], c[i], f[i] - коэффициенты канонического представления системы уравнений с трехдиагональной матрицей:` a[i] := a[1]/sqrt(h): b[i] := 2*a[1]/sqrt(h)+1/tau: c[i] := a[1]/sqrt(h): f[i] := -Tn[i]/tau: #` alfa[i], betta[i] - прогоночные коэффициенты` alfa[i]:=(a[i])/((b[i]-c[i]*alfa[i-1])): betta[i]:=(c[i]*betta[i-1]-f[i])/((b[i]-c[i]*alfa[i-1])): od: #`Определяем значение температуры на границе фазового перехода:` T[k] := Tz: #`Используя отношение (7) определяем неизвестное поле температуры в промерзшей зоне грунта:` for i from k-1 by -1 to 1 do T[i] := T[i+1]*alfa[i]+betta[i]: od: #`Определяем начальные прогоночные коэффициенты для талой зоны грунта на основе условия на границе раздела двух сред:` alfa[k] := 0: betta[k] := Tz: #`Цикл с параметром для определения прогоночных коэффициентов по формуле (8):` for i from k+1 to N-1 do a[i] := a[2]/sqrt(h); b[i] := 2*a[2]/sqrt(h)+1/tau; c[i] := a[2]/sqrt(h); f[i] := -Tn[i]/tau; alfa[i] := a[i]/(-alfa[i-1]*c[i]+b[i]); betta[i] := (betta[i-1]*c[i]-f[i])/(-alfa[i-1]*c[i]+b[i]) ; end do: #`Определяемзначение температуры на правой границе на основе граничного условия, используя соотношение (21):` T[N] := (2*a[2]*tau*betta[N-1]+sqrt(h)*Tn[N])/(2*a[2]*tau*(1-alfa[N-1])+sqrt(h)): #`Используя соотношение (7) определяем неизвестное поле температуры в талой зоне грунта:` for i from N-1 by -1 to k do T[i] := T[i+1]*alfa[i]+betta[i]: od: #`Определяем максимум модуля разности между соответствующими значениями температуры на текущей и предыдущей итерации:` Max:=abs(T[1]-Ts[1]): for i from 2 to N while Max>eps do if (Max Max:=abs(T[i]-Ts[i]): end if: od: tn:=tn+tau: print(tn); od: Результат: -3.276170445 10^6 -1.247835529 10^14 1.247835465 10^14 |