Диссертация (1149881), страница 8

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 8)

Сетка построена с учётом задания периодических47граничных условий: если совместить верхнюю и нижнюю границы блока,то границы сетки будут совпадать "узел в узел".Рис. 2.1: Вид расчётной области.Рис. 2.2: Вид расчётной сетки, слева − около роторной лопатки, справа −увеличенная область около передней части лопатки.482.3.2Численный алгоритмДля численного интегрирования исходных уравнений используется метод конечного объёма. Рассмотрим численный алгоритм на примере уравнений Навье−Стокса. Для уравнений Рейнольдса процедура расчёта аналогична. Интегрируя уравнения (2.1) по объёму Ωi ячейки Ci и применяя теорему Остроградского−Гаусса, получим интегральную форму записи уравнений Навье−Стокса для данного объёма (рис.

2.3):ZIII∂Q dΩ + FdS = GdS + HdΩ,∂tΩiSiSi(2.30)ΩiF = Fx nx + Fy ny , G = Gx nx + Gy ny .Здесь Si − поверхность, ограничивающая объём Ωi . В случае плоскоготечения в качестве объёма Ωi рассматривается площадь i-й ячейки сетки.При этом величина поверхности Si будет равна длине контура, ограничивающего данную ячейку.Рис. 2.3: Элемент расчётной сетки.При вычислении интеграла в (2.30) по объёму Ωi используем теорему о среднем значении. Рассмотрим в качестве примера четырёхугольную ячейку.

Контур Si представляется в виде суммы четырёх отрезков Sij49(j = 1, 2, 3, 4), соответствующих граням ячейки. Интегралы по Si вычисляются как сумма интегралов по этим граням. Обход контура осуществляется против часовой стрелки. Обозначим верхними индексами n и n + 1значения параметров газа соответственно на n-м и (n + 1)-м шаге по времени. Значения векторов Fj и Gj вычисляются в центре соответствующихграней.

Тогда аппроксимацию 1-го порядка по времени для интегральногосоотношения (2.30) для i-й ячейки можно записать в виде:4Qn+1i=Qni4∆t X∆t XFij Sij +Gij Sij + ∆tHi .−Ωi j=1Ωi j=1(2.31)В данной работе векторы Fij ("невязкие" потоки консервативных переменных через j-ю грань i-й ячейки) вычисляются по схеме Роу [109] сиспользованием процедуры энтропийной коррекции, предложенной Хартеном (Yee H.C., Harten А., 1987) [142].Опишем кратко схему расчёта в применении к произвольно ориентированной грани в плоском случае. Введём в рассмотрение матрицу ЯкобиA = ∂F/∂Q. Её собственные значения λ(1) , λ(2) , λ(3) , λ(4) получаются изрешения характеристического уравнения det(A − λI) = 0 и в рассматриваемом случае соответственно равны vn , vn , vn +c и vn −c, где vn = vx nx +vy ny− проекция вектора скорости газа на внешнюю нормаль к стороне ячейки(nx и ny − проекции единичной нормали на оси x и y).

Обозначим какR и L матрицы, столбцами и строками которых соответственно являютсяправые и левые собственные векторы матрицы A. Эти матрицы имеют вид:1ρ2c0ρ2cρρvxρny(vx + cnx )(vx − cnx )2c2cR=ρρvy−ρnx(vy + cny )(vy − cny )2c2c 1ρρ(vx2 + vy2 ) ρvτ(H + cvn )(H − cvn )22c2c50,L=γ−1 2vxvyγ−1M(γ − 1) 2(γ − 1) 2− 22cccny−nxvτ0−ρρρ.c γ − 1 2 vn1vx1vyγ−1 M −nx − (γ − 1)ny − (γ − 1)ρ 2c ρcρcρcc γ − 1 2 vn1vx 1vy γ − 1 M +− nx + (γ − 1)− ny + (γ − 1)ρ2cρcρcρc1−Здесь H = cp T + (1/2)(vx2 + vy2 ), vτ = vx ny − vy nx , M − число Маха,γ = cp /cv − показатель адиабаты.Рассмотрим вычисление потока для j-й грани Fij в правой части (2.31).Опустим для краткости записи индекс ячейки i.

При обходе ячейки вдольеё границы против часовой стрелки, содержимое ячейки всегда находитсяс левой стороны от её грани, а параметры в любой соседней ячейке c правой стороны. Обозначим индексами l и r параметры газа на j-й грани врассматриваемой ячейке (слева от j-й грани) и на той же грани в соседней ячейке (справа от j-й грани). Тогда вектор Fj в (2.31) вычисляетсяследующим образом:Fj = 0.5(F(Qr ) + F(Ql ) − RΛL(Qr − Ql )).Λ=(1)ψ(λj )0000(2)ψ(λj )0000ψ(λj )0,000ψ(λj )(3)(4)Элементы матриц R, L и собственные числа λ(k) (k = 1, 2, 3, 4) относятся к грани между ячейками и вычисляются по осреднённым параметрамгаза между ячейками. Осреднение в методе Роу происходит по следующимформулам (Гильманов А.Н., 2000) [12]:√√√√vxl ρl + vxr ρrvyl ρl + vyr ρrvx =,vy =,√√√√ρl + ρrρl + ρr√√qHl ρl + Hr ρrH= √,c = (H − (vx2 + vy2 )/2)(γ − 1),√ρl + ρr√ √ρ = ρl ρr ,51функция ψ(λ) задаётся соотношениями (Yee H.C., Harten А., 1987) [142]:(|λ|,|λ| ≥ ,ψ(λ) =(λ2 + 2 )/2,|λ| < .Здесь − малый положительный параметр, который в расчётах полагался равным 0.1.

Последнее соотношение представляет собой энтропийнуюкоррекцию по Хартену. Таким образом происходит коррекция собственногозначения λ, когда оно принимает значение близкое или равное нулю, чтоможет привести к возникновению осцилляций в численном решенииДля построения численной схемы первого порядка аппроксимации попространству векторы консервативных переменных Ql и Qr вычисляютсяна основе средних значений параметров в ячейке.

Для построения схемывторого порядка используется линейная реконструкция параметров газав ячейках φ = {ρ, vx , vy , p}. Распределение каждого из этих параметроввнутри ячейки представляется с помощью билинейной функции:φ(x, y) = φi + ai (x − xi ) + bi (y − yi ),(2.32)где φi − среднее значение газодинамического параметра в ячейке, ai и bi− составляющие вектора градиента данной функции. Рассмотрим нарядус i-й ячейкой любые другие две соседние ячейки, которые являются такжесоседними между собой. Зная параметры газа в центрах этих трёх ячеек,найдём коэффициенты ai и bi , и таким образом построим для выбранногопараметра билинейное восполнение в пределах контрольного объёма.

Повторим эту процедуру для всех возможных наборов из троек соседних ячеек (для четырёхугольных ячеек таких наборов будет четыре) и вычислим(m)(m)набор коэффициентов ai , bi(m = 1, 2, 3, 4) (рис. 2.4). Далее на основечетырёх полученных наборов составляющих вектора градиента рассматриваемого параметра вычислим окончательные значения ai и bi с использованием ограничителя "minmod" (Yee H.C., Harten А., 1987) [142]:(m)ai = minmod(ai ),(m)bi = minmod(bi ),52m = 1, 2, 3, 4,(m)minmod(fi)=(1)(4)(1)(1)(4) min(|fi |, ..., |fi |), sign fi , если signfi = ... = signfi0,в противном случае. 1, x > 0sign(x) =0, x = 0 −1, x < 0Используя найденные значения ai и bi вычислим параметры газа в центре j-й грани i-й ячейки по формуле (2.32) и по ним вычислим все компоненты вектора консервативных переменных Ql .

Аналогично вычисляетсявектор Qr на той же грани со стороны соседней ячейки. В областях гладкого течения данная процедура всегда даёт наименьший градиент из четырёхвозможных вычисленных вариантов, в точках же экстремума газодинамического параметра она даёт нулевое значение градиента и тем самым обеспечивает подавление возможных осцилляций решения.Рис. 2.4: Шаблон для процедуры реконструкции газодинамических параметров.Источниковый член H в (2.31), описывающий обратное влияние твёрдой53примеси на несущий газ определяется следующим образом [137]:0NpXfHpxpppiH = −hH i, Hi = , f  , Hi = −Ωi py i=1apгде N − количество частиц примеси в рассматриваемой ячейке, Hpi − вектор параметров для каждой отдельной i-й частицы, fpx и fpy − компонентысилы, действующей на частицу со стороны газа, ap = fp · vp + lp · ! p + qp− скорость изменения полной энергии i-й частицы вследствие взаимодействия с газом, fp − сила аэродинамического сопротивления частицы, lp −аэродинамический момент, действующий со стороны газа на частицу, vp ,!p,qp − поступательная и угловая скорости частицы, полный тепловойпоток от газа к частице.Число Нуссельта рассчитывается по формуле Кавано−Дрейка [44], которая учитывает влияние сжимаемости газа на теплообмен:qp = 2Nup πrp λ(T − Tp ),Nu0p,Nup =1 + 3.42Nu0p Mp /(Rep Pr)0.33.Nu0p = 2 + 0.459Re0.55p PrЗдесь Nu0p − число Нуссельта для течения несжимаемой сплошной среды,√Mp = |v − vp |/ γRT − число Маха частицы, Rep = 2rp ρ|v − vp | − числоРейнольдса при обтекании частицы.Рассмотрим процедуру для вычисления "вязких" потоков Gij в (2.31).Векторы Gx и Gy содержат компоненты τxx , τyy , τxy , qx , qy .

Для их вычисления мы должны знать значения следующих параметров на грани:vx , vy ,∂vx ∂vy ∂vx ∂vy ∂T ∂T,,,,,.∂x ∂x ∂y ∂y ∂x ∂yЗначения переменных vx , vy , T на грани можно вычислить с помощьюосреднения (рис. 2.3):fij = (1/2)(fci + fcj ),54где f − любой газодинамический параметр, fci и fcj − значение параметрана грани со стороны i−й или j−й ячейки.Для вычисления производных ∂f /∂x и ∂f /∂y вводится локальная система координат (ξ, η). Вычисление производятся по следующим формулам:fcj − fci ∂ffn2 − fn1∂f=,=,∂ξ∂ξ∂η∂ηгде fn − значение газодинамического параметра в узле. Оно вычисляетсяс помощью процедуры интерполяции и осреднения по всем примыкающимк узлу ячейкам.Производные ∂f /∂x и ∂f /∂y вычисляются с помощью соотношений:∂f ∂ξ ∂f ∂η ∂f∂f ∂ξ ∂f ∂η∂f=+,=+,∂x∂ξ ∂x ∂η ∂x ∂y∂ξ ∂y ∂η ∂y∂y −1 ∂ξ∂x∂η∂y∂η∂x −1∂ξ=J ,= − J −1 ,= − J −1 ,=J ,∂x ∂η∂y∂η∂x∂ξ∂y∂ξJ=2.3.3∂x ∂y ∂x ∂y−.∂ξ ∂η ∂η ∂ξГраничные и начальные условияГраничные условия для уравнений Эйлера и Навье−СтоксаВвиду одинакового шага решёток и предполагая отсутствие глобальныхотрывов, распространяющихся на несколько лопаток решёток, на верхнейи нижней границах каждого блока ставятся условия периодичности решения по пространству.

На границе между роторным и статорным блокамииспользуется процедура согласования решений (см. п.2.2.4). При постановке граничных условий на входной границе роторного блока имеется особенность, связанная с тем, что рассматриваемое течение дозвуковое и этаграница близка к решётке. В этом случае возмущения распространяются израсчётной области вверх по потоку и на входной границе нельзя задаватьпараметры невозмущённого потока.

На входной границе характеристика55Римана, соответствующая скорости распространения возмущений vn + c, иэнтропийная характеристика (соответствует скорости потока vn ) "входят"в расчётную область, а характеристика Римана, соответствующая скоростиvn − c, "выходит" из неё. Поэтому на этой границе следует задавать два параметра, а один экстраполировать из расчётной области. В качестве задаваемых параметров обычно выбираются такие, которые подвергаются наименьшему влиянию со стороны обтекаемого объекта. В данной работе наγ p∞ V∞2+входной границе задаются значения полной энтальпии h0 =γ − 1 ρ∞2γи энтропийной функции s = p/ρ в невозмущённом течении.

Третий параметр − давление − экстраполируется из расчётной области. Тогда, считаяпоток перпендикулярным к входной границе, остальные параметры определяются следующим образом:ρ=pss 1γ , vx = 2 h0 −γ pp, vy = 0, T =.γ − 1ρρRАналогичные рассуждения для выходной границы приводят к тому, чтопри дозвуковой скорости на этой границе следует задавать один "внешний"параметр газа, а два экстраполировать из расчётной области, что и получалось в расчётах. В данном исследовании, исходя из опытных данных,на выходной границе задаётся давление в 1.2 раза большее, чем статическое давление в невозмущённом натекающем потоке, остальные параметрыэкстраполируются из расчётной области: ρ, vx , vy .На контурах профилей для уравнений Эйлера ставится условие непротекания vn = 0.

Для уравнений Навье−Стокса ставится условие непротекания vn = 0 и прилипания vτ = 0, а также задаётся температура стенкиTw = 288.15 К.Граничные условия для модели турбулентности k − ω SST.Для граней на поверхности лопатки задаются условия для функций kи ω:800µt,(∆y)2где ∆y − расстояние от рассматриваемой поверхности стенки до центраk = 0, ω =56прилегающей ячейки (Menter F.R., 1994) [95]. Такое граничное условиесправедливо если относительный размер прилегающей ячейки удовлетворяет условию y + ≤ 3.Для сеток с более крупными ячейками вблизи стенки используютсяспециальные пристеночные функции [43], [96]:spµVwlog1/4subsub 4log 4 0.25kw ,, uτ = Cµuτ = ((uτ ) + (uτ ) ) , uτ =ρ∆yuτ ρ∆y1y+ =, u+ = ln(y+ ) + B,µκrutauµsub 4log 4 0.25subCτ = ((Cτ ) + (Cτ ) ) , Cτ = uτ, Cτlog = + ,ρUw ∆yu 26ρ uτ1 ρ u2τsub 4log 4 0.25sublogωw ((ωw ) + (ωw ) ) , ωw =,, ωw = p++βµyµyκCµ kwτwPk =, τw = ρCτ Vw ,ωwκµy +где Vw , kw , ωw - скорость газа и значение величин k и ω в пристеночнойячейке, B = 5.5, Cµ = 0.09, β = 0.075, κ = 0.41, V∞ − скорость натекающего потока, Pk − член ответственный за генерацию турбулентности впристеночной ячейке.На входной границе расчётной области задаётся интенсивность турбулентности I = 5% и турбулентная вязкость µt = (1...10)µ.

Характеристики

Список файлов диссертации