Главная страница
Навигация по странице:

  • 13.1. Общие принципы численного моделирования

  • 13.2. Математическая постановка задачи

  • 13.3. Конечно-разностная аппроксимация и метод решения

  • Задача о распаде разрыва.

  • 13.4. Аппроксимация и устойчивость численного решения

  • 13.5. Течения с разрывами параметров

  • 13.6. Нефизическое поведение решения

  • 13.7. Контрольные вопросы

  • Гидравлический пресс. Механика жидкости и газа


    Скачать 5.98 Mb.
    НазваниеМеханика жидкости и газа
    АнкорГидравлический пресс
    Дата01.06.2022
    Размер5.98 Mb.
    Формат файлаpdf
    Имя файлаTPM_Zezin_mehanika_zhidkosti_gaza.pdf
    ТипУчебное пособие
    #561398
    страница17 из 22
    1   ...   14   15   16   17   18   19   20   21   22
    ГЛАВА 13. ПРИМИНЕНИЕ ЧИСЛЕННЫХ МЕТОДОВ В
    МЕХАНИКЕ ЖИДКОСТИ И ГАЗА
    Многие задачи механики жидкости и газа, с которыми приходится иметь дело инженерам и исследователям, не поддаются аналитическому решению в силу нелинейности описывающих эти задачи систем уравне- ний
    *
    . Именно поэтому при решении рассмотренных выше и важных для практики задач использовались упрощения при их математической поста- новке, например, допущение об одномерности потока, несжимаемости жидкости и пр., а для описания отдельных факторов, влияющих на реше- ние, применялись экспериментальные данные, например, по величине ко- эффициента гидравлического трения и др. И даже при таких упрощениях, например, такая, казалось-бы тривиальная задача определения расхода в простом трубопроводе при заданном перепаде давлений не имеет аналити- ческого решения. Поэтому единственной возможностью для теоретическо- го анализа сколько-нибудь сложных задач механики жидкости и газа явля- ется получение их численного решения.
    В данном разделе кратко поясняются основные идеи численного моде- лирования, дается общее представление о разностных схемах, отражаются некоторые особенности, которые необходимо учитывать при выполнении расчетов и интерпретации полученных результатов.
    Для более подробного изучения вопросов численного моделирования отправляем интересующихся к обширной специализированной литературе, в качестве которой может служить, например [6,…,12].
    13.1. Общие принципы численного моделирования
    Численное решение задач МЖГ, описываемых дифференциальными уравнениями, основано на их сведении к задачам решения систем алгеб- раических уравнений. Для этого расчетная область разбивается на ячейки
    (формируется расчетная сетка). Непрерывное распределение гидродина- мических параметров
    **
    заменяется дискретным со значениями, привязан- ными к узлам сетки, либо к центру ячеек. Распределение гидродинамиче- ских параметров по времени также заменяется дискретным со значениями, привязанными к каждому временному слою. Пространственная сетка при переходе с одного временного слоя на другой может оставаться как неиз- менной, так и трансформироваться, адаптируясь под текущее поле гидро-
    *
    Более того, доказательство существования и гладкости решения уравнений Навье-
    Стокса относится к одной из шести нерешенных до сих пор математических проблем тысячелетия, сформулированных в 2000 г институтом Клэя (Кембридж).
    **
    Здесь имеется в виду, что в расчетной области отсутствуют скачки уплотнения.
    При наличии скачков возможны подходы, когда они не учитываются при построении пространственной сетки и тогда скачки «размазываются», либо положение скачков специально вычисляется и ячейки сетки «привязываются» к поверхности скачка.

    206 динамических параметров. Таким образом производится пространственно- временная дискретизация поля искомых переменных.
    Производные, входящие в математическое описание исследуемого про- цесса заменяются разностными аналогами. В результате система диффе- ренциальных уравнений заменяется системой алгебраических уравнений, искомыми переменными в которой являются сеточные функции – дис- кретные значения поля гидродинамических параметров исходной задачи
    МЖГ. Описанная процедура носит название конечно-разностной аппрок-
    симации. Полученная система алгебраических уравнений решается каким- либо численным методом.
    Имеется несколько методов получения конечно-разностных моделей: разложение функций в ряд Тейлора, интерполяция функций полиномами, интегральный метод, метод контрольного объема. Конечно-разностная мо- дель и метод ее решения представляют собой разностную схему решения исходной системы дифференциальных уравнений.
    В настоящее время разработано множество разностных схем, отличаю- щихся точностью, с которой они аппроксимируют систему дифференци- альных уравнений, и методами решения получившейся системы алгебраи- ческих уравнений. Если искомые значения сеточных функций на следую- щем временном шаге выражаются явно через значения этих функций на текущий и предшествующий моменты времени, то такие разностные схе- мы называют явными. В противном случае разностная схема называется
    неявной. Неявные схемы менее требовательны к ограничениям на шаг ин- тегрирования, но требуют выполнения больше операций для вычисления сеточных функций на новом временном слое.
    Выбор той или иной разностной схемы зависти от типа решаемой зада- чи и располагаемых ресурсов ЭВМ.
    Рассмотрим основные этапы численного моделирования и некоторые особенности, которые могут возникнуть при этой работе на примере реше- ния простой задачи распространения ударной волны в сжимаемой жидко- сти в горизонтальном жестком трубопроводе постоянного сечения.
    13.2. Математическая постановка задачи
    На данном этапе на основе фундаментальных уравнений МЖГ и других физических законов разрабатывается система уравнений, описывающих рассматриваемый физический процесс. При этом обосновываются и при- нимаются допущения, упрощающие постановку и решение задачи (при не- обходимости).
    В рассматриваемом нами примере течение для простоты будем предпо- лагать одномерным с равномерным распределением гидродинамических параметров по поперечному сечению. Система уравнений для данной за- дачи, получается путем применения указанных выше допущений к общим уравнениям динамики (4.9), (4.21) и имеет вид:

    207

    уравнение неразрывности
    1 d
    0
    d
    w
    t
    x
     




    ;
    (13.1)

    уравнение Навье-Стокса
    2 1
    2
    w
    w
    p
    a
    t
    x
    x








     




     


    ,
    (13.2) где w,

    , p – средние по сечению трубопровода скорость, плотность и дав- ление соответственно.
    В уравнении (13.2) через а

    обозначена отнесенная к единице массы жидкости сила вязкого сопротивления, то есть осредненная по площади сечения канала величина
    2
    (
    )
    x
     
    u
    . Выразим ее через напряжение трения на стенках канала

    w следующим образом d
    d
    w
    w
    x
    a
    S x
    S

     
     




    , где

    , S – периметр и площадь сечения трубопровода соответственно. С учетом (8.18) последняя формула для круглого трубопровода диаметром d
    примет вид
    2 2
    w w
    a
    d w



    (13.3)
    Множитель w/|w| добавлен в последней формуле для того, чтобы сила вязкости действовала противоположно направлению скорости w. Для оп- ределения коэффициента гидравлического трения будем использовать формулы раздела 8.5.
    Преобразуем уравнения (13.1), (13.2) к удобному для решения виду.
    Для учета сжимаемости среды представим первое слагаемое в (13.1) в сле- дующем виде
    1 d
    1 d d
    1 d
    1
    d d
    d d
    p
    p
    p
    p
    w
    t
    p
    t
    E
    t
    E
    t
    x


















    ,
    (13.4) где Е = d
    / d
    p


    – модуль упругости жидкости. Тогда, учитывая, что
    2
    /
    E
    a
     
    (а – скорость распространения волн давления в жидкости) урав- нение (13.1) запишем следующим образом
    2 0
    p
    p
    w
    w
    a
    t
    x
    x




     




    (13.5)
    Упростим уравнение (13.2), учитывая, что сжимаемость реальных жид- костей, характеризуемая модулем упругости Е, как правило, незначитель- на. По крайней мере, в большинстве случаев справедливо соотношение
    E
    p
    
    . Тогда имеем

    208








    2 2
    /
    /
    /
    /
    1
    p
    p
    p
    p
    p
    p
    p
    p
    p
    p
    x
    x
    x
    x
    p x
    x
    E x
    x









    
     








     



     

     



    С учетом последнего соотношения и (13.3) уравнение (13.2) принимает вид
    2 2
    0.
    2 2
    w
    w
    p
    w w
    t
    x
    d w
















    (13.6)
    Итак имеем систему двух уравнений (13.5), (13.6) относительно двух неизвестных: давления р и скорости w. В качестве начальных условий за- даем


    0,
    0
    w t
    x


    ,


    0 0,
    p t
    x
    p


    , то есть условия покоящейся жидкости.
    Для формирования граничных условий при х = 0 будем считать, что трубо- провод стыкуется с большим резервуаром, давление в котором постоянно и равно р
    0
    (сам метод задания левых граничных условий описан ниже). На правом конце трубопровода х = L зададим переменные по времени гранич- ные условия: при
    0
    t
    t

    ,


    н
    ,
    p t x
    L
    p


    (где р
    н
    – давление окружающей среды); при
    0
    t
    t

    ,


    ,
    0
    w t x
    L


    . То есть считаем, что при времени
    0
    t
    t

    трубопровод в выходном сечении мгновенно перекрывается.
    13.3. Конечно-разностная аппроксимация и метод решения
    Для конечно-разностной аппроксимации рассматриваемой нами задачи используем интегральный метод. Для этого построим сетку с равномерны- ми шагами по координате
    1
    j
    j
    x
    x
    h



    и времени

    , рис. 13.1.
    Рис. 13.1. Схема расчетной сетки
    Параметры на границах ячеек обозначим прописными буквами. Они имеют целочисленный индекс. Их считаем постоянными на временном ин- тервале интегрирования

    . Параметры в середине ячейки (обозначенные строчными буквами и имеющие дробный индекс) считаем не зависящими от координаты в пределах шага сетки h. Эти параметры, относящиеся к нижнему временному слою, обозначаются с нижним расположением ин- декса, а относящиеся к верхнему временному слою – с верхним располо- жением, см. рис. 13.1. Истинные значения скорости и давления в трубо- проводе будут аппроксимироваться значениями параметров в серединах ячеек (с дробными индексами).
    t
    1 0
    N
    x
    j
    j+1
    j

    1
    t+

    t
    1 2
    1/2 1/2
    ,
    j
    j
    w
    p


    j+1/2 1/ 2 1/ 2
    ,
    j
    j
    w
    p


    ,
    j
    j
    W P

    209
    Проинтегрируем по объему ячейки составленные уравнения (по х от 0 до h, по t от 0 до

    ). В частности, для уравнения (13.5) имеем
    0 2
    0 0
    d d
    h
    p
    p
    w
    w
    a
    x t
    t
    x
    x







     








    
    (13.7)
    Рассмотрим интеграл от первого слагаемого (13.7)

     

    1/ 2 1/ 2 1/ 2 1/ 2 0 0 0 0 0
    d d d d d
    h
    h
    h
    j
    j
    j
    j
    p
    x t
    p
    p
    p
    p
    h
    x
    t
    p
    x













    
    

    . (13.8)
    При интегрировании второго слагаемого скорость w принимаем равной своему значению в центре ячейки на нижнем временном слое. Тогда имеем




    1/ 2 0
    1 1/ 2 1
    0 0
    d d d
    j
    h
    j
    j
    j
    j
    j
    p
    w
    x t
    w
    P
    P
    t
    w
    P
    P
    x













    

    (13.9)
    Интегрирование третьего слагаемого дает




    0 0 2
    2 0
    2 1
    1
    d d d
    j
    j
    j
    j
    h
    w
    a
    x t
    a W
    W
    t
    a W
    W
    x








     




    

    (13.10)
    С учетом (13.8)…(13.10) уравнение (13.7) принимает вид




    1/ 2 2
    1/ 2 1/ 2 1
    1
    j
    j
    j
    j
    j
    j
    j
    p
    p
    w
    P
    P
    a W
    W
    h





     




     



    (13.11)
    Аналогично интегрируется и уравнение (13.6)
    2 2
    2 1/ 2 1/ 2 1/ 2 1/ 2 1/ 2 1/ 2 1
    2 2
    2
    j
    j
    j
    j
    j
    j
    j
    j
    w
    w
    W
    P
    W
    P
    w
    w
    h
    d
    w














     






















    .(13.12)
    Таким образом, формулы (13.11) и (13.12) позволяют вычислить ско- рость и давление в центрах ячеек сетки на верхнем временном слое (при
    t = t +

    ) по известным своим значениям на нижнем временном слое (при
    t = t), если известны параметры течения на границах ячейки. Для опреде- ления этих величин используется решение задачи о распаде разрыва гид- родинамических параметров.
    Задача о распаде разрыва. В рассматриваемой разностной схеме ис- пользуется кусочно постоянная аппроксимация параметров потока по дли- не трубопровода. То есть на начало каждого рассчитываемого интервала времени на границах ячейки имеется разрыв параметров. Это разрыв физи- чески не может существовать и должен распасться. Таким образом, если удастся определить, как изменяются давление и скорость при распаде их разрыва, то мы сможем найти необходимые «большие» величины P и W.
    Вновь рассмотрим уравнения (13.5) и (13.6). Оценим порядок их сла- гаемых. Так как, в соответствии с (8.76) изменение давления при гидро- ударе имеет порядок p
    a w

     
    , то справедлива оценка

    210
    p
    w
    w
    aw
    x
    x





    Следовательно, можно пренебречь вторым слагаемым уравнения (13.5) по сравнению с третьим. Аналогично оценивая члены уравнения (13.6), получим, что можно пренебречь кинетической энергией потока под знаком производной. Кроме того, при решении задачи распада разрыва не будем учитывать действие сил вязкости, что позволит получить простые алгеб- раические выражения для «больших величин». После упрощения система
    (13.5), (13.6) примет вид
    2 0;
    1 0.
    p
    w
    a
    t
    x
    w
    p
    t
    x



     

     


    





     
    
    (13.13)
    Умножим первое уравнение (13.13) на 1/ a

    . Полученное уравнение сложим со вторым. Затем из второго уравнения вычтем первое, умножен- ное на 1/ a

    . В результате получается два следующих уравнения
    0
    p
    p
    w
    a
    w
    t
    a
    x






















    ;
    (13.14)
    0
    p
    p
    w
    a
    w
    t
    a
    x






















    (13.15)
    Введем новые переменные
    p
    Y
    w
    a
     

    ,
    p
    Z
    w
    a
     

    . С учетом этих обо- значений уравнения (13.14), (13.15) примут вид
    0,
    0
    Y
    Y
    Z
    Z
    a
    a
    t
    x
    t
    x












    (13.16)
    Эти уравнения можно рассматривать, как производные от Y вдоль на- правления d / d
    x
    t
    a

    и Z – вдоль направления d / d
    x
    t
    a
     
    . Параметры Y и Z
    носят название Римановых инвариантов, а направления d / d
    x
    t
    a
     
    назы- ваются характеристиками системы уравнений (13.16). Равенство нулю производных вдоль указанных направлений свидетельствует о том, что
    Римановы инварианты сохраняют свои значения на соответствующих ха- рактеристиках.
    Рассмотрим две соседние ячейки расчетной сетки, рис. 13.2.
    В момент времени t = 0 гидродинамические параметры на их границах имеют разрыв, который распадается. Условия постоянства Римановых ин- вариантов на характеристиках имеют вид:

    на характеристике d / d
    x
    t
    a

    1/ 2 1/ 2
    j
    j
    j
    j
    p
    P
    w
    W
    a
    a







    ;
    (13.17)

    211

    на характеристике d / d
    x
    t
    a
     
    1/ 2 1/ 2
    j
    j
    j
    j
    p
    P
    w
    W
    a
    a







    (13.18)
    Рис. 13.2. Распад разрыва на границе ячеек
    Из (13.17) и (13.18) находим искомые «большие величины»




    1/ 2 1/ 2 1/ 2 1/ 2 1
    2 2
    j
    j
    j
    j
    j
    a
    P
    p
    p
    w
    w









    ,
    (13.19)




    1/ 2 1/ 2 1/ 2 1/ 2 1
    1 2
    2
    j
    j
    j
    j
    j
    W
    w
    w
    p
    p
    a









    (13.20)
    Таким образом, получаем следующий алгоритм расчета распределения скорости и давления по длине трубопровода:
    1) По известным значениям параметров в середине ячеек на предыду- щем временном слое
    1/ 2
    j
    w

    ,
    1/ 2
    j
    p

    ,
    1,
    j
    N

    по формулам (13.19), (13.20) находим параметры течения на границах ячеек
    ,
    j
    j
    W P
    ,
    1,
    1
    j
    N


    2) По формулам (13.11), (13.12) вычисляем гидродинамические пара- метры течения на верхнем временном слое (при t = t +

    )
    1/ 2
    j
    w

    ,
    1/ 2
    j
    p

    ,
    1,
    1
    j
    N


    Приведенный алгоритм численного решения известен, как метод
    С.К. Годунова [6].
    3) Для того, чтобы найти параметры течения в середине первой и по- следней ячеек на верхнем временном слое
    1/ 2
    w
    ,
    1/ 2
    p
    ,
    1/ 2
    N
    w

    ,
    1/ 2
    N
    p

    необ- ходимо задать граничные условия. Для формирования граничных условий на левой границе канала используем соотношение на характеристике d / d
    x
    t
    a
     
    , исходящей из ячейки с индексом ½, и соотношение, следую- щее из уравнения Бернулли, записанное для сечений в резервуаре и на входе в трубопровод:
    2 0
    0
    вх
    0 2
    p
    P
    W


     
    (13.21)
    Решая совместно (13.21) и (13.18) (при j = 0) находим выражения для
    «больших» параметров на левой границе трубопровода


    2 0
    0 1/ 2 1/ 2 2
    2
    /
    W
    a
    a
    p
    p
    w a
      


     
    ,
    (13.22)
    t
    j
    j+1
    j

    1
    t+

    t
    1/ 2 1/ 2
    ,
    j
    j
    w
    p


    1/ 2 1/ 2
    ,
    j
    j
    w
    p


    d d
    x
    a
    t

    d d
    x
    a
    t
     
    ,
    j
    j
    W P

    212 2
    0 0
    вх
    0 2
    P
    p
    W


     
    (13.23)
    Для формирования граничных условий на выходе из трубопровода до- бавим к расчетной сетке справа одну фиктивную ячейку. Параметры в ней при открытом клапане задаем следующим образом
    1/ 2 1/ 2
    N
    N
    w
    w



    ,
    1/2
    н
    N
    p
    p


    , что имитирует истечение жидкости в атмосферу. После за- крытия клапана эти параметры принимают значения
    1/ 2 1/ 2
    N
    N
    w
    w


     
    ,
    1/2 1/2
    N
    N
    p
    p



    , чем моделируется жесткая стенка. Искомые величины
    1/ 2
    w
    ,
    1/ 2
    p
    ,
    1/ 2
    N
    w

    ,
    1/ 2
    N
    p

    теперь могут быть вычислены с использованием опи- санного в п. п. 1) и 2) алгоритма
    *
    13.4. Аппроксимация и устойчивость численного решения
    Исходя из метода вычисления «больших величин» можно понять, что если вести интегрирование системы (13.5), (13.6) с шагом большим, чем
    *
    /
    h a
     
    , мы получим нефизическое решение. Действительно, из схемы рис. 13.2 понятно, что в момент времени
    /
    t
    t
    h a
      
    пересекутся характе- ристики d / d
    x
    t
    a
     
    , исходящие не из соседних ячеек с номерами
    1 / 2
    j

    и
    1 / 2
    j

    , а из ячеек, располагающихся левее первой и правее второй. В ре- зультате мы получим неверные значения «больших величин» и тем самым будут нарушены законы сохранения (13.6), (13.7).
    Численное решение, полученное при таком шаге интегрирования, в лучшем случае, не будет соответствовать истинному протеканию модели- руемого физического процесса и приведет к неверной интерпретации ре- зультатов расчетов. В худшем случае решение потеряет устойчивость. Ог- раничение на шаг интегрирования
    *
    /
    h a
       
    носит название условия
    Куранта. Выполнение этого условия сопряжено либо с увеличением шага расчетной сетки по пространству h, либо уменьшением шага интегрирова- ния по времени

    . Первое чревато потерей точности решения, второе – увеличением времен, потребного на решение задачи. Число Куранта, оп- ределяемое в вычислительной газовой динамике, как
    *
    Заметим, что если для задания граничных условий на входе в трубопровод доба- вить к нему слева фиктивную ячейку (как это сделано при формировании правых гра- ничных условий) и положить в ней параметры, равные параметрам в резервуаре, ре- зультаты расчета будут существенно отличаться от известных экспериментальных дан- ных. Дело в том, что в такой ячейке нельзя задавать скорость равной нулю, так как фак- тически в прилегающем к трубопроводу объеме резервуара происходит разгон жидко- сти. Как следует из (13.19), такой подход приведет к существенному занижению давле- ния в начальном сечении трубопровода и, как следствие, к значительному уменьшению скорости течения в нем и уменьшению давления гидроудара.

    213


    Ku min
    /
    h
    u
    a







    ,
    (13.24) где операция минимум берется по всей расчетной области, является важ- ной характеристикой эффективности численной схемы. Для явных разно- стных схем условие устойчивости решения имеет вид
    *
    Ku 1

    (13.25)
    В качестве примера, на рис. 13.3 приведены результаты расчетов рас- сматриваемой модельной задачи. Результаты соответствуют одним и тем же моментам времени, но если результаты расчетов давления на рис. 13.3,
    а (Ku = 0,5) качественно верно описывают процесс распространения удар- ной волны вверх по трубопроводу и обратной волны разрежения, то рис. 13.3, б (Ku = 1) не поддается никакому физическому трактованию.
    Причина осцилляций заключается в ошибках округления при работе про- цессора ЭВМ, что приводит к некоторому нарушению физических законов сохранения в результатах вычислений, так как численная схема находится на границе устойчивости. При Ku > 1 амплитуда осцилляций возрастает и решение полностью теряет устойчивость.
    а)
    б)
    Рис. 13.3. Изменение давления при гидроударе,
    /
    t
    ta L

    :
    1 0
    t

    ;
    2 0,065
    t

    ;
    3 0,365
    t

    ;
    4 0,665
    t

    ;
    5 0,965
    t

    ;
    6 1,265
    t

    ;
    7 1,565
    t

    ;
    8 1,865
    t

    ;
    а – Ku = 0,5; б – Ku = 1
    Кроме естественного требования устойчивости численного решения разностная схема, применяемая для решения задачи, должна обладать свойством аппроксимации. То есть при уменьшении шагов сетки по време- ни и пространству разность между дифференциальным оператором исход- ной системы уравнений и ее разностным аналогом, иначе говоря, погреш-
    ность аппроксимации должна стремиться к нулю. Оценим эту погреш- ность для нашей схемы. Для этого разложим разностное решение в ряд
    Тейлора в окрестности точки (
    1/ 2
    j
    w

    ,
    1/ 2
    j
    p

    ). Для скорости имеем выраже- ния
    *
    Неявные схемы сохраняют устойчивость решения и при Ku > 1.
    0 1
    2 3
    4 5
    6 0
    0,2 0,4 0,6 0,8 1
    0
    p
    p
    /
    x L
    1
    t
    2
    t
    3
    t
    4
    t
    5
    t
    6
    t
    7
    t
    8
    t
    -1 0
    1 2
    3 4
    5 6
    0 0,2 0,4 0,6 0,8 1
    0
    p
    p
    /
    x L
    1
    t
    2
    t
    3
    t
    4
    t
    5
    t
    6
    t
    7
    t
    8
    t

    214
     
    2 2
    1/ 2 2
    1/ 2 2
    2
    j
    j
    w
    w
    w
    w
    O
    t
    t







     




    ,
    (13.26)
     
    2 2
    2 1/ 2 1/ 2 2
    2
    j
    j
    w
    w h
    w
    w
    h
    O h
    t
    x










    ,
    (13.27)
     
    2 2
    2 3/ 2 1/ 2 2
    2
    j
    j
    w
    w h
    w
    w
    h
    O h
    t
    x










    ,
    (13.28) где
       
    2 2
    ,
    O
    O h

    – члены разложения второго порядка малости.
    Подставим (13.26)…(13.28) и аналогичные разложения для давления в уравнения разностной схемы (13.11), (13.12), (13.19), (13.20). В результате получим
    2 2
    2 2
    2 2
    2 2
    w
    w
    p
    w
    h
    w
    p
    a
    t
    x
    a
    t
    x






     




     
















    2 2
    1
    w p
    w
    a
    O
    h
    a x x
    x

     




     

      
     
    ,
    (13.29)


    2 2
    2 2
    2 2
    2 2
    2
    p
    p
    w
    p
    h
    w
    p
    w
    a
    a
    a
    O
    h
    t
    x
    x
    t
    x
    x





     



     
     




     










    . (13.30)
    Из (13.29), (13.30) видно, что погрешность аппроксимации разностной схемы имеет первый порядок малости


    O
    h
     
    . При измельчении сетки система (13.29), (13.30) переходит в (13.5), (13.6). Описанная разностная схема называется схемой первого порядка аппроксимации. Существуют схемы второго и более высоких порядков, то есть аппроксимирующих ис- ходную систему дифференциальных уравнений с меньшей погрешностью.
    Вместе с тем, необходимо отметить, что система (13.29), (13.30), соот- ветствующая разностной схеме, все таки отличается от исходной системы уравнений
    *
    . Следовательно, численному решению будут присущи свойст- ва, отсутствующие в рассматриваемом физическом процессе. В частности, как видно из рис. 13.3, а, фронт ударной волны в численном решении по- логий («размазан»), хотя теоретически он должен иметь вид скачка так как разностная схема содержит алгоритм вычисления гидродинамических па- раметров на скачке (распад разрыва). Выясним причину этого факта.
    13.5. Течения с разрывами параметров
    Используя исходную систему уравнений (13.5), (13.6), заменим произ- водные по времени в правых частях (13.29), (13.30) производными по ко- ординате. После преобразований получаем
    *
    В исходной системе правые части равны нулю.

    215 2
    2 2
    / 2 2
    2
    a
    w
    w
    p
    w
    a
    t
    x
    t x
    t









     













     









    2 2
    2 2
    2 2
    1 1 Ku
    2
    ha
    w p
    w
    p
    w
    O
    h
    x x
    a
    a
    x
    x


     







     


     






    ,
    (13.31)
    2 2
    2 2
    2 2
    p
    p
    w
    p
    w
    p
    w p
    p
    w
    a
    a
    w
    w
    t
    x
    x
    x
    x
    x x
    x











     

     



     
     















     

     















    2 2
    2 2
    2 2
    2 2
    / 2 1 Ku
    2
    a
    w
    ha
    w
    p
    a
    w
    O
    h
    x
    x
    x
    x











    






     















    . (13.32)
    По аналогии с уравнением Навье-Стокса члены уравнений (13.31),
    (13.32), включающие вторые производные по координате, можно считать силами вязкости, а множитель


    1 Ku
    2
    s
    ha
     

    (13.33) следует считать некоторой эффективной вязкостью. Величину
    s

    называ- ют схемная вязкость. Именно наличие схемной вязкости приводит к раз- мазыванию фронта волны.
    Как видно из (13.33), схемная вязкость убывает при Ku
    1

    . На рис. 13.4 приведены результаты расчета гидроудара при Ku = 0,995. Ос- тальные данные те же, что и в расчетах, соответствующих рис. 13.3, а.
    Рис. 13.4. Изменение давления при гидроударе, Ku = 0,995
    Приведенные результаты показывают, что, действительно, схемная вяз- кость исчезла, фронт ударной волны в момент ее возникновения скачкооб- разный. По мере развития гидроудара фронт волн несколько выполажива- ется, но это связано с действием коэффициента гидравлического трения

    , имеющегося в математической модели процесса.
    Вообще необходимо отметить, что результаты численного моделирова- ния гидрогазодинамических процессов, где возможно образование разрыва параметров (ударных волн, скачков уплотнения), во многих случаях пред-
    0 1
    2 3
    4 5
    6 0
    0,2 0,4 0,6 0,8 1
    0
    p
    p
    /
    x L
    1
    t
    2
    t
    3
    t
    4
    t
    5
    t
    6
    t
    7
    t
    8
    t

    216 ставляют реальную картину в искаженном виде с «размазанными» разры- вами. Причин здесь может быть несколько. Это и грубая сетка, и примене- ние схем «сквозного счета», базирующихся на выполнении интегральных законов сохранения и не «обрабатывающие» разрывы, и схемная вязкость.
    Кроме того, в разностные схемы второго и более высоких порядков специ- ально вводится искусственная вязкость (псевдовязкость) для подавления осцилляций, так как данным схемам присуща немонотонность решения.
    Поэтому, если задачей расчетов является идентификация разрывов па- раметров, то для этого должны выбираться и соответствующие численные методы.
    13.6. Нефизическое поведение решения
    При численном моделировании возможны ситуации, когда поведение решения не соответствует ожидаемому (исходя из качественных представ- лений о процессе), либо вообще противоречит физическим законам. При- чина этого может, в частности, заключаться в том, что примененная для расчетов математическая модель стала некорректной при реализовавшихся в расчете параметрах. В качестве примера на рис. 13.5 показаны результа- ты расчетов распространения волны разрежения после прохождения удар- ной волны в трубопроводе по математической модели, описанной выше.
    Рис. 13.5. Распространение волны разрежения по трубопроводу:
    9 2,165
    t

    (остальные моменты времени те же, что и на рис. 13.3)
    Как видим в момент относительного времени
    9 2,165
    t

    на конечном участке трубопровода реализовалось отрицательное давление, что в невоз- можно в классической термодинамике. Причина этого кроется в том, что использованная математическая модель описывает только сплошную сре- ду, а в данном расчете в наиболее интенсивной области волны разрежения сплошность нарушается вследствие кавитации жидкости. Для корректного описания процесса математическая модель должна учитывать возможность появления перемещающейся свободной поверхности жидкости с поста- новкой граничных условий на ней, следующих из модели парового объема
    -3
    -2
    -1 0
    1 2
    3 4
    5 6
    0 0,2 0,4 0,6 0,8 1
    0
    p
    p
    /
    x L
    1
    t
    2
    t
    3
    t
    4
    t
    5
    t
    6
    t
    7
    t
    8
    t
    9
    t

    217 кавитационной каверны. Сами уравнения движения жидкости также долж- ны быть изменены для учета потерь массы и количества движения при фа- зовом переходе.
    В настоящей главе затронута только малая часть вопросов, с которыми приходится сталкиваться исследователю при численном моделировании процессов гидрогазодинамики.
    Методы численного моделирования интенсивно развиваются уже более
    60 лет благодаря существенному прогрессу в области создания ЭВМ. В на- стоящее время имеется множество специализированных коммерческих программных продуктов, предназначенных для моделирования гидрогазо- динамических процессов, например, Ansys Fluent, Ansys CFX, Flow Vision,
    Gas Dynamics Tool и др. Задача современного инженера заключается в грамотном применении этих инструментов для решения практических за- дач, для чего и необходимо быть знакомым с принципами получения чис- ленного решения.
    В заключение заметим, что особенности, присущие численным мето- дам, некоторые из которых отмечены выше, свидетельствуют, что, приме- нение несомненно мощных и полезных инструментов численного модели- рования гидрогазодинических параметров требует вдумчивого и аккурат- ного подхода. Начиная работу по численному моделированию надо четко представлять цель проведения расчетов, выполнить приближенные оценки диапазона изменения параметров процесса, базируясь на которых следует выбирать математическую и постановку и численный метод, в наибольшей степени отвечающей задаче моделирования. Особенно аккуратно следует подходить к моделированию малоизученных процессов, к которым можно отнести, например, турбулентные течения в областях сложной геометрии, течения с химическим реакциями и др.
    13.7. Контрольные вопросы
    1. Дайте определение понятию конечно-разностная аппроксимация системы дифференциальных уравнений.
    2. Дайте определение понятию разностная схема.
    3. Что означает термин порядок аппроксимации?
    4. Дайте определение понятиям явная и неявная разностная схема.
    5. Дайте определение понятию число Куранта.
    6. Каково условие устойчивости решения при использовании явной разностной схемы?
    7. В чем причина появления схемной вязкости, каковы формы ее про- явления?
    8. Что означает понятие псевдовязкость?
    9. Что может явиться причиной получения численного решения не со- ответствующего физическим закономерностям?

    218
    1   ...   14   15   16   17   18   19   20   21   22


    написать администратору сайта