Н.Н. Калиткин - Численные методы (1133437), страница 94
Текст из файла (страница 94)
Уравнение (48) выражает закон сохранения импульса и во всех численных расчетах используется именно в такой форме. Уравнение изменения энергии (49) обычно удобнее преобразовать. Если подставить в него Йуе, определенную из уравнения (47), то получим особенно простую форму записи: — 'уе+)т-'-('-) =О дт вг , р г (50а) Если умножить уравнение (48) на о и прибавить к уравнению (49), то получим другую форму — закон сохранения полной энергии: р лг( + о'о )+п)ч(Ро) =0 (50б) В том, что указанные уравнения являются законами сохранения, нетрудно убедиться. Например, проинтегрируем (боб) по объему д)г, занятому некоторой массой вещества, второй интеграл преобразуем к поверхностному, а в первом интеграле заменим р о)г на г)т, после чего производную по времени можно вынести за знак интеграла.
Тогда получим — ~ (е+-- оэ) г)т+ с~ рв гВ=О. Здесь первый интеграл есть полная (внутрещгяя и кинетическая) энергия данной массы газа; второй равен работе в единицу времени сил давления на поверхности, ограничивающей данную массу газа. Действительно, это закон сохранения энергии. Одномерные задачи бывают трех типов: с плоской, цилиндрической или сферической симметрией. Введем показатель симметрии у, равный для этих случаев соответственно О, 1 или 2. Масса слоя толщиной г(» в этих случаях равна (рис.
97) г(пт = рг' г(г (51) Уравнение неразрывности (47) выражает закон сохранения массы. Его тоже обычно преобразуют, но уже после приведения урав- нений к одномерной записи. 441 ОДНОМЕРНЫЕ УРАВНЕНИЯ ГАЗОДИНАМИКИ с точностью до численного множителя, равного 1, 2п или 4п. При помощи соотноп1ения (51) введем массовую координату данной материальной точки: с т(г) =~ р(9)$'с(9. (52) О По закону сохранения вещества массовая координата материальной точки не меняется со временем; поэтому такая координата позволяет легко следить за каждой частицей вещества и, в частности, за границей раздела слоев. Преобразуем уравнения газодинамики в одномерном случае к лагранжевой форме.
В качестве первого уравнения возьмем определенье скорости: дс — = О. дг 'Рл (53) Уравнение неразрывности (47), выражаРис. 97. ющее закон сохранения массы, заменим имеющим тот же смысл соотношением (51), записывая следующим образом: д „„а У+1 дт р (54) В уравнениях импульса (48) и энергии (50б) перейдем к производной по массовой координате, то есть в одномерных выражениях д игад = —, дг Йу=г-' — (г' ...) дг заменим дг на дт при помощи соотношения (51), и получим дй „дР др + г~с-1 =О (55) — 1е+ — О') + — (г'Ор) = О.
дГ~ 2 ) дт (56а) Уравнение энергии можно взять также в форме (50а): (56б) Система (53) — (56) является лагранжевой формой записи уравне- ний одномерной газодинамики. В большинстве численных расче- тов используется эта именно форма. 442 гнпвгволнчвскив тгхвнвния 1гл. хги 2. Псевдовязкость. Уравнения (53) — (56) составляют гиперболическую квазилинейную систему. Из курса газодинамики известно, что среди ее решений есть сильные разрывы — ударные волны. В главе 1Х мы видели, что для разностного расчета таких решений надо изменять уравнения, вводя в них специально подобранные днссипативные члены. В газодинамнке такие члены удается найти из физических соображений. Дело в том, что уравнение газодинамики сравнительно грубо описывают поведение газа. Эти уравнения выводятся из кинетического уравнения Больцмана для функции рас« пределения молекул.
Если при выводе учесть эффекты диффузии молекул, то в уравнениях газодинамики появятся так называемые вязкие члены. Например, уравнение импульса (48) примет вид (см. 1191) р -;- = — дгас1 р+т) Лп+~йгас) с)1ч п, ~) -8 т) О, (57) дэ 1 где т1 и Ь являются коэффициентами физической вязкости. Учет физической вязкости приводит к изменению качественного характера решения: плоские ударные волны превращаются в аналитические решения, в которых скачки сглажены и имеют эффективную ширину порядка длины свободного пробега молекул.
Качественно это легко понять на примере плоского течения, где уравнение (57) принимает форму 1 д7 (1+ -)дхт дх ' да Уо др (58) напоминающую уравнение теплопроводности; видно, что вязкий член должен сглаживать разрывы решения. Обычно в численных расчетах довольствуются только вторым вязким членом в уравнении (57) и считают коэффициент с" ,слабо меняющимся. Тогда этот член можно объединить с давлением: — ягас( р+ ь ягас) 61ч и -- — ягас1 (р — ь с)1ч е), (59) и рассматривать величину спс = — ь с(1ч и (60) 81 (т-' 1) и 8 (6!) как вязкое давление.
При этом в уравнение энергии вместо обычного давления также ставят величину р+сп„называя ее полным давлением. Вязкость са, называется линейной. Она приводит к «размазыванию» ударной волны со скачком скорости бо на интервал ОДНОМЕРНЫЕ УРАВНЕНИЯ ГАЗОДИНАМИКИ 443 где у — показатель политропы вещества. Физический коэффициент вязкости очень мал и дает ничтожно малое сглаживание. Для расчетов по разностным схемам необходимо сглаживание на несколько интервалов сетки.
Поэтому в численных расчетах величину Ь приходится увеличивать на много порядков по сравнению с ее физическим значением. Для численных расчетов необходимо введение вязкости лишь в окрестности ударных волн, Но вязкий член в (69) присутствует во всех точках пространства. Поскольку в численных расчетах коэффициент ь выбирается много больше физического коэффициента, то наличие псевдовязкостн, помимо положительного эффекта — сглаживания разрывов, — приводит к отрицательному— внесению заметной погрешности. Чтобы уменьшить эту погрешность, Нейман и Рихтмайер [72) предложили выбирать коэффициент псевдовязкости большим в окрестности скачков скорости бо н малым в зонах гладких течений, где скорости соседних точек близки. Для этого они положили ь=с,(т((ч (62) где ьв — коэффициент, небольшой по величине.
Такая псевдовязкость называется квадратичной, потому что она приводит к вязкому давлению: таз = — Ьа (5!йп п!ч е1) (Й!ч т1)а. Переписывая (62) в виде ь ьо(бо/бг) и подставляя в (61), нетрудно убедиться, что квадратичная псевдовязкость сглаживает скачок бо любой интенсивности на один и тот же интервал: бг 1у —. аьо У (у+1)р (64) Обычно коэффициент псевдовязкости ье выбирают так, чтобы бг равнялось 2 — 3 шагам разпостной сетки. Линейная вязкость приводит к монотонным (нли почти монотонным) разностным решениям, так как ей соответствуют аналитические точные решения, которые хорошо аппроксимнруются разностными схемами; зато фронты скачков при этом сильно сглажены.
Квадратичная вязкость приводит к более крутым фронтам; но ей соответствуют точные решения с разрывами первой или второй производной, поэтому разностное решение немонотонно вблизи слабых и сильных разрывов. Нередко используют комби- В главе Х, 5 2 было проведено строгое исследование квадратичной псевдо- вязкости иа примере простейшего квазилинейного уравнения (10.44); при этом для зоны сглаживания сильного разрыва было получено выражение (10.51)— (10,о2), сходное с (64).
444 ГИПЕРБОЛИЧЕСКИЕ УРАВИЕИИЯ 1гл. хи1 нацию линейной и квадратичной вязкости !и==-а,!и!+и!сил с экспериментально подобранными коэффициентами. Поскольку б(у и = — р 1(др,'дг), то вязкое давление положительно при сжатии и отрицательно при разрежении вещества. Сильными разрывами являются только ударные волны, поэтому для сглаживания разрывов можно вводить вязкость только при сжатиях. При (др!д!) (О присутствие псевдовязкости не обязательно и даже уменьшает точность расчета. Поэтому обычно по- лагают — ь1д!Уп+сл(61чп)' при 61ча(О, Ы~ =- (65) О л при ЙУ и ==- О. Этот вид псевдовязкости независимо предложен рядом авторов (см.
(27!). 3. Схема «крест». Это наиболее простая н довольно точная однородная разностная схема счета газодинамики. Ее шаблон приведен иа рис. 96; значения радиусов приписываются узлам сетки, значения скорости — границам пространственных интервалов на полуцелых слоях, а значения плотности, давления и внутренней энергии — серединам интервалов на целых слоях. Построение схемы напоминает акустический лкресть. Для простоты записи выберем равномерные по массе и времени шаги л! и т и аппроксимируем систему (53)— (56б) следующими разностными уравнениями: л л л l ил Рилли л ти л ии х х илм л ил — ил+ т (кл.1 кл)1лг Я=Р+ !и~ гл Ел+тол~ (У+1) т Рл= -ля 1 -,+! гл.)- 1 л ! - !1 1; л= Бл+ й (Ил+Ил) ~ 1. (66а) гл (66б) (66в) и„х Рис.
98, (66г) Эти уравнения записаны в том порядке, который удобен для вычислений. Обсудим разностное выражение для вязкого давления (65). Чтобы выполнить предельный переход от разностной схемы к уравнениям газодинамики, надо сначала устремить т и и! к нулю при фиксированном коэффициеите вязкости, а затем построить серию таких предельных решений для неограниченно уменьшающихся значений ь. Но это очень трудоемко. Поэтому на практике объединяют этн предельные переходы в один об1ций, полагая йл=- = рлр(бг)' и ь,=р„рсбг, хотя законность такой процедуры пе ОДНОМЕРНЫЕ УРЛВНЕНИЯ ГЛЗОДННЛМИКИ 445 доказана (плотность введена в формулу для того, чтобы коэффициенты р» были безразмерны).
Таким образом, вязкое давление (65) принимает вид Шл = РЯР» (О„„— Ол)' — Р,спл (О,, — Пл) ПРИ Олы — О„< О, (67а) шл= О пРи Ели — О„) О. (67б) где с- Р'др!др — скорость звука. Выражение (67) написано для плоского случая; но обычно им пользуются при любой симметрии задачи. Аппроксимация. Из вида шаблона на рис. 98 и симметричного написания схемы (66) нетрудно заметить, что на течениях без сжатий, когда псевдовязкость (67) обращается в нуль, схема лкрест» имеет локальную аппроксимацию 0(т»+и»).
На течениях со сжатиями (в том числе — с ударными волнами) псевдовязкость отлична от нуля. Правда, квадратичный член в (67а) имеет величину 0 (й»); но линейный член имеет величину 0 (»1) и, тем самым, ухудшает порядок аппроксимации. Кроме того, вязкие члены записываются не вполне симметрично по времени. В итоге аппроксимация ухудшается до 0(т+Л). Нахождение разностного решения. Схема (66)— явная; вычисления по ней проводятся следующим образом.