книги / Основы математического моделирования и численные методы
..pdfЛ т " |
(*‘ ) + Т \ т " Ы |
+ 0 ( ь 7 ) - 2 Т Ы + |
+ Т(хк)-ЬТ'(хк) Л г ( х к) - |
||
- ^ Т' ( х‘ ) +^ |
Т'Г(** ) - ^ тГ |
(** ) + о(А7 ) |
|
J _ |
|
|
h2 ^ г ( *‘ )+ 7 Г г “’ (* * )+ 0 (А‘ ) |
|
т; - Г Ы -Л т»{хк)+0(ъ<)--0{*).
Видно, что порядок аппроксимации в этом случае равен двум.
2.5. Построение разностных схем
Рассмотрим одномерное уравнение теплопроводности
Э77дт=Э277Эх2 |
(2.12) |
в области G = {0< x< d,0< x< t}. Пусть функция температуры Т(х)
и ее производные до 4-го порядка по переменной х и до 2-го порядка по времени т непрерывны. Введем пространственно-временную сетку
=-kh , k = 0,N,h =d/N ;т, = tr\,i = 0,M,r\ = t/M } .
Для производной по времени выберем правое разностное соот ношение в точке (/, к):
т;=(т1"-т;)/п. (2.13)
Для второй производной по координате х на г-м временном слое запишем разностное соотношение (2.11):
2^.* = (г;., - 2 2 ;' + |
) / А2. |
(2.14) |
на временном слое i + 1
= ( С - 2 7 f ' + r;:l ) / А2 |
(2.15) |
Заменяя в исходном дифференциальном уравнении (2.12) про изводные выражениями (2.13) и (2.14), получим явную разностную схему; используя выражения (2.13) и (2.15) - неявную разностную схему:
(П*}- п ) _ ( т ^ - 2 т ;+ т Ц
(2.16)
Ч А2
(2.17)
T] |
h2 |
Используя среднее значение для аппроксимации второй произ водной относительно уравнений (2.16) и (2.17), получим еще одно разностное уравнение, называемое схемой Кранка - Николсона на шеститочечном шаблоне:
. (2.18)
Л |
2 |
А2 |
+ |
А2 |
Шаблоны разностных схем приведены на рис. 2.8.
/Л |
|
|
|
/А |
» |
* |
|
J+ 1 |
|
|
|
7 + |
1 |
||
|
|
|
|
|
|||
j |
|
|
|
j |
|
|
|
i - 1 |
----------------- i |
/+1 |
> |
|
|
/ - 1 |
-------- > |
x |
|
|
/ + 1 x |
||||
Рис. 2.8. Шаблоны разностных схем
Разностная схема (2.16) называется явной, поскольку неизвест ная функция Т‘к х на (/ + 1)-м временном слое может быть явно вы ражена через известные функции Т‘к на /-м временном слое. При ис пользовании схемы (2.17) явно выразить неизвестные невозможно и для вычисления Т'к х необходимо решить систему N - 1 линейных алгебраических уравнений на каждом временном слое. Такие разно стные схемы называют неявными. Таким образом, все разностные схемы по этому принципу (как выражены неизвестные) делятся на явные и неявные. Достоинством явных разностных схем является простота вычислений. Однако явные схемы условно устойчивы, т.е. существует условие, при выполнении которого решение дискретной задачи может быть получено, при этом достигается сходимость и ус тойчивость метода. Такое условие накладывает ограничение на вы бор шагов по времени и координате. Иначе говоря, не всякое соот ношение этих шагов приведет к реальному, адекватному результату решения задачи.
Каждому виду уравнений соответствует, если это возможно, свое условие сходимости вычислительного процесса. Для уравнения
дТ I Эт = ад2Т / дх2 при использовании явной разностной схемы по лучено условие
i\ =h2/ 2 l a . |
(2.19) |
Условие (2.19) накладывает существенное ограничение на шаг по времени.
Неявные схемы (2.17), (2.18) являются безусловно устойчивы ми, на них не накладывается какого-то ни было ограничения по ша гам сетки (возможно любое соотношение шагов). Это достоинство неявных методов. Недостатком является необходимость решения систем алгебраических уравнений с большим числом неизвестных на каждом временном слое.
Погрешность аппроксимации для первой и второй схемы -
O (TI + /*2), для схемы Кранка - Николсона —О(ri2 +й2), т.е. в по
следнем случае исходное дифференциальное уравнение аппрокси мируется со 2 -м порядком точности по времени Г|.
В общем разностные схемы могут быть получены на основе ли нейных комбинаций уравнений (2.16), (2.17).
2.6. Интегроинтерполяционныи метод, метод баланса
Для численной реализации краевых задач возникает необходи мость при построении разностной схемы соблюдать законы сохра нения (энергии, массы, количества движения и т.п.), на которых ос нована постановка задачи в дифференциальной форме (обычно дифференциальное уравнение отражает какой-либо физический за кон сохранения). Разностные схемы, удовлетворяющие этим требо ваниям, называются консервативными; схемы, которые не следуют законам сохранения, - неконсервативными. Одним из наиболее эффективных методов построения консервативных разностных схем является интегроинтерполяционный метод [2 ].
Таким образом, получение консервативной разностной схемы исходит из уравнения баланса, записанного для элементарного объ ема (ячейки) сеточной области. Производные и интегралы в уравне ниях баланса заменятся приближенными разностными соотноше ниями на ячейках сетки. В результате записывается однородная раз ностная схема.
Рассмотрим уравнение теплопроводности с переменным коэф фициентом на отрезке [0 , 1 ]:
(Хи\х))' = 0 , м(0 ) = 1, «(1) = 0 . |
(2 .2 0 ) |
Продифференцировав уравнение (2.20), получим
Хи +Х'и = 0 .
Если применить метод формальной замены производной, полу чим следующую разностную схему:
и ., - и.
|
fy+l |
^1-1 | ^l+l |
^*1-1 ty+1 Щ-\ _ Q. |
(2.21) |
|
' |
h2 |
2A |
2h |
||
|
0 < i< N , щ = 1, % = 0.
При использовании неконсервативной разностной схемы (2.21) решение расходится и результат получается неверным. С точки зрения вычислительного процесса происходит накопление ошибки, с физической точки зрения - не выполняется закон сохранения энергии (возникают дополнительные источники притока или стока тепла).
Рассмотрим метод баланса для стационарного уравнения теп лопроводности на том же отрезке [0, 1]:
^ =qv(x) +k(x)T, 7X0) = 7], Т(1) =Т2, к(х)> 0, À(x)>0; (2.22)
ах
q = -X{x)— =- X r |
(2.23) |
dx |
|
Введем равномерную сетку на отрезке с шагом h. Определим баланс потоков тепла на отрезке xf_0 5 < х < xi+0 5, проинтегрировав
уравнение (2.23) в пределах этого отрезка:
|
*/+0,5 |
*/+0,5 |
|
“ й+0.5 +?,-о.5 |
= J <?v(*)<&+ { k(x)T(x)dx, |
(2.24) |
|
|
*/-0,5 |
*/-0,5 |
|
где #;_05 - количество |
тепла, |
подводимое к сечению |
JC;_05 на |
отрезке длиной h; qi+05 - количество тепла, отводимое в сечении
*/+0,5
х|+05 этого отрезка; J qv(x)dx - количество тепла выделившегося
*/-0,5
на отрезке х,_05 < х < х,+0 5 за счет внутренних источников тепла;
*i+0,5
J k(x)T(x)dx - количество тепла, отведенного во внешнюю
* 1- 0,5
среду.
Предположим, что на каждом отрезке [x,_0j5,xJ+05] значение температуры постоянно, т.е. Г = const = Tt , тогда
■*1+0,5 1 -*i+0,5
iГ *(х)Г(х>й: = /ВД., dt = - |
Jf k(x)dx, |
*i-0,5 |
*i-0,5 |
где à, - среднее интегральное значение к(х) на отрезке длиной И.
Проинтегрировав выражение для потока тепла (2.23) на отрезке
х;], |
получим |
|
|
|
|
|
|
|
|
jrp |
*/\dT=-*i г |
d x . |
(2.25) |
||
|
|
X dx |
■* |
|
■*... M*) |
|
|
Предположим, что на отрезке |
[ |
x , |
поток тепла неизменен |
||||
и равен значению в центральном узле отрезка - |
q = const = qf_o s . |
||||||
Тогда, взяв интеграл левой части (2.25), будем иметь |
|
||||||
|
т _ т |
г |
dx |
_ _ |
п |
(2.26) |
|
|
h |
Vi - ^-0.5 J Т7Т“ |
^-0.5- |
||||
|
|
|
* 1-1 |
Ux) |
|
|
|
|
п-1 |
|
|
|
|
|
|
где а, = |
1 % dx |
. Определим из формулы (2.26) #,_0 5 : |
|||||
L |
1М д |
|
|
|
|
|
|
|
|
4 i - 0.5 |
= - « , |
|
|
|
(2.27) |
Рассмотрев отрезок [х, ,х/+1] и проинтегрировав на нем выраже ние (2.23), получим
Поток тепла q = const = qi+05 на отрезке [х(,дг,+|]. Аналогично
предыдущим действиям будем иметь
нение (2.24) и получим консервативную разностную схему
Нужно отметить, что одним из достоинств метода баланса яв ляется возможность его использования для уравнений с разрывными коэффициентами, поскольку интегральная интерпретация записи за кона сохранении (в частности, энергии) позволяет получить адек ватное решение.
2.7. Устойчивость и сходимость разностных схем
Найдем решение задачи теплопроводности в стержне единич ной длины, описываемой следующим уравнением и краевыми усло виями:
Э77 Эт = Э2Г / Эх2,0 < т < 1,0 < х < 1 ;
Г(х,0) = 0,
[дТ /дх]х=0= 0,
[9 г / Н - . =1-
Введем |
прямоугольную |
сетку: |
xk =kh,(k |
= 0,1,—,Ю> |
т, =гт],(/ = 0,1,...,М), где h =\ / N,r\ =\ / М. |
Используя |
различные |
||
варианты записи разностного соотношения для производной по вре мени, получим три разностные схемы:
- вариант 1 :
(г;*1- п ) / ч =щ (гы - гтк + ii_ ,)i h2
|
(2.29) |
(к = 1,2, |
1, i = О,1 , М -1); |
вариант 2:
(г ;-2 -;-1) / л = а » ( ^ | - 2 Гк +т;_,)/н2
(2.30)
(* = 1,2.....Д'-l, 1 = 1,2
вариант 3:
{ П " - T f ' ) / (2Т1)= а, (Гы - 2Гк + 2 t , ) / h 2
(2.31)
(к = 1 , 2 , N - 1 , 1 = 1,2,..., М -1).
Ккаждому из выражений (2.29)-(2.31) добавляются начальные
играничные условия, записанные в разностном виде:
тк =о>(г ; - г„')/а = о,
( ^ - r ; . , ) / A = 1, 0 = 1,2
Можно показать, что погрешность аппроксимации исходного дифференциального уравнения разностными схемами (2.29) и (2.30)
определяется как o[h2 +Г|), а схема (2.31) - 0 (/г2 +т|2). Поскольку
аппроксимация краевых условий имеет 1-й порядок относительно Л, все варианты разностных схем имеют один и тот же порядок ап проксимации исходной задачи по времени.
Первая и третья схемы представляют собой явные разностные схемы. Решение задачи по схеме 3 начинается с вычисления значе-
ний неизвестной на первом временном слое через известные значе ния на нулевом (в начальный момент времени), например с помо щью разложения в ряд Тейлора:
Г{=7*+л(Э Г/Э т)2+0(||г) (*-1,2,..,Л Г -1), |
(2.32) |
при этом
( т Х = а°лтЛ - |
<2-зз> |
Подставив формулу (2.32) в формулу (2.33), определим Т \. Да
лее вычисления такие же, как и для схемы 1. Варианты 1 и 3 с точки зрения вычислительного процесса наиболее просты.
Казалось бы, что использование варианта 3 наиболее целесооб разно, поскольку исходное уравнение аппроксимировано по сравне нию с вариантами 1 и 2 более точно. Вычислим поля значений температуры, используя три варианта решения. Для анализа полу
ченных решений введем величину (3 = Г|/й2, определяющую соот
ношение шагов сетки. Расчет поля температур исходной задачи по трем вариантам (2.29)-(2.31) показал, что второй вариант позволяет получить удовлетворительные результаты, а расчеты по вариантам 1 и 3 приводят к абсурдным результатам при (3 = 1 , т.е. h = 0 ,1;
ri = 0,01 уже на пятом временном слое т5 =0,05. В случае явных разностных схем 1 и 3 к такому результату приводит так называемая неустойчивость разностных схем, когда в процессе счета вычисли тельные погрешности, связанные с округлением чисел или их огра ничением, быстро растут, приводя к нереальным результатам. Такие схемы называются неустойчивыми.
Если же возникающие в процессе расчетов вычислительные по грешности имеют тенденцию убывать (по крайней мере, не возрас тают), то разностная схема называется устойчивой, т.е. малому из менению входных данных соответствует малое изменение решения.
Уменьшение шага по времени |3 = 1/2, TJ = 0,005 позволило
решить задачу и получить близкие значения по температуре для ва
риантов 1 и 2. Поскольку решение по варианту 2 получено при усло вии выполнения соотношения шагов Р = 1/ 2иг| = 0,005, то исполь зованная в этом случае схема является условно устойчивой.
Разностные схемы, которые при достаточно мелких шагах по независимым переменным и любом сколь угодно малом Р > 0 неус тойчивы (вариант 3), называются абсолютно неустойчивыми.
Если же факт устойчивости имеет место при любом соотноше нии шагов по различным переменным, лишь бы они были достаточ но малы, то схема называется безусловно устойчивой (вариант 2).
В общем случае разностное решение или его устойчивость мо жет зависеть от правой части, граничных и начальных условий. Исхо дя из этого различают три вида устойчивости: по правой части, гра ничным условиям и начальному условию. Для исследования устойчи вости разностной схемы используют различные методы: разделение переменных, принцип максимума, разложение в ряд Фурье и т.д.
2.8. Экономичные схемы
Несмотря на то, что при решении нестационарных краевых за дач в двухмерных и трехмерных постановках построение разностно го аналога не вызывает затруднений, возникают сложности в реали зации, так как существенно возрастает количество неизвестных в разностных уравнениях. Пусть для всех пространственных перемен ных выбран шаг h, тогда на каждом временном шаге необходимо решать систему алгебраических уравнений с ~КР неизвестными, где р - число измерений. В связи с этим значительно возрастает объем физических операций, необходимых для решения системы разност
ных уравнений. Так, если в прямоугольнике (0<х, <d, 0 < х 2 <è)
задана краевая задача с граничными условиями первого рода, |
|
||
дТ |
д*Т | |
i2«Л |
|
э т |
|
||
Эт = а |
иЛ|2 |
Э*2 J ср |
|
0<*,<*/, |
0 < дt2 < b, 0 < т < Г ; |
(2.34) |
|