Уравнения в частных производных. Численное решение уравнений параболического типа. Понятие о методе конечных разностей (МКР). Основные определения и конечно-разностные схемы

Уравнения в частных производных Многие задачи механики и физики приводят к исследованию дифференциальных уравнений с частными производными второго порядка. Их изучение дает возможность построить теорию широкого круга физических явлений и решить ряд физических и технических задач, которые встречаются в различных областях науки и приложениях. Уравнения в частных производных получаются огромным числом способов и возникают в самых разнообразных физических задачах, тем не менее, такие уравнения и системы уравнений очень часто принимают одинаковую или близкую форму. Это весьма существенно, так как множество принципов в физике, особенно в классической, сформулировано в виде законов сохранения. Масса не создается и не уничтожается, импульс, энергия, поток магнитного поля, полный электрический сохраняются. Например, уравнение диффузии часто появляется при описании переноса частиц импульса или энергии. В частности, оно описывает распределение температуры в твердом теле, где в классическом случае энергия переносится посредством теплопроводности. По закону сохранения энергии скорость изменения энергии в объеме V должна равняться потоку энергии q через поверхность S ограничивающую объем V плюс мощность источников энергии действующих в этом объеме и имеющих плотность J.  (1)  x , t dV   q  n dS  JdV       t  V S V Применяя теорему Гаусса–Остроградского   q  n dS     qdV S V    (2)    q  Q dV  0   t  V Так как объем V может быть произвольным то (2) будет выполняться, только если подынтегральное выражение для любого момента времени в каждой точке пространства равно нулю  (3)    q  J . t Плотность тепловой энергии среды T    cdT (4) q  kT , (5) где  - плотность (кг/м3), с удельная теплоемкость (Дж/кг·град). Из эксперимента известно, что поток тепла прямо пропорционален градиенту температуры (закон Фурье) где k – коэффициент теплопроводности (Вт/м·град). После подстановки выражений (4) и (5) в (3) получим уравнение диффузии температуры или уравнение теплопроводности   T T (6)   c    k T  J t T t t К аналогичным по форме уравнениям приводят рассмотрение и ряда других физических явлений таких, например, как процессы диффузии материалов или диффузии электромагнитного поля. Дифференциальные уравнения типа (6), являются нестационарными дифференциальными в частных производных второго порядка параболического типа и представляют краевую задачу математической физики которая может быть решена при задании начальных и граничных условий. Порядок уравнения определяется максимальным порядком производных входящих в уравнение. В декартовой системе координат x, y, z оператор  (набла) имеет вид     i j k x y z где i, j, k – единичные векторы направленные вдоль координатных осей. Градиент скалярной функции f и дивергенция векторной функции v  [vx , vy , vz ]T запишутся так f f f grad f  f  i j k x y z vx vy vz divv    v    x y z и соответственно раскрывая дивергенцию градиента в (6) получим c T T   T    k T  c  k t t x  x    T   k  y  y    T   z  k z     J  Численное решение уравнений параболического типа. Понятие о методе конечных разностей (МКР). Основные определения и конечно-разностные схемы Схема дискретизации. В течение последних десятилетий МКР является одним из широко используемых численных методов решения уравнений с частными производными (уравнений математической физики). Его популярность объясняется прежде всего относительной простотой математической идеи дискретизации. Под дискретизацией понимается аппроксимационная процедура, при которой непрерывная область заменяется сеткой из отдельных точек, и неизвестные физические переменные отыскиваются не в произвольных точках области, а именно в этих конкретных точках. Частные производные при таком подходе могут аппроксимироваться с помощью конечных разностей несколькими способами. Все эти аппроксимации получаются путем разложения соответствующих функций в ряд Тейлора. Порядок точности аппроксимации зависит от ошибок, возникающих при усечении ряда. Рассмотрим задачу расчета температуры в бесконечно протяженной в двух измерениях (y и z) пластину (или стену) конечной толщины l в направлении x. Температура не зависит от y и z, а меняется только в направлении x. Первоначальное распределение температуры в стенке известно, а дальнейшее изменение этого распределения описывается нестационарным уравнением теплопроводности T  2T T  2T (1) c  , or a t t x 2 x 2 здесь предполагается отсутствие внутренних источников тепла и конвективного теплопереноса и полагается, что температуропроводность a   / c . Введем безразмерные координату ξ и время τ которые определим соотношениями   x / xs ,   t / ts и положим xs  l , ts  l 2 / a  cl 2 /  . Привести температуру Т к безразмерному виду несложно, обычно эту процедуру проводят по-разному в зависимости от задаваемых граничных условий. Приведенную к безразмерному виду величину Т обозначим через θ=T/Ts (Ts выбранный масштаб для температуры), тогда окончательно будем иметь уравнение:   2  2 0 (2)   Введем пространственно-временную сетку с равномерным шагами ξ=h по координате и  = k по времени, ее можно изобразить графически: τM τj k τ1=0 ξ1=0 ξi h ξN-1 ξN=1 Точки пересечения вертикальных и горизонтальных линий сетки называются узлами пространственно-временной сетки. Задача, которая решается численными методами – это найти приближенные значения температуры в узлах пространственновременной сетки, т.е. в точках тела ξi = (i –1)∙h в моменты времени j = (j - 1)∙k. Температуры в узлах пространственно-временной сетки называются сеточной функцией. В отличие от точных значений температуры θ для сеточной функции вводят обозначения ui,j; θ(ξi,j) – точные значения температуры в узлах сетки; ui,j – приближенные значения температуры в тех же узлах. Множество всех узлов сетки, имеющих одну и ту же временную координату, называется временным слоем. Узлы сетки, принадлежащие нулевому временному слою (1 =0), а также границам области ξ1=0 и ξN=1, называются граничными, остальные –внутренними. Значения сеточной функции в граничных узлах определены начальными и граничными условиями. Для определения ui,j во внутренних узлах сетки нужно получить уравнения, исходя из исходной 2задачи (2). Для этого точные дифференциальные операторы  и  2 заменяют их   разностными аналогами (приближенными выражениями), используя значения сеточной функции в соседних узлах. При этом возникает вопрос точности аппроксимации, которая может быть различной в зависимости от выбранного способа аппроксимации. Приведем некоторые из них. Некоторые способы аппроксимации дифференциальных операторов Производная        2  2   Разностная аппроксимация ui1, j  uij h uij  ui1, j h ui 1, j  ui 1, j (вперёд) 0(h) (назад) 0(h) (центральная) 2h ui 1, j  2 u ij  ui 1, j h2 ui , j 1  ui , j k Порядок точности (вперёд) 0(h2) 0(h2) 0(k) Результаты представленные в таблице следуют из разложения функции θ (ξ,) в ряд Тейлора как функции ξ (при фиксированном ) или как функции  (при фиксированном ξ).  1  2 (3) 3  (  h, )   ( , )  h h 2  O( h )   2!  2   1  2 3  (  h, )   ( , )  h  h 2  O( h )   2!  2  (4)  (  h, )   ( , )  Из (3)   O ( h) разность вперед имеем h  Из (4)  ( , )   (  h, )  разность назад   O ( h) имеем h  Вычитая из разложения (3) разложение (4) найдем  (  h, )   (  h, )  центральная разность  O( h 2 ) 2h  Суммируя (3) и (4) получим аппроксимацию для второй производной   (  h, )  2 ( , )   (  h, ) h 2  2  2  O( h 2 )  (5) Поскольку (5) – приближенное выражение, использующее значения функции в узлах сетки, вместо  в (2) надо подставить значения сеточных функций  2 ui 1, j  2ui, j  ui 1, j Итак, (6)   2 h2  ui , j 1  ui, j Аналогично, (7)   k Определим теперь, какое количество точек участвует в разностной аппроксимации уравнения (2). Во-первых, для аппроксимации  мы  должны взять два соседних временных слоя : j и j+1. 2 Для аппроксимации   на j – том слое требуется использовать 3 точки:  2 i-1, i, i+1. Итого, мы использовали 4 соседние точки. i Множество точек, используемых для j+1 аппроксимации дифференциального уравнения его разностным аналогом, j называется шаблоном. В данном случае – 4-х i-1 i i+1 точечный шаблон. Используя (6) и (7) найдем конечно разностную аппроксимацию уравнения (2): ui, j 1  ui, j ui 1, j  2ui, j  ui 1, j  2 k h откуда ui, j 1  s  ui 1, j  ui, j (1  2s)  s ui 1, j , (8) (9) где s = k/h2 и индексы принимают значения i = 2,3,…N-1, j = 1,2,…M. Полученная схема решения (9) является явной схемой. Действительно из (9) следует, что если на предыдущем (j –том) слое нам известны значения ui во всех точках пространственной сетки, то на последующем слое j+1 мы их с легкостью определим, прогоняя i от i = 2 до i = N-1 – см. шаблон. Далее, при j=1 (нулевой временной слой) нам заданы все значения температуры из НУ ( ui,1.= θ(ξi,0) ). Стало быть, во всех внутренних точках последующего временного слоя мы определим значения температуры явно по итерационной формуле (9). Для определения температуры в граничных узлах (i = 1 и i = N) необходимо использовать граничные условия . Граничные условия для уравнения теплопроводности ГУ 1-го рода – задается распределение температуры на поверхностях (или границах) области в течение всего времени процесса : Тгр = φ(t); ГУ 2-го рода – задается плотность теплового потока q(t) на поверхностях (или границах) области. ГУ 3-го рода – задан закон теплообмена с окружающей средой на поверхностях (или границах) области: qгр = (Tгр-Тср); температура внешней среды Тср= const, -коэффициент теплоотдачи, заданная величина. Однако, чаще всего задают ГУ смешанного типа, т.е. типы ГУ на различных частях границ области в которой ищется решение не совпадают.   2 Пример решения уравнения   0,   2 на пространственно временной сетке ξi = (i –1)∙h, i = 1,..N, τj = (j –1)∙k, с начальным условием  ( ,0)  4 (1   ) и граничными условиями u1, j  uN , j  0 Используя метод конечных разностей, аппроксимируем производные:  ui 1, j  ui , j  ,  h  ui , j 1  ui , j  ,  k  2 ui 1, j  2ui , j  ui 1, j  . 2 2  h Тогда аппроксимацию уравнения можно записать так ui, j 1  ui, j k  ui 1, j  2ui, j  ui 1, j h2 или ui, j 1  sui 1, j  (1  2s)ui, j  sui 1, j , s  k / h2 function lab6 N = 11; % число пространственных точек ksi = linspace(0,1,N)'; h = ksi(2)-ksi(1); % шаг пространственной сетки u0 = 4*ksi.*(1-ksi);% начальное распределение температуры s = 0.5; % максимально возможное значение k = s*h^2; % шаг повремени M = 31; % число временных слоев (шагов по времени) tauend = (M-1)*k; tau = linspace(0,tauend,M); U = zeros(N,M); % массив для записи решения U(:,1) = u0; for j=2:M u1(2:N-1) = s*u0(1:N-2)+(1-2*s)u0(2:N-1)+s*u0(3:N); U(:,j) = u1; u0 = u1; end plot(ksi,U,'b-','LineWidth',2) set(gca,'FontSize',12) xlabel('\it\tau','FontSize',14) ylabel('\itu','FontSize',14) Распределение температуры в области в различные моменты времени рассчитанные по явной схеме при s = 0,5. Расчет распределения температуры в пластине при s = 0,7. % % % % % % % % % % % % % % % % % % % % % % % u05 = U(:,6); dT = 30; Tp = 20; T05 = l = 0.5; % wall thickness lambda1 = c1 = ro1 = timescale1 = t1 = tau*timescale1; lambda2 = c2 = ro2 = timescale2 = t2 = tau*timescale2; plot(t1,T05,'r-',t2,T05,'b-','LineWidth',2) set(gca,'FontSize',12) legende('--','--') xlabel('time, s','FontSize',14) ylabel('temperature, grad','FontSize',14) Метод конечных разностей. Неявная схема. Шаблон для неявной схемы j+1 i-1 i i+1 В этом случае аппроксимация уравнения (1) для пространственных узлов i = 2,..,N-1 сетки ui , j 1  ui , j k j i-1 i i+1  ui 1, j 1  2ui , j 1  ui 1, j 1 Уравнения (4) образуют алгебраических уравнений h 2 систему u1, j 1  (2  s )u 2, j 1  u3, j 1  su 2, j  u 2, j 1  (2  s )u3, j 1  u 4, j 1  su 4, j  h2   , s  k u  (2  s )ui , j 1  ui 1, j 1  sui , j  i 1, j 1    uN  2, j 1  (2  s )uN 1, j 1  uN , j 1  suN , j . (4) линейных (5) Система уравнений (5) содержит N-2 уравнения. Очевидно, что для определения всех N значений сеточной функции u на временном слое j+1 ее необходимо дополнить еще двумя уравнениями используя граничные условия. Запишем систему уравнений (5) в матричном виде, добавив уравнения для узлов i=1 и i=N.  b1   c2 0  0 0  a1 b2 a2 c3 b3 a3 cN 1 bN 1 aN 1 cN bN   u1, j 1   d1      u  d   2, j 1   2     u3, j 1    d3       u d   N 1, j 1   N 1   u      N , j 1   d N  (6) Для системы (5) и значения коэффициентов в (6) есть c2  c3    cN 1  1, a 2  a3    aN 1  1, b2  b3    bN 1  (2  s ), di  sui , j , i  2,..., N  1. Коэффициенты b1, a1, d1 необходимо определить исходя из граничных условий на левой границе (ξ = 0). Коэффициенты сN, bN, dN из граничных условий на правой границе (ξ = 1). Заметим, что в матрице системы (6) отличны от нуля только коэффициенты расположенные вдоль трех диагоналей. Такая система уравнений эффективно может быть решена используя метод прогонки. Решение системы линейных уравнений с трехдиагональной матрицей методом прогонки Для узла i уравнение системы (6) имеет вид ciui 1  biui  aiui 1  di положим ui 1  i 1ui  i 1 и подставим в (7), получим отсюда (cii 1  bi )ui  aiui 1  di  ci i 1 ai d c  ui  ui 1  i i i 1 . (cii 1  bi ) (cii 1  bi ) (7) (8) (9) Это выражение будет иметь вид аналогичный (8) если положить i  d c  , i  i i i 1 . (cii 1  bi ) (cii 1  bi ) ai ui  iui 1  i (10) (11) Полагая с1=0, aN=0 и начиная с α1=-a1/b1 , β1=d1/b1 все значения αi, βi для i =2,…N можно найти используя (10). Нахождение этих коэффициентов называют прямым ходом метода прогонки . Учитывая αN=0, из (11) имеем uN=βN . Все остальные значения ui для i = N-1,…,1 находятся с использованием (11). Эту процедуру называют обратным ходом метода прогонки. Пример реализации метода прогонки в MatLab function u = marching(a, b, c, d) % marching - прогонка N = length(b); % число точек сетки по x u = zeros(N,1); % вектор для записи решения alpha = zeros(N,1); alpha(1) = -a(1)/b(1); beta = zeros(N,1); beta(1) = d(1)/b(1); % прямой ход for i = 2:N g = 1/(c(i)*alpha(i-1)+b(i)); alpha(i) = -a(i)*g; beta(i) = (d(i)-c(i)*beta(i-1))*g; end u(N) = beta(N); % обратный ход for i = N-1:-1:1 u(i) = alpha(i)*u(i+1)+beta(i); end Коэффициенты b1, a1, d1 и сN, bN, dN для ГУ 1-го, 2-го и 3-го рода Первое и последнее уравнения системы (6) имеют вид b1u1  a1u2  d1 , cN uN 1  bN uN  dN . Необходимо определить входящие в эти уравнения коэффициенты исходя из заданных граничных условий и выбранных масштабов. ГУ 1-го рода. На границах x=0 или x =l заданы температуры T x 0  T0 (t ), T x l  T1(t ) в безразмерном виде при заданных масштабах для времени ts , температуры Ts   0  T0 ( ts ) / Ts ,   1  T1( ts ) / Ts таким образом для значений сеточных функций на границах имеем b1u1  T0 ( ts ) / Ts , bN uN  T1(ts ) / Ts и следовательно b1 = 1, a1 = 0, d1 = T0(τ ts)/Ts = φ0(τ) и сN = 0, bN= 1, dN = T1(τ ts)/Ts = φ1(τ) ГУ 2-го рода. На границах x=0 или x =l заданы тепловые потоки.  T x x 0  q0 (t ),  T x x l  q1(t ), В безразмерном виде q ( t )l q ( t )l     0 s  S 0 ( ),   1 s  S1( )   0 Ts x  1 Ts Для сеточных функций получим u1  u2  S0 ( )h , uN 1  uN  S1( )h помним первое и последние уравнения системы (6) b1u1  a1u2  d1 , cN uN 1  bN uN  dN . и следовательно a1  1, b1  1, d1  S0 ( )h , cN  1, bN  1, dN  S1( )h   2  2  0 с начальным условием  ( , ) Решаем уравнение 0  0    и граничными условиями  1,   1  0.   0 В соответствии с заданными граничными условиями находим b1  1, a1  1, d1  h , cN  0, bN  1, dN  0. Распределение температуры в пластине в различные моменты времени при ее нагреве постоянным тепловым потоком ГУ 3-го рода. На границах x=0 или x =l задан закон теплообмена поверхности со средой с постоянной температурой T T   (Tcp1 T x  0),     (T x  l Tcp2 ), x x  0 x x  l Положим T = Tsθ+Tc и соответственно θ = (T- Tc)/ Ts тогда в безразмерном виде эти условия запишутся так  (Tcp2 Tc )    l  (Tcp1 Tc )  l        0  ,      1  .    0   Ts    T s     1  l   (T T ) / T ,   (T T ) / T Обозначим Bi  , cp1 cp1 c s cp2 cp2 c s , тогда  для сеточной функции граничные условия дадут   u1  u2  Bi h cp1  u1 ,   uN 1  uN  Bi h uN  cp 2 , (1  Bi h )u1  u 2  Bi h cp1, uN 1  (1  Bi h )uN  Bi hcp 2 что позволяет определить коэффициенты b1u1  a1u2  d1 , cN uN 1  bN uN  dN . b1  (1  Bi h ), a1  1, d1  Bi h cp1 cN  1, bN  (1  Bi h ), dN  Bi h cp2 2     Решаем уравнение  2  0 с начальным условием  ( , )  0  1     и граничными условиями  0,   Bi  1 ,   0   1 Используя метод конечных разностей, аппроксимируем производные:  ui 1, j 1  ui , j 1  ,  h  ui , j 1  ui , j  ,  k  2 ui 1, j 1  2ui , j 1  ui 1, j 1  . 2 2  h неявная схема Тогда аппроксимации уравнения для узлов i = 2,..N-1 будет иметь вид ui 1, j 1  (2  s)ui, j 1  ui 1, j 1  sui, j , s  h2 / k. Из граничного условия на левой границе (ξ=0) u1, j 1  u2, j 1  0. Граничное условие на правой границе (ξ=1) дает uN 1, j 1  (1  Bih )uN , j 1  0 В соответствии с заданными граничными условиями находим b1  1, a1  1, d1  0, cN  1, bN  (1  Bi h ), dN  0. Распределение безразмерной температуры в пластине а различные моменты времени. Красные кривые Bi = 10, синие кривые Bi = 0.5. Для первого и последнего уравнений системы (6) b1u1  a1u2  d1 , cN uN 1  bN uN  dN , В зависимости от типа заданных граничных условий значения коэффициентов в этих уравнениях представлены в таблице ГУ 1-го ГУ 2-го ГУ 3-го a1 -1 -1 b1 1 1 (1+ Bi h) d1 φ0(τ) S0(τ)h Bi h θcp1 сN 1 1 bN 1 -1 -(1+Bi h) В таблице используются обозначения 0 ( )  (T0 ( ts ) Tc ) / Ts , 1( )  (T1( ts ) Tc ) / Ts S 0 ( )  q 0 ( ts )l , Ts S1 ( )  cp1  (Tcp1 Tc ) / Ts , Bi   l / . q1 ( ts )l , Ts cp2  (Tcp2 Tc ) / Ts , dN φ1(τ) S1(τ)h -Bi h θcp2   2 Повышение точности численной аппроксимации уравнения  2 0   Для аппроксимации производной по времени использовали выражение  ui , j 1  ui , j  , которое имеет первый порядок точности по шагу сетки k=∆τ.  k i ui 1, j  2ui , j  ui 1, j  2 Явная схема (1) j+1  2 2   i h   j j j+1 i-1 i-1 i i i+1 i+1 i j+1 i  i+1 j i  i   j 1   2  i-1 2 ui 1, j 1  2ui , j 1  ui 1, j 1 h2 Аппроксимации второй производной по ξ имеет второй порядок точности по шагу сетки h= ∆ξ. Вычислим вторую пространственную производную для момента времени    j    . Получим j i-1  2 Неявная схема (2) i+1 2  i   j  k  (1   )  ui 1, j  2ui , j  ui 1, j h2 ui 1, j 1  2ui , j 1  ui 1, j 1 h2 (3) Неявная схема (метод) Кранка – Никольсона Схема (3) при ε=0 дает явную схему (1), при ε=1 – неявную схему (2). При ε=0,5 получаем сему второго порядка точности по временному шагу . В этом случае для точек сетки i=2,…,N-1 получим уравнения ui 1, j 1  (2  p )ui , j 1  ui 1, j 1  pui , j  (ui 1, j  2ui , j  ui 1, j ), (4) или в общей записи ci ui 1, j 1  bi ui , j 1  ai ui 1, j 1  di где p= 2h2/k Соответственно значения коэффициентов используемых в программе реализующей решение системы уравнений (4) методом прогонки c2  c3    cN 1  1, a 2  a3    aN 1  1, b2  b3    bN 1  (2  p ), di  pui , j  (ui 1, j  2ui , j  ui 1, j ), i  2,..., N  1. Повышение порядка точности аппроксимации граничных условий Граничные условия 2-го и 3-го рода включают первую пространственную производную искомой функции θ . В частности для границы ξ=0 имеем  S ( )   ,   Bi(cp    0 ) (5)   0 h x  1 Очевидно использование аппроксимации первой производной через значения сеточной функции типа u u   2 1,   0 h  имеющей первый порядок точности по шагу сетки h не позволит повысить общую аппроксимацию решаемой задачи и в случае использовании схемы Кранка – Николсона. Как повысить точность? Запишем разложение θ в ряд Тейлора на границе ξ=0  ( , ) 1  2 ( , )  ( , )   (0, )      0 2  2 ( ) 2  ...,  0 (6) но согласно исходному уравнению должно выполнятся равенство  2 ( , )  2   0  ( , ) .   0 (7) Пренебрегая членами разложения выше второго порядка выразим из (6) первую пространственную производную θ с учетом (7). Получим  ( , )  ( , )   (0, ) 1  ( , )    .    2    0  0 Через значения сеточной функции выражение (8) на временном слое j+1 запишется так u2, j 1  u1, j 1 h (u1, j 1  u1, j )  ( , )   ,  h 2 k  0 Соответственно граничные условия (5) на границе ξ=0 примут вид u1, j 1 (1  p / 4)  u2, j 1  S ( )  pu1, j / 4, u1, j 1 (1  p / 4)  u2, j 1  Bi h (cp  u1, j 1 )  pu1, j / 4. (8) (9) Аналогичные выкладки можно проделать и для ГУ 2-го и 3-го рода заданных на границе ξ=0. Коэффициенты для первого и последнего уравнений системы (4) b1u1  a1u2  d1 , cN uN 1  bN uN  dN . ГУ 2-го рода b1  (1  p / 4), a1  1, d1  S 0 ( )h  u1, j p / 4; cN  1, bN  (1  p / 4), dN  S1( )h  uN , j p / 4; ГУ 3-го рода b1  (1  Bi h  p / 4), a1  1, d1  Bi h cp1  pu1, j / 4 cN  1, bN  (1  Bi h  p / 4), dN  Bi h cp2  puN , j / 4