Диссертация (Моделирование пространственных течений в газовых трактах с использованием адаптивных сеток), страница 8

величина y   30,400 , чтодостигаетсялокальнымперестроениемсетки(измельчениемилиукрупнением).3.1.3. Построение численной схемы повышенного порядка точностиРешение системы дифференциальных уравнений, моделирующихтечение невязкого газа, осуществляется при помощи метода типа Годунова. Вметоде, предложенном С.К. Годуновым, при решении системы уравненийдля вычисления потоков используются значения, получающиеся прирешении газодинамической задачи о распаде произвольного разрыва сначальными условиями: , x  0,0   1 2 , x  0(3.49)где 1 и  2 - значения в соседних ячейках.Численная схема первого порядка точности для пространственногослучая, основанная на методе Годунова, имеет вид:Vii  Fnj S j i  0, i  1,.., N ,tj(3.50)где N – число контрольных объёмов, i - среднее по контрольномуобъёму значение консервативной величины, Fn – поток, нормальный граниячейки, вычисляемый из «больших переменных», полученных решениемзадачи Римана на грани, Sj – площадь грани, Vi – объём ячейки. jn 1   jn   , t  F , n dx  G, n dyS(3.51)где F – вектор потоков, n – вектор нормали к ячейке.

Согласно методуГодунова для расчёта величин на гранях ячейки, необходимых для расчёта56потоков через грани, используется решения задачи распада разрыва (задачаРимана).Согласноработечисленный[72]потокf ABчерезграньаппроксимируется со вторым порядком точности формулой:f AB  f  x1 , y1 , t y AB  g  x1 , y1 , t x AB .(3.52)Для аппроксимации потока через грань с порядком точности большевторого необходимо учитывать больше число членов в разложении,например для третьего порядка точности необходимо второй член, K=2:f AB 1 f  x1 , y1 y AB  g  x1 , y1 x AB  21 f  x2 , y 2 y AB  g  x2 , y 2 x AB 2x  xBy  yBx1  A x AB / 2, y1  A y AB / 2 .22x  xBy  yBx2  A x AB / 2, y 2  A y AB / 222  1/ 3(3.53)Таким образом, для повышения порядка точности расчётной схемы попространству до второго достаточно изменить представление данных внутриячейки, не меняя аппроксимацию интеграла вдоль контура.

Такое повышениеточности достигается при использовании кусочно-линейной аппроксимацииискомых величин внутри ячейки [64], в отличие от кусочно-постоянного ихпредставления в оригинальном методе Годунова [12].Запишем разложение функции в ряд Тейлора второго порядка точностив окрестности точки центра масс ячейки: x    x0   x  x0  x0   O x  x02.(3.54)Таким образом, для интерполяции значений функции внутри ячейкивторого порядка точности можно задать внутри ячейки распределение вида:57 t n , x   mn  mn r  rm  , где  mn - градиент функции в центре масс ячейкис порядком точности не меньше первого.Для аппроксимации градиента можно воспользоваться формулами,описанными выше, однако в этом случае возможно возникновениенефизических осцилляций. На практике важным свойством численных схемявляется способность получать численные решения не осциллирующиевблизи разрывов решения. Таким свойством обладают монотонные схемы.Схема называется монотонной, если выполнено условие:Hi k ,...,i k   0, j : j  k  j  i  k , j(3.55)H i k ,...,i k  - оператор перехода на новый временной слой.Годунов (11) доказал утверждение, что линейные схемы на постоянномшаблоне, обладающие свойством монотонности, не могут иметь порядокточности выше первого.

Соответственно, для повышения порядка точностирасчётной схемы она должна обладать определённой нелинейностью.Для обеспечения монотонности решения зададим интерполяциюпараметров внутри ячейки в виде: t n , x   mn  mn r  rm  ,(3.56)где λ – нелинейный коэффициент, позволяющий уменьшать влияниеградиентов вблизи областей сильного изменения величин, тем самымуменьшая порядок схемы до первого. В этом случае численная схема,полученная модификацией схемы Годунова с использованием кусочнолинейных распределений, будет иметь второй порядок точности на гладкихрешениях и первый в областях сильных разрывов.58Для определения параметра λ используются функции-ограничители.Наиболее часто используемой является функция minmod, описанная в работе(12):  min moda, b 1signa   signbmin  a , b .2(3.57)На практике также часто используются следующие функции:min mod a,2b , a  bsup erbeea, b   ; min mod2a, b , a  b(3.58)vanleer a, b  (3.59)2ab.abТакже для определения величиныf  qmn возможно использоватьметод основанный на тесселяции Фонга.

Основная идея тесселляции Фонгадля треугольника состоит в следующем: через каждую вершину треугольнойячейки vi проводится касательная плоскость, и точка P' проектируется накаждую из трех получившихся плоскостей. После этого получившиеся точкиумножаются на веса u, v, w и складываются давая искомое положениевершины, соответствующее данным барицентрическим координатам.   P*  u 0 P '  v1 P '  w 2 P ' ,(3.60)где u, v, w находятся как координаты точки P* в барицентрическойсистеме координат треугольной ячейки.Нормаль в этой точке находится по следующей формуле:n* un0  vn1  wn2.un0  vn1  wn2(3.61)Для центра масс ячейки барицентрические координаты центра массячейки u=v=w=1/2.59При использовании алгоритма Фонга для определения линейногораспределения величины q внутри ячейки рассматриваются три плоскости,проходящие через центр ячейки и центры соседних ячеек.

В итоге получаетсянесколько градиентов величины для одной точки. Для ограничения величиныградиента можно выбрать градиент, обладающей наименьшей нормой, нолучший результат получается в результате использовании суперпозиции всехтрёх градиентов. Наиболее простой способ предложен в работе (13). Приэтом способе выбирается такая норма градиента, что значение, полученноепри интерполировании в центры соседних ячеек, не будет превосходитьзначения в этих ячейках. Это достигается при вычислении величины φ поформуле: maxqi  q0 ,0 ,rrqi00 minqi  q0 ,0 φ, ri  r0 q0 1ri  r0 q0 maxqi  q0 ,0 ri  r0 q0 minqi  q0 ,0  .(3.62)Приведённая выше модификация схемы Годунова, основанная нааппроксимации кусочно-линейными функциями параметров внутри ячейки,относится к TVD-классу, т.е.

классу схем с невозрастающей полнойвариацией:TVn 1  TVnNTVn   U in1  u in .(3.63)i 1Оригинальный метод Годунова для решения уравнений газовойдинамики имеет первый порядок точности по времени. Повышения порядкаточности схемы по времени производится при помощи метода Рунге-Кутты.Схема третьего порядка точности по времени имеет вид:60q 0   q n  q 1  q 0   t F q 0 q 2  .(3.64)3 0  1 1 tq  q F q 1444  В общем виде N-стадийный алгоритм представляется в виде:Q 0   Q  n Q m   Q 0   m tR Q m 1 , m  1,..., N ,Q n 1Q(3.65)N где λ – коэффициенты Рунге-Кутта.Шаг по времени выбирается на основании критерия КурантаФридриха-Леви: Lt  Ku min   ii vin  ci ,(3.66)где i – номер ячейки, vin - нормальная к грани ячейки компонентаскорости в i-й ячейке,ci - скорость звука в i-й ячейке, Li - расстояние между центром ячейкии ближайшим ребром, Ku – число Куранта (0,1).3.1.4. Учёт изменения теплофизических параметров газаВ общем случае термодинамические параметры воздуха являютсяфункцией двух термодинамических величин, но для давлений меньших1МПа в качестве уравнения состояния можно взять уравнение МенделееваКлапейрона, а все термодинамические параметры считать функцией толькотемпературы.

При таком приближении погрешность вычисления показателяадиабаты воздуха не превышает 3% (на основе сравнения данных дляреальных газов из [69], [70]). При вычислении удельных теплоёмкостей припостоянном давлении и объёме погрешность не превышает 1.5% .61Реальнаятермодинамикагазамоделируетсяспомощьюполиномиальной интерполяции приведённого потенциала Гиббса [69].Алгоритм точного решения задачи Римана приводится в [68]. Самойресурсоёмкой частью решения является определение интегралов Римана,поэтому они вычисляются приближённо с использованием кубическойинтерполяции изоэнтропы [77]. Отличие используемого в работе метода от[77] заключается в предварительной табуляции значений интегралов оттемпературы с целью уменьшения интервалов интегрирования.Важнымвопросомявляетсянеобходимостьучётаизменениятеплофизических параметров газа при решении задачи Римана о распадепроизвольного разрыва. Для оценки влияния изменения теплофизическихсвойств на точность решения задачи Римана рассмотрим следующиетестовые задачи (Таблица №1).В тесте №1 происходит столкновение двух сверхзвуковых потоков, врезультате которого возникают две ударные волны.

Тест №2 отличается оттеста №1 различными температурами слева и справа от контактного разрыва.В тесте №3 образуется волна разряжения, контактный разрыв и ударнаяволна. Задача Римана о распаде произвольного разрыва ставится наинтервалеx   0.5,0.5, t  0, t  .Параметры слева и справа от точкиконтактного разрыва xo   0.5,0.5 в начальный момент времени задаютсяпостоянными и равными  l , ul , pl  на интервале x   0.5, xo  и  r , u r , pr  наинтервале x  xo , 0.5 .№123 l , кг/м³ u l , м/с1.20.51.0100010000p l , Па·10511010 r , кг/м³ u r , м/с1.22.00.125-1000-10000p r , Па·105 t , мc x0, м1100.110.30.2000Таблица №1. Тестовые задачи.Решение данных задач производится в предположениях: газа спостоянным показателем адиабаты во всём пространстве(1), газов с62различными показателями адиабаты по обеим сторонам начального разрыва,но постоянными в каждом полупространстве (2), точное решение для газа сучётом зависимости показателя адиабаты от температруры газа(3).Решение газодинамической задачи Римана для теста №1 состоит издвух ударных волн.

Поскольку температуры газа слева и справа отконтактного разрыва в начальный момент равны, решения, полученные спомощью методов (1) и (2) совпадают. Погрешность в вычислении скоростейударных волн составляет 6%, погрешности плотности, температуры идавления в точке x=x0 составляют 5%, 6% и 1% соответственно.Тест №2 отличается от теста №1 различными первоначальнымитемпературами с разных сторон от контактного разрыва. Структура течениясостоит из двух ударных волн и контактного разрыва. Погрешностиопределения скоростей ударных волн для методов (1) и (2) составляют 5% и4% соответственно.

В таблице №2 представлены погрешности в определенииплотности и давления в зонах между левой ударной волной и контактнымразрывом (УВ-КР), и контактным разрывом и правой ударной волной (КРУВ).Метод (1)УВ-КРКР-УВε(ρ), %92.7Метод (2)ε(p), %1.61.6ε(ρ), %7<0.1ε(p), %1.21.2Таблица №2. Погрешности определения параметров в тесте №2.В результате решения газодинамической задачи теста №3 образуетсяволна разряжения, контактный разрыв и ударная волна. Погрешностипредставлены в таблице №3.Метод (1)Метод (2)ε(ρ), % ε(u), % ε(p), % ε(ρ), % ε(u), % ε(p), %ВР x=-0.22.57.63.40.01.70.0ВР x=-0.10.12.33.70.10.20.2ВР x=01.80.64.70.10.20.5ВР-КР x=0.13.11.74.42.11.61.9КР-УВ x=0.256.21.74.43.11.61.963Таблица №3.

Погрешности определения параметров при решении задачи №3.Для определения влияния точности решения задачи Римана наточность решения газодинамической задачи в целом рассмотрим две задачи:обтекания сверхзвуковым потоком уступа в полу бесконечном пространствеи течения газа в плоском канале с сужением.Вид расчётной области для решения задачи обтекания уступа приведённа рисунке 3.1.