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

Исходная геометрия расчётной областизадавалась в эти программы в виде полигона, ограниченного кусочнолинейнойграницей,криволинейныхвсвязиэлементов.Восотсутствиемвсехслучаяхвозможностивыполнениязаданияпрограммполучаются схожие минимальные углы триангуляции и параметры b. Врезультате работы программы gmsh образуется несколько большее числоэлементов,связанноесточнымвыполнениемкритерияневозможностью ограничения минимального размера элементов.44Делонеи3.

Моделирование течения вязкого газа с использованиемнеструктурированных сеток3.1.Численная схемаВ случае высокоскоростных течений при высокой разрешающейспособности сетки величина схемной вязкости оказывается сравнимой свеличиной физической вязкости. Схемная вязкость зависит от размера ячееки самого вида расчётной сетки, т.е.

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

Также к плюсам этих сеток относитсявозможность повышения разрешающей способности в интересующихобластях без введения подсеток и точное совпадение границ расчётнойобласти с заданной геометрией.Длярешенияуравненийвчастныхпроизводныхнанеструктурированных сетках применяется метод контрольных объёмов [7].ОсновнойидеейэтогометодаявляетсяиспользованияформулыОстроградского-Гаусса для построения разностной схемы:  dV   n  vdS ,V(3.1)Vгде n – внешняя нормаль к границе элемента.45В общем виде дифференциальное уравнение переноса величины φзаписывается в виде: ui       , txi xi xi (3.2)где:φ – физическая величина (скорость, плотность, температура, энергия);γ – коэффициент переноса (коэффициент вязкости, диффузии);σ – силы, не входящие в конвективную и диффузионную части.При интегрировании этого уравнения по конечному объёму получаетсяинтегральное уравнение вида: t  dV   x u  dV   xiVПрименяяiViVтеорему   dV   dV . xi VОстроградского-Гаусса(3.3)интегралпообъёмузаписывается через поверхностный интеграл: t  dV   x u  dV   xiViVi   dV   dV . xi V(3.4)Элементарные конечные объёмы, попарно непересекающиеся друг сдругом, строятся вокруг точек, в которых рассчитываются неизвестные.

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

К томуже при этом в центрах масс контрольных объёмов мы имеем осреднённыезначения по контрольному объёму, которые представляют истинное значениев центре масс со вторым порядком точности, в то время как применениеконтрольных объёмов построенных вокруг вершин даёт лишь приближение спервым порядком точности.3.1.1.

Аппроксимация уравнений Навье-СтоксаСистема уравнений движения вязкого сжимаемого газа, описанная в(1.1 – 1.7) может быть записана в консервативном виде [65]: F G H0,t x y z(3.5)где:u  vwe(3.6)- вектор консервативных переменных, F,G,H – вектора потоковuu  p   xxuv   yx2F(3.7),uw   zxu e  p    xx u   xy v   xz w  q x47vGvw   zyuv   xyv 2  p   yy,(3.8).(3.9)v(e  p)   yx u   yy v   yxz w  q ywuw   xzvw   yzHw   zzwe  p    zx u   zy v   zz w  q z2Где ρ, p, E – плотность, давление и полная энергия соответственно, u,v– декартовые (полярные) компоненты вектора скорости.Для замыкания данной системы уравнений можно использоватьуравнение состояния в виде:e  e , p  .(3.10)Для идеального газа уравнение состояния упрощается до вида:u 2  v 2  w2p    1 E  2E  e(3.11)γ – отношение удельных теплоёмкостей, в общем случае зависит оттемпературы и давления.Компоненты тензора вязких напряжений в предположении Стоксаопределяются из уравнений: xx 2  u v w  2   3  x y z  xy   yx      y x  uv  yy 2  v u w  2   3  y x z  xz   zx    zz 2  w v u  2 3  z y x  yz   zy     z y  w u   x z  v48w (3.12-3.18)Вектор тепловых потоков определяется из соотношений:q x  kTTT, q y  k, q z  kxyz(3.19-3.22)Данная система уравнений позволяет описывать как ламинарные так итурбулентные течения газа.

При ламинарном течении эффективная вязкостьполагается равной молекулярной и вычисляется из закона Сазерленда:T   *   T* 3/ 2T*  S 0T  S0Для(3.23)воздуха*  1.68  10 5значенияконстант:кг, T*  273.15K , S 0  100.5K .мсДля сверхзвукового течения в ряде случаев вязкостью газа можнопренебречь. Тогда система уравнений значительно упрощается, векторапотоков приобретают вид:uu  puv ,uwu e  p 2FvuvG  v 2  p ,vwv(e  pHwuwvw.(3.25-3.28)wwe  p 2При моделировании турбулентных течений эффективная вязкостьвычисляетсякаксуммамолекулярнойитурбулентнойвязкости,определяемой с использованием модели турбулентности.Плоскопараллельный случай получается из приведённой системы, еслипринять в ней w=0. В этом случае четвёртое уравнение системы становитсявырожденным, и система упрощается до вида:49t  f x    g y    0, u    , v   e  u  u  u 2  p  uv f   , g     2 uv  v  p e  p u e  p v (3.29)Система уравнений Эйлера в цилиндрической системе координат имеетвид:1  ru r  u z    r rz  1  ru z u r   u z2 p u   z  rrzz   0 .2t  1  ru r   u r u z   p ur  r rzr 1  ru r H  u z H  e    rrz (3.30)При умножении обоих частей каждого из уравнений на r, получаетсясистема,имеющаясхожийвидссистемойуравненийдляплоскопараллельного течения.

Данная система получается из системыуравнений Навье-Стокса () заменами y=r, uz=u, ur=v и введением вектораправой части S  0 0 p 0 T . Перепишем систему уравнений Навье-Стоксадля плоскопараллельного и осесимметричного течения в общем виде: t  f    g    S ,x u    , v   e yu2u  p   xx,f q  uv   yxe  p u   xx u   xy v 00 u v v S2 vp    y x y y 0  1     y.50uuv   xy,g q  v 2  p   yye  p v   yx u   yy v (3.31)где  – величина, отвечающая за выбор системы координат исоответственнотипатечения.При =0моделируетсядвумерноеплоскопараллельное течение,  =1 – трёхмерное осесимметричное.Рассмотрим область       R n , на которой построена сетка из Nэлементов: Ti    R n , int Ti   int T j   0, i  jNi 1Другими словами сетку, из элементов произвольной формы, которыепопарно не пересекаются, имеют общие границы и покрывают всюрассматриваемую область.Для каждого элемента контрольного объёма система уравненийзаписывается в интегральной форме: f g h dxdydz     dxdydz  0.t x y z V (3.32)После применения к первому слагаемому теоремы о среднем ипреобразования второго слагаемого по формуле Грина получим:ST  F ,.n dS  0, F  f , g , h .t (3.33)Система уравнений Навье-Стокса является инвариантной относительноповорота системы координат.

Таким образом, вместо вычисления отдельновекторов потоков f, g и h можно вычислить только вектор потока Fn  F ,.n  .Построим систему координат, связанную с границей расчётной областиизаданнуютремянеn  n x , n y , n z , l  l x , l y , l z , m  mx , m y , mz.ортогональнымиВкачествевекторавекторамиnпримемнаправление внешней нормали к границе расчётной области. Вектора l и m впространственном случае выбираются любые, перпендикулярные n. Вектор l51получается в результате векторного произведения вектора n и любогонеколлинеарного вектора t. Для уменьшения ошибки, связанной с машиннойточностью вычислений, вектор t выбирается по следующему алгоритму.Рассматривается тройка векторов i, j, k, образующих исходную системукоординат. В качестве вектора t используется вектор с минимальнойпроекцией на вектор n, где проекция определяется как скалярноепроизведение:    in  i , n ; in   j , n ; k n  k , n .(3.34)Тогда вектора l и m получаются из соотношений:  l  n , t ;  m  n, l (3.35)Вычислим нормальную составляющую вектора потока в системекоординат границы:Fn  n x F  n y G  T 1 F   , где100000nxnynz0T 0lxlylz00 mxmymz0001001(3.37)00000 nxlxmx0T 1  T T  0 n ylymy00 nzlzmz000010(3.36)(3.38)52u nu n2  pu n vF   u n wu n e  p (3.39)После преобразований получим:u nu n u  pn xFn  Fn x  Gn y  Hnz  T F TU   u n v  pn y ,u n w  pn zu n e  p 1гдекомпонентыскоростивповёрнутой(3.40)системекоординатопределяются из выражений:u n  un x  vlx  wm xv  un y  vl y  wm y(3.41)w  un z  vlz  wm zСледует отметить, что в двумерном случае вектора, образующиесистему координат, связанную с границей, однозначно вычисляются поформулам:n  nx , n yl   n y , nx(3.42)3.1.2.

Модель турбулентностиВ качестве модели турбулентности используется модель СпалартаАллмараса [18]. В этой модели предполагается, что мелкие вихри, размеркоторых не превышает размера ячеек, являются изотропными. Коэффициентвязкости представляет собой эффективный коэффициент, включающий всебя молекулярную вязкость газа и турбулентную вязкость: eff    t .53Турбулентная вязкозть  t определяется как t  v~f v1 .

Уравнение длявихревой вязкости:2~~ ~v~v~1  ~~ v~ ~  v   cb 2 v vuj cb1 1  f t 2 S v  cw1 f w   vvtx jxk   xk xk d   xk (3.43)1/ 6 1  cw6 3  6f v1  3,f1,fgv2w6 1  f v1  cv31 g  c w3 3v~v~6  , g  r  cw2 r  r , r  ~ 2 2vSk d~S STijv~, S  2Tij Tijk 2d 2 f v2- тензор завихренности, d – расстояние до стенки.Используемые в модели константы:1  cb 2  , c  0.3, c  2, k  0.41c2cb1  0.1355, cb 2  0.622, cv17.1 ,   , cw1  b21 w2w33kВ случае, если размер ячеек в пристеночной области больше толщиныпограничного слоя, модели турбулентности не позволяют верно вычислитьповедение течения вблизи стенки.

В этом случае для улучшения точностимоделированияиспользуетсяметодпристеночныхфункций.Впредположении о полностью развитом турбулентном характере течениясправедливы следующие утверждения:- профиль скорости перпендикулярный стенки имеет универсальныйхарактеринезависитотчислаРейнольдса;- градиентом давления вдоль потока можно пренебречь;- величина турбулентного сдвигового напряжения постоянна вдольслоя и равна касательному трению на стенке;54Для обезразмеривания уравнений газовой динамики в пристеночнойобласти вводится динамическая скорость:u w.(3.44)В этом случае безразмерная скоростьU Uu(3.45)и безразмерное расстояние до стенкиy yu(3.46)Пристеночная область разбивается на три зоны:U   y  , y   0,5y   5,30  1U  k ln Ey , y  30,150(3.47)1) вязкий подслой y   0,5 с линейным профилем скорости U   y 2) переходный подслой y   5,303) турбулентный логарифмический слой y   30,150 с логарифмическойзависимостью скорости от расстояния до стенки U   ln Ey   , где1kk=0.41 – константа Кармана, E =8.4.Из определения динамической скорости и безразмерной скоростикасательное напряжение на стенке:w  y U.U y(3.48)55При использовании пристеночных функций первый узел сетки долженнаходится в логарифмическом подслое, т.е.