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

Работа выполнена при финансовой поддержке гранта РНФ № 23-21-00269,

В исследованиях, посвященных математическому моделированию процесса гашения электрической дуги в спутном газовом потоке, были проведены как аналитические, так и численные решения различных математических моделей гиперболического уравнения теплопроводности [1], выведенного на основе обобщения гипотезы Фурье [2]. Дифференциальное уравнение переноса в частных производных, полученное из классического закона Фурье, в двухмерном пространственном контексте принимает следующий вид:

∂u    ∂2u ∂2u

∂t   a ∂x2   ∂y2 .

Формула (1) является уравнением в частных производных параболического типа, и её точные решения демонстрируют известный характер бесконечно высокой скорости распространения температурного возмущения. Каждый локальный скачок температурной величины провоцирует немедленные изменения во всех точках окружающей среды, независимо от удаленности от отправной точки. Это расходится с теорией относительности и всеми известными принципами передачи тепла. Основанием данной ошибки в логической последовательности является то обстоятельство, что уравнение теплопередачи имеет первый порядок по временной переменной. Второй недостаток, связанный с этим — это ограничение задачи c начальными значениями. В задаче Коши допускается определение температуры T исключительно на начальном этапе времени, то есть при t=0. В ходе экспериментов возможно адаптировать производную температуры по времени ∂T /∂t в момент времени t=0 согласно условиям эксперимента. Это достижимо в гиперболическом уравнении теплопроводности второго порядка во времени. Эксперименты, проведенные с использованием второго звука в твердом гелии и других кристаллических твердых телах при крайне низких температурах и короткой продолжительности, наглядно продемонстрировали, что тепло распространяется в виде затухающей волны. В случае практически бездефектной (идеальной) кристаллической структуры и соблюдения условий для второго звука, после импульсного нагрева можно наблюдать образование горбовидной области с повышенной температурой, которая перемещается с постоянной скоростью через материал [3]. Волна отражается от границ, постепенно теряя свою энергию в процессе движения. Затухающая тепловая волна описывается гиперболическим уравнением в частных производных, которое было впервые выведено Максвеллом, а затем предложено Вер-ноттом и Каттанео [2]. Используя обобщенный закон Фурье, мы можем получить гиперболическое уравнение теплопроводности для двухмерного пространственного случая:

∂2u   ∂u    ∂2u ∂2u

τ ∂t2    ∂t   a ∂x2 ∂y2 .

1 Постановка точной начально-краевой задачи

Вносим изменения в математическую модель, изучая уравнение теплопроводности параболо-гиперболического типа:

∂ ² u

b(x, y, t) dt 2 + a{^x,y,t)

ju =A(dxU + ddyu)+c(x,y,t')u+f(x,y,t), (3)

в области G = [0,X] x [0, Y] x [ — T, T] , T>0. При этом, b(x,y,t) = 0, при t <  0 ; b(x,y,t) > 0,t >  0 ; a(x,y,t)>0, для всех (x,y,t) из области G, т. е. при t <  0 уравнение (3) — параболическое, а при t>0 — гиперболическое.

Начально-краевая задача. Для выполнения конкретного численного расчета в программе Mathcad-15 температурного поля в однородной мембране c X=Y= π и промежутком времени [-Т,Т] используются следующие значения коэффициентов: b(x,y,t)=0 при t <  0 ; b(x,y,t)=1 при t>0; a=1, c(x,y,t)=const<0, A = const >  0 и переменный источник тепла f(x,y,t).

Таким образом уравнение (3) примет вид:

• b (x, y, t)utt + a(x,y,t)ut = A(uxx + Uyy) + c(x,y,t)u + f (x,y,t).    (4)

Начальные условия:

u(x, y, t)\t=T = uo(x, y), uo = 10sinx • siny u(0, y, t) = 0, u(n, y, t) = 0, u(x, 0, t) = 0, u(x, n,t) =0, 0 < x,y < n, —T < t < T.

• 2    Численное решение поставленной задачи по явной конечно-разностной схеме

Мы применим метод конечно-разностной аппроксимации [4] к формулировке начально-краевой задачи (4-6) для получения дискретного решения функции u(x,y,t). Дискретная функция Uik,j , будет определена в узле, заданном координатами: Xi = ih,yj = jh,tk = k.

Рассмотрим явную разностную схему для уравнения (3):

b(x,y,t)

k +1        k k - 1

i,j            i,j i,j

τ2

+ a(x,y,t)

k+1     k i,j          i,j

τ

k    kk    kk i+1,j        i,j i-1,j i,j-1 i,j+1

+ c(^x,y,t)U ^kj + f (x^t).

Для уравнения (4) получим:

k+1      k k-1 k+1 k b i,j          ij + i,j + i,ji,j

τ2

. U i +1 j - -U -j + Uh, j + U kj - 1 + U -j+1

Ah

⁺ c^{U -}j ^{+ f} (x i , y j , t k +1 ).

Из предыдущего уравнения далее выводим расчетную формулу для параболического уравнения при t <  0 :

uj = (1-;  +ст )Uj+h2 (U-+ij+Ui-I,j +Ukj-1+Ukj+1)+Tf (x,yj ,tk+i).

Итоговые формулы для расчета в среде программирования для гиперболического уравнения при t>0 принимают следующий вид:

k +1 i,j

1 + т

(2 + т

-

-Ат 2 h2

+ ст^U-j -

¹ U^-xl+ 1+ T ^j

τ ²

1+т ^f(x i ,y j ^,t - +1 ⁾.

+ (1^ (U-+1j + U-—1j + j + U-+1) +

Проводим дискретизацию первых краевых условий:

Ui-0 = 0,U-n =0,U-j = 0,U-j = 0.

Метод решения данного уравнения с использованием явной разностной схемы в среде MathCAD-15 описан в работах [10, 11, 12].

• 3    Численное решение поставленной задачи по неявной конечно-разностной схеме

Применяя метод конечно-разностной аппроксимации [4] к точной формулировке начально-краевой задачи (4-6), рассмотрим дискретную реализацию по неявной разностной схеме для уравнения (4) с использованием локально-одномерного метода [5, 6]:

иk+1 _ kUk ■ + U-T1 uk+1 - uk ■ i,j            i,j + i,j     . i,j          i,j

+

T2                  T

U k+1

= AU ±1 i

- wj + u k+j + u k+i + u k+i

h2

+ cU^1 + f (xi,yj ,tk+i).

Тогда, формулы для параболического уравнения при t <  0 имеют вид:

Uk+1 _ ^к- i,j          i,j

τ

иk+1

= X^+J

- 4- + U + + j + j

h2

+cUk+1 + f (xi,yj ,tk+1).

Приводим к канонической форме для определения коэффициентов метода прогонки [4].

X (Uk+1 I Uk+1    <1-1- 4Xt k+7k+1 I X (nk+1  I Uk+1 —

h2 (Ui+1,j + Ui-1,j) - (~ + ~h? - c)Ui,j + h2 (UiJ-1 + UiJ+1) =

—    Uij — f (xi,yj ,tk+1)

τ

Обозначим коэффициенты в этом уравнении для дальнейшего использования в программировании:

A — C —

h2

B — -τ

+

4Xt h ² - c,

D — — TUk,j - f(xi,yj ,^+1).

Прогонку проводим сначала в направлении оси Ox, потом в направлении оси Оу. Формулы для гиперболического уравнения при t>0 так же приводим к канонической форме. Домножаем правую и левую часть уравнения на T ² и приводим подобные слагаемые:

Рис. 1. График смещения поля температуры

Uj - 2^^ + Uk-1 + T (Uj - Uj) =

_ __/ rrk +1       k +1 k +1      k+1     k+1      2т fk +1 , 2

= h? (Ui+1,j-4Ui,j + Ui-1,j + Ui,j-1 + Ui,j+1)+cT Ui,j +T f (xi,yj ,tk+1),

_T (kr^1     k+1             4_t      2 k+1 , _T /t++^     k+1 \ _

h2 (Ui-1,j + Ui+1,j) — (1 + T +  h2  — cT )Ui,j + h2 (Ui,j-1 + Ui,j+1)~

— — (2 + T)Uk,j + Uk- 1 — T2f (xi,yj ,tk+1).

Тогда, получаем коэффициенты для использования в системе програм- мирования:

A — C — _T2, h2

B — 1 + T +

4_t 2 h2

-

cτ

D — — (2 + T)Uk,j + Uk- 1 — T2f (xi,yj, tk+1).

Заключение

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