Бураго Н.Г. Вычислительная механика (1185926), страница 30
Текст из файла (страница 30)
То есть, варианты дифференциальнойформулировки с условием несжимаемости и с уравнением Пуассонадля давления эквивалентны.Дискретизированные уравнения движения для скоростей иуравнение Пуассона для давления обеспечивают выполнениедискретизированного условия несжимаемости только приближенно,причем в процессе пошагового интегрирования по временипогрешность в условии несжимаемости растет.17.4. Метод коррекции давленияВ работах (Chorin,1968), (Temam, 1969), (Fortin, 1971),(Гущин, Щенников, 1974) были предложены варианты метода176Глава 17.
Расчет несжимаемых теченийкоррекции давления (pressure correction method, projection method),который позволяет обеспечить выполнение условия несжимаемостина дискретном уровне. Метод коррекции давления можно трактоватькак расщепление по физическим процессам (Гущин, Щенников,1974). Для реализации метода на каждом шаге по времени расчетпроводится в три этапа.1-й этап. Определение новой скорости uɶ n +1 без учетадавления:uɶ n +1 − u nµVρn nnnn+ u ⋅∇u = ∇ ⋅ ( ∇u ) + g∆t nρ0ρ0где (важно!) на скорость uɶ n +1 граничные условия не накладываются.“Скорость” uɶ n +1 это просто обозначение суммы членов уравнениядвижения за исключением члена с градиентом давления.2-й этап. Определение давления p n +1 из решения краевойзадачи для уравнения Пуассона:∇ ⋅ uɶ n +11= ∇ ⋅ ( ∇p n +1 )∆t nρ0с учетом главных граничных условий для давления (если давлениезадано на какой-то части границы Sp ⊂ S )x ∈ Sp :p n +1 = p* (x, t n +1 )и естественных граничных условий для давления получаемыхпроектированием уравнения движения на нормаль к границе:x ∈ S \ Sp : (u n +1 − uɶ n +1 ) n ⋅∇p n +1 = −n ⋅ ρ0∆t nВ граничных условиях для давления надо учитывать, что:а) на входной границе должна быть задана или проекция граничнойскорости на нормаль к границе, или давлениеx ∈ S \ Sp = Sun : u n +1 ⋅ n = u* (x, t n +1 ) ⋅ n или p n +1 = p* (x, t n +1 )б) на стенках задана нормальная скорость и либо касательныескорости, либо соответствующие касательные проекции граничных177Глава 17.
Расчет несжимаемых теченийсил трения.в) на открытых границах условия для градиента давленияполучаются проектированием векторного уравнения движения нанормаль к границе с учетом условий продолжения решения(“мягких” граничных условий для .искомых функций).г) на свободных границах давление уравнивается с давлением вовнешней среде.3-й этап. Определяется новая скорость u n +1u n +1 − uɶ n +11= − ∇p n +1∆t nρ0с учетом главных (кинематические) граничные условия дляскоростейt > 0 , x ∈ Su : u n +1 = u* (x, t n +1 )Уравнения 1-го и 3-го этапов в сумме аппроксимируютуравнения движения.
Уравнение 2-го этапа является результатомскалярного умножения уравнения 3-го этапа на операторпространственного дифференцирования с учетом условиянесжимаемости ∇ ⋅ u n +1 = 0 .Очевидно, имеется бесконечное множество способовдискретизации и последующего решения уравнений описанногометода коррекции давления. Oбсуждение нюансов и особенностейконкретных реализаций, присущих различным вариантамизложенного метода, здесь опущены.Заметим, что можно включить в уравнение первого этапачлен с градиентом давления ∇p n с n-го временного слоя, тогда навтором и третьем этапах вместо нового давления на (n+1)-мвременном слое будет фигурировать поправка к давлению∆p n +1 = p n +1 − p n .
При этом граничные условия 2-го этапа такжеформулируются для приращений давления.17.5. Переменные “функция тока – вихрь”В переменных “функция тока – вихрь” скоростьопределяется через векторную функцию тока Ψuu = ∇×Ψблагодаря чему условие несжимаемости ∇ ⋅ u = 0 удовлетворяется178Глава 17. Расчет несжимаемых теченийтождественно. Определение вихряω = ∇×uи выписанные выше уравнения движения в скоростях черезфункцию тока и вихрь переписываются в виде уравнений дляфункции тока и вихря соответственно∇ × (∇ × Ψ ) = ω∂ω+ u ⋅ ∇ω − ω ⋅∇u = ∇ ⋅ (ν V ∇ω) + ∇ × (ρg / ρ0 )∂tгде ν V = µ V / ρ0 - кинематическая вязкость.
Эти два уравненияиспользуются в формулировке “функция тока – завизренность”вместо уравнения неразрывности и уравнения движения.Заметим, что давление в данной формулировке не входит вразрешающую систему уравнений для функции тока и вихря и оноопределяется отдельно после определения поля скоростей изуравнения Пуассона для давления∇ 2 p = −ρ0∇ ⋅ (u ⋅∇u) + ∇ ⋅ (∇ ⋅ (µ V ∇u)) + ∇ ⋅ (ρg )которое выводится с учетом условия несжимаемости скалярнымумножением уравнения движения на оператор пространственногодифференцирования.
При этом граничные условия для давления натой части границы, гле оно неизвестно, определяютсяпроектированием (скалярным умножением) векторного уравнениядвижения на нормаль к границеn ⋅∇p = n ⋅ [−ρ(∂ t u + u ⋅∇u) + ∇ ⋅ (σV ) + ρg ]Граничные условия для функции тока и вихря имеют видограничений на граничные значения этих функций или ихнормальных к границе производных. Эти условия должны бытьсогласованы с распределением граничных скоростей, то есть награнице должно выполняться соотношение∇ × Ψ = u*которое позволяет переформулировать ограничения, накладываемыена граничные скорости, в виде ограничений на функцию тока.Граничные условия для вихря получаются аппроксимациейвыражения ω = ∇ × u с учетом граничных условий для скорости и179Глава 17. Расчет несжимаемых теченийее производных.
Важно заметить, что во избежание счетныхколебаний решения во времени целесообразно использоватьрелаксацию (ослабленную форму) граничных условий для вихряx ∈ Sω :ω n +1 = ω n + α1 (∇ × u n +1 − ω n )где α1 ≈ 0.25 - коэффициент релаксации. Граничные условия длявихрей задаются на участках границы, служащих источникамивихрей, а именно на стенках и входных границах, а также награницах симметрии, где вихри, как правило, зануляются.
Наостальных границах ставятся “мягкие” условия продолжениярешения, а именно, зануляются нормальные к границе производныеот вихря.Формулировка граничных условий для функции тока и вихрязависит от используемых формул дискретного представленияпроизводных. Неоднозначность разностных представлений означает,что не существует единого подхода к записи граничных условий иони всегда являются приближенными (подробное описаниераспространенных вариантов дано в книге Роуча (1980)).Плотность в несжимаемом течении.В неоднородноймногофазной несжимаемой среде плотность может быть переменнойв точке пространства в зависимости от того, какая фаза занимаетданную точку пространства в данный момент времени. Подстановкаусловия несжимаемости в уравнение закона сохранения массы даетдля плотности несжимаемой многофазной среды транспортноеуравнение∂ t ρ + u ⋅∇ρ = 0Транспортное уравнение показывает, что значение плотностинесжимаемой среды постоянно вдоль лагранжевых траекторий илитраекторий материальных частиц ( dρ/dt=0 ).
Для однофазной средыпостоянные значения плотности принимаются для всего объемазанятого несжимаемой средой и уравнение переноса плотностиудовлетворяется тождественно.Однако, при рассмотрении многофазных несмешивающихсясред (в частности, течений со свободными границами) плотностьнесжимаемой среды может иметь разные значения для различныхфаз и, поэтому, может претерпевать скачок на межфазных границах.В этом случае для расчета движения фаз и идентификации ихвзаимного положения необходимо интегрировать транспортноеуравнение и отслеживать межфазные подвижные границы.
В этомслучае плотность, определяемая решением транспортногоуравнения, может в численном решении принимать значения,180Глава 17. Расчет несжимаемых теченийотличные от номинальных значений для данной фазы. Несмотря наконсервативность дискретных аппроксимаций и выполнениеусловия несжимаемости полем скоростей при этом баланс фазовыхмасс может нарушаться из-за погрешностей в определенииположения и диффузного “размазывания” межфазных границ. Этотребует контроля закона сохранения массы и введения в алгоритмырешения корректировок, обеспечивающих консервативность.Безвихревые течения.
Задачи для рассматриваемой системыуравнений сильно упрощаются в случае идеальной (невязкой) средыпри отсутствии источников вихрей. В этом случае решается толькоуравнение для функции тока при заданных на границах областирешения граничных условиях Дирихле или Неймана.17.6. Переменные “функция тока –завихренность”В случае двумерных течений отличными от нуля являютсяодин компонент вектора вихряωθ = ω =∂u ∂w−∂z ∂rназываемый завихренностью, и один компонент векторной функциитока ϕθ = ϕ , с помощью которого определяются компонентыскорости.
В этом разделе для обозначения пространственныхпеременных двумерной задачи использованы обозначения r и z ,которые в случае осесимметричных течений (параметр геометрииα = 1 ) означают радиальную и осевую координаты, а в случаеплоских течений ( α = 0 ) отвечают координатам x и y . Выражениядля скоростей имеют вид:1 ∂ϕr α ∂z1 ∂ϕuz = w = − αr ∂rur = u =благодаря чему уравнение неразрывности удовлетворяетсятождественно.Будучиподставленнымивопределениезавихренности, выражения для скоростей дают двумерное уравнениедля функции тока181Глава 17. Расчет несжимаемых течений∂ 1 ∂ϕ ∂ 1 ∂ϕ + =ω∂r r α ∂r ∂z r α ∂z Двумерное уравнение для завихренности имеет вид∂ωωu∂ v2 ∂T− α α − α α − gβ=∂tr∂z r ∂r=1 ∂ α ∂ω ∂ ∂ω ω νr + ν − αν 2αr ∂r ∂r ∂z ∂z rДругой способ определения функции тока имеет вид:ψ = ϕ / rαтогда скорости определяются соотношениямиur = u =∂ψψ∂ψ−αи uz = w = −∂rr∂zи уравнение для функции тока выглядит так1 ∂ α ∂ψ ∂ ∂ψ ψr+ −α 2 = ωαr ∂r ∂r ∂z ∂z rВопрос о том, какой из описанных способов введения функции токаприменять на практике, однозначного ответа не имеет.
В литературевстречаются примеры применения обоих способов.17.7. Методы в переменных “функция тока– вихрь”Для реализации начально-краевых задач в переменных“функция тока – вихрь”, как правило, применяются вариантыпроекционных методов (сеточных или бессеточных). Конвекциявихрей рассчитывается чаще всего по какой-либо явной схеме, а дляфункции тока и для диффузионных членов в уравнении для вихряприменяются неявные аппроксимации. На каждом шаге по временив нестационарных задачах или на каждой итерации в стационарных182Глава 17.