Тензорное представление алгоритмов вычислительной гидромеханики

The tensor notation of the computational fluidmechanics algorithms. V.N. Khramushin, A.V.Fain (Special Design Office of the Means of Marine Research Automation, FEB RAS, Yuzhno-Sakhalinsk).

Рождение математических правил построения пространственных объектов и методов анализа их пространственно-временных взаимодействий может быть связано с работами механиков XVII в. Рене Декарта и Исаака Ньютона, определивших законы векторных взаимодействий в абсолютных системах отсчета классической механики.

Декарт, продолжая математизацию и геометризацию физики, определил меру локального движения одного объекта относительно другого: «Если из двух человек один движется с кораблем, а второй стоит неподвижно на берегу ... то нет никакого преимущества ни в движении первого, ни в покое второго». Это означало введение единых векторных законов для описания физических объектов и процессов в абсолютной и в локальных системах отсчета, формально связанных с различными подвижными объектами.

Ньютон вносит основные правила векторной геометрии в механику: «Абсолютное пространство, по самой своей сущности, безотносительно к чему бы то ни было внешнему, остается всегда одинаковым и неподвижным. Относительное - есть его мера или какая-либо ограниченная подвижная часть...». Получается, что абсолютное пространство всегда условно, а физические объекты могут быть математически строго связаны с производными от них локальными системами отсчета. Ньютон доопределяет и главную мировую скалярную величину – время: «Все движения могут ускоряться или замедляться, течение же абсолютного времени изменяться не может»

«Время и пространство составляют как бы вместилища самих себя и всего существующего. Во времени все располагается в смысле последовательности, в пространстве же – в смысле порядка положения». Этих определений пока вполне достаточно для построения новых компьютерных вариантов математического описания элементарных объектов, операций и алгоритмов вычислительной гидромеханики.

Однако мир компьютерных моделей не совсем однозначно отображается в современной математике, представляющей законы движения жидкости в форме систем дифференциальных уравнений в частных производных. Нельзя ли конечноразностные объекты и операции, присущие цифровым ЭВМ, использовать как независимый математический инструмент? Де-факто математические основы такого подхода хорошо разработаны Ю.М. Давыдовым [3] как метод крупных частиц, иногда называемый методом конечного объема.

В завершенном виде вычислительные эксперименты в гидромеханике нуждаются в комплексном представлении информации о больших полях течений жидкости, к которым строго однозначно должны применяться как законы механики жидкости, так и различные алгоритмы для адаптации численных решений к реальным условиям. В такой постановке проектирование вычислительных экспериментов должно опираться на теоретические основы «функционального программирования», в котором исходные данные управляют как применимостью математических функций, так и структурой вычислительных алгоритмов в целом.

Традиционные конечноразностные аппроксимации уравнений движения жидкости, после их алгоритмической формализации, с определенной долей условности могут рассматриваться уже в качестве некой другой гидромеханической модели, качественно отличающейся от исходной записи законов механики в виде системы дифференциальных уравнений в частных производных. Оптимизация же численных моделей для повышения скорости расчетов, контроль физических критериев корректности решений и применение адаптивных методов для расширения области применимости инженерных расчетов нередко выполняются непосредственно в понятиях алгоритмических операций для особых численных объектов внутри ЭВМ. Особо актуален и наиболее сложен такой анализ численных моделей при согласовании физических свойств жидкости в особых условиях, нередко возникающих как внутри расчетной области, так и на свободных границах. Практический опыт реализации вычислительных экспериментов, связанный с решением инженерных задач гидромеханики, требовал практической формализации тех действий, которые необходимы как для ускорения расчетных алгоритмов (включая оптимизацию циклов вычислений), так и для унификации процедур, моделирующих законы гидромеханики на внутренних разрывах и внешних границах расчетной области (обеспечении функциональной невидимости границ). Естественным математическим инструментом для построения элементарных вычислительных объектов стало классическое тензорное исчисление [5], точнее – его малая часть, моделирующая линейные интерполяционные операции с элементарными физическими объектами в размерном виде (не выше тензоров II ранга – для 3-мерного пространства). При этом свертывание численных схем к тензорной форме приводит к своеобразной канонической записи исходных уравнений, которые наиболее быстрым и эффективным образом реализуются на цифровых ЭВМ. Как показал практический опыт проектирования и реализации вычислительных экспериментов, новые модели в тензорном представлении несут в себе очень глубокий физический смысл и вполне могут быть применены для решения определенного круга задач механики сплошных сред, в том числе в аналитической форме.

И все же важнейшим элементом рассматриваемой работы является поиск универсального математического инструмента для построения вычислительных экспериментов в гидромеханике, который имел бы более глубокие и полные аналитические возможности и одновременно позволял автоматически применять различные физические законы в зависимости от критериев развития процессов в сплошных средах, допускал бы сквозное распараллеливание всех вычислений вплоть до каждой элементарной частицы жидкости.

Вводные определения. Вычислительные эксперименты в механике сплошных сред опираются на собственный математический аппарат конечных разностей, в основе которого лежат аппроксимации непрерывных функций ограниченным количеством опорных величин, распределенных в дискретных ячейках памяти цифровой ЭВМ [3]. Дифференциальное исчисление, основанное на анализе бесконечно малых величин и справедливое для аналоговых моделирующих устройств, уже не в полной мере применимо для цифровых моделей дискретно аппроксимированных физических полей и разорванных расчетными интервалами времени кинематических процессов.

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

Возможность численного решения прикладных задач гидромеханики практически всегда усложнена аппроксимационными ограничениями или невыполнимыми критериями, корректный учет которых возможен только при использовании специальных адаптивных решений и гибридных счисленных схем, основанных на эвристических, эмпирических или экспериментально поверенных зависимостях. Для сохранения принципиальной возможности использования гибридных алгоритмов, а также с целью сквозного контроля качества выполняемых расчетов при проектировании вычислительных экспериментов очень важно поддерживать физическую непротиворечивость всех вычислительных объектов и численных операций с ними [6]. Тогда: 1) элементарные объекты гидромеханики должны определяться только в размерной форме, в строгом соответствии с заданным набором свойств жидкости; 2) все вычислительные операции должны соответствовать конкретным законам гидроаэромеханики, а использование промежуточных аналитических преобразований для упрощения исходных математических моделей - недопустимо; 3) расчетные алгоритмы должны автоматически отслеживать принадлежность численных компонент к единой абсолютной и связанным локальным системам координат, обеспечивая необходимые переходы между ними, с безусловным контролем соответствия всех физических размерностей.

Соблюдение столь жесткого комплекса логических правил заметно ограничивает использование новых алгоритмов для прямого воспроизведения многих традиционных задач вычислительной гидромеханики, и особенно – основанных на дифференциальных уравнениях высокого порядка. Однако эти алгоритмы не содержат принципиальных ограничений для прямого численного моделирования задач механики сплошных сред любой сложности. А некоторые из реализованных вычислительных экспериментов иногда лучше сравниваются с натурными наблюдениями, так как аналитические решения и традиционные численные методы пока не позволяют получать аналогичных качественных результатов.

Основные предпосылки для построения проекта вычислительного эксперимента. Традиционная физическая интерпретация, так же как и численное решение задач гидромеханики, опирается на два взаимосвязанных и существенно различающихся по форме геометрических представления о течении жидкости: это известные дифференциальные формализации течений по Эйлеру – на фиксированной сетке и изменяющихся в ее ячейках значениях гидродинамических полей, и по Лагранжу – с фиксацией законов сохранения внутри каждой из движущихся частиц жидкости, которые тоже образуют своеобразную подвижную и нерегулярную сетку.

Наиболее эффективные вычислительные эксперименты в гидромеханике основаны на методах разделения решений по физическим процессам, в которых этап моделирования с использованием Лагранжевых свободных частиц жидкости сменяется процессом согласования физических полей на неподвижных Эйлеровых узлах расчетной сетки. Формально это два этапа интегрирования первого порядка, которые совместно воссоздают решение традиционного дифференциального уравнения механики второго порядка.

Лагранжевы свободные частицы в пределах собственных геометрических размеров должны вести себя подобно энергоемким объектам, которые во взаимодействии с окружающей средой частично растрачивают энергию собственного движения, частично аккумулируют ее внутри себя (рис. 1).

Рис. 1. Элементарный объект вычислительной гидромеханики – крупная частица жидкости, будучи всегда окруженной непрерывной средой с подобными свойствами, реагирует на суммарное внешнее воздействие с учетом предыстории ее несвободных движений и внутренних деформаций

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

Прописными буквами обозначаются величины, измеренные в абсолютных системах отсчета (R_k. – вектор координат для точки в пространстве, Т - отсчет времени от начала вычислительного эксперимента). Строчными - измеренные внутри локальных (подвижных) базисов (r^k – вектор локальных координат точки, определенной в связанной системе отсчета внутри конкретной частицы жидкости). Тогда в пространственных обозначениях прописными буквами показываются большие по величине отсчеты, соизмеримые с масштабом расчетной области, строчными - конечные разности, которые всегда используются в качестве линейных дифференциалов.

Левые верхние индексы будем отмечать текущие отсчеты во времени, к которым относятся вычислительные объекты или процессы. Индекс + - это следующий отсчет времени. Если время в вычислительном эксперименте связано с расчетными тактами k, то: ⁺F = ^k+1F, что соответствует отсчету времени: ⁺T=^k+1T=⁰T+(k+1)·t, где t - расчетный интервал времени. Левые нижние индексы - номера или индексы пространственных ячеек сеточной области, которые можно трактовать как неподвижную Эйлерову, или абсолютную, систему координат. Исходное местоположение крупной частицы жидкости размечаем левыми нижними индексами прописными буквами { J,K }, пространственные метки подвижных частиц жидкости - строчными индексами { j,k }. Нижним индексом + обозначим некоторую смежную расчетную ячейку.

Все правые индексы являются тензорными идентификаторами и обозначают принадлежность к системе отсчета. Нижние индексы привязывают объект к абсолютной системе координат, верхние, соответственно, к локальному базису. Векторные стрелки могут заменять индексы, при этом стрелка вправо ® отмечает принадлежность к абсолютной системе отсчета (нижний индекс), влево ¬ к локальной (верхний индекс). Тензорные величины отмечаются двойными индексами или уголками: Ù - абсолютная; Ú - локальная; >,< и ´ - смешанные системы отсчета.

Местная система координат математически формализуется в виде локального базиса, образуемого некомпланарными векторами с условно единичной длиной: = _i= r_ij, который основан на смежных узлах расчетной сетки, образующих базисные векторы элементарного пространственного симплекса. Новый пространственный объект удобно называть геометрическим или базисным тензором второго ранга, или просто – тензором.

По аналогии, в местной системе отсчета единственным образом могут быть представлены проекции единичных ортов абсолютной системы координат, что образует дуальный базис типа: = ^j = r ^{ij -1}, где компоненты тензора могут быть сформированы простыми геометрическими построениями свободных единичных отрезков абсолютного базиса в отсечках по сторонам косоугольных параллелограммов внутри локального базиса (рис. 2).

Рис. 2. Локальный базис местной системы координат образуется тройкой базисных векторов – ортов условно единичной длины

Смешанные тензоры образуются аналогично, и обычно они проявляются в операциях определения законов механики сплошной среды. К сожалению, при построении тензорных величин для гидродинамических объектов невозможно воспользоваться привычным аппаратом матричной алгебры, ее описательных возможностей не хватает для определения смешанных тензоров. Все правила тензорной алгебры [4] справедливы в полной мере.

Правило суммирования. Если один и тот же индекс встречается в явно одночленном выражении дважды – один раз как верхний, а другой как нижний, то предполагается суммирование соответствующих выражений при пробегании этим индексом значений 1, 2, 3. Такой индекс называют немым.

Проверка на согласованность. Любой индекс, не являющийся немым, должен встречаться во всех членах уравнения либо как верхний либо как нижний индекс. Любое уравнение, содержащее индексы отличные от немых, должно быть интерпретировано как система уравнений, когда каждый индекс принимает значения 1,2,3.

Заметим, что верхние индексы, встречающиеся в знаменателе одночленного выражения, эквивалентны нижним индексам этого выражения.

В индексной форме операция полного тензорного умножения: = · может представлена в следующей форме: a_i =r_ij · a^j, где индексы, стоящие при сомножителях справа на разных уровнях, называют немыми, а их появление в операциях произведения подразумевает суммирование векторных компонент по этому индексу:
a_i = е_j (r_ij · a^j). (9’)

Операции умножения могут выполняться с повышением или понижением ранга результирующего произведения. Так, например, в результате умножения вектора на тензор: (r_ij · a_k) формируется тензор третьего ранга - с тремя индексами. Такие операции не имеют особого физического смысла и потому в настоящей работе не рассматриваются. Более того, они запрещены к использованию в обычных операциях умножения даже в качестве промежуточного результата алгебраических преобразований. Умножение с понижением ранга тензора, типа скалярного произведения векторов: q=a_i · bⁱ (проекция одного вектора на другой), обычно используется только при анализе свойств вычислительных объектов и также запрещено к использованию в операциях умножения еще и потому, что понижение ранга результата означает потерю информации об исходных векторных объектах. Запреты на изменение ранга при выполнении любых операций произведения иногда позволяют избежать логических ошибок при неудачной расстановке тензорных и векторных сомножителей в одночленном выражении.

Жесткие ограничения на применимость операций сложения и умножения для рассмотренных в этой работе векторных и тензорных объектов необходимы для построения специальных вычислительных алгоритмов, которые по формальному признаку принадлежности вычислительных объектов к конкретным системам отсчета позволят автоматизировать все координатные преобразования как при непосредственном решении задач механики, так и при представлении их результатов на экране ЭВМ.

Левые индексы удобно использовать для привязки вычислительных объектов (крупных частиц жидкости): в пространстве – левый нижний индекс, и во времени – левый верхний индекс. Отсутствие левых индексов подразумевает, что определение относится к текущему моменту времени и ко всем элементарным ячейкам расчетной области, соответственно.

_Ω^T – координаты узловой точки, где R указывает на отсчет от нуля абсолютной системы координат, стрелка вправо ® на принадлежность к ней, W определяет местоположение узла в расчетной области; T – момент времени от начала проведения вычислительного эксперимента.

₊^t– ссылка на смежную точку, относительно _Ω^t  точку, смещение в сторону (+) с учетом базисной точки элементарного объема жидкости Ω (крупной частицы) за время t. Стрелка влево ¬ показывает, что «длинный» вектор абсолютных отсчетов  переведен в проекцию локального базиса, но при этом он не был связан с текущим местоположением крупной частицы жидкости.

Также в виде конечных разностей (рис. 3) вводится динамическое понятие приращения скорости: = ^D = ⁺ - ⁰ . Реологические свойства жидкости в вычислительном эксперименте всегда проявляются в виде тензорных описателей внутренних свойств элементарных расчетных объектов, для которых уже в ходе вычислений должно обеспечиваться сохранение реологических параметров течения жидкости, параллельно с проведением контроля критериев существования самого численного решения.

Построение базовых операций вычислительной гидромеханики. Основные законы гидромеханики в тензорном виде позволяют определить не только кинематику элементарных частиц жидкости, но и их внутреннюю деформацию в строгом соответствии с напряжениями от внешних силовых воздействий. Пусть в качестве исходных данных для построения вычислительного эксперимента в гидромеханике будут определены следующие размерные поля в абсолютной (глобальной) системе координат:{}[м] – поле координат узловых точек; {} [м/с] – поле скорости; {}[кг] – тензорное поле внутренних свойств для каждой из частиц жидкости.

В определенной таким образом сеточной области вводятся локальные вычислительные объекты - как элементарные крупные частицы жидкости, в которых физические параметры течения жидкости определяются в той же абсолютной системе координат: [м³] – тензор формы крупной частицы; [м³/c] – тензор локальных (конвективных) скоростей; [Н·м²] – тензор напряжений.

Рис. 3.Тензор локальных скоростей образуется деформационными смещениями базисных векторов крупной частицы жидкости за расчетный интервал времени

Законы механики для деформируемой жидкой частицы в тензорном виде получаются не сложнее обычных арифметических операций, которые удивительно похожи по форме и по сути как для внешних сил, так и для внутренних деформаций. Это позволяет на всех этапах численного решения оперировать довольно простыми аналитическими зависимостями, которые к тому же учитывают эффекты аккумуляции и последующего излучения механической энергии, полученной в результате внешних взаимодействий между смежными расчетными ячейками.

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

1. Векторный закон Ньютона для деформируемой частицы: = · =· ·. (1)
2. Тензор вязких напряжений для Ньютоновой жидкости^**: _Н = ·_H = · _H·. (2)
3. Тензор упругих напряжений для твердого тела Гука: _Г = ·(+ _r·t)· =·( + _r·t), (3)
где тензор локальный скоростей определен алгоритмически как:   = _+i - _0i (рис.3); а главные реологические константы формально представлены как тензоры, для которых справедливы выражения типа: = M _i ^j = · = r_ij·r ^jk - тензор инерции в проекции на абсолютную систему отсчета, и определения: - тензор плотности и аккумуляции внутренней энергии; , - тензоры динамической вязкости и жесткости жидкости.

Используем разложение тензора конвективных скоростей на девиаторный ^* и шаровой тензоры как: = ₀ + ^*. Тензор шарового сжатия ₀ представляется диагональной матрицей, все элементы которой равны между собой, в тензоре конвективных скоростей он характеризует скалярную дивергенцию (сжатие) потоков жидкости и связан с первым реологическим параметром жидкости: – коэффициентом динамического сжатия.

Девиаторный тензор ^* имеет нулевой линейный инвариант:tr(^*)=0, что формально свидетельствует о свойствах несжимаемости оставшегося потока жидкости. Тензор вязких напряжений получается выделением кососимметричной части из матрицы девиаторного тензора *, что формально соответствует вращению элементарных частиц жидкости и связанному с ним внутреннему напряжению: = · _H = · ( ^*- ^{* T})/2. Оставшийся симметричный тензор будет определять чистую деформацию элементарной частицы жидкости, которая может быть связана с  упругими напряжениями: _Г = · _Г · t = ·( ^*- ^{* T})·t/2. Соответственно, полный тензор внутренних напряжений: = (· ₀+ · _Г)·t + · _H.

Обычно тензорные определения реологических параметров применяются только при анизотропных свойствах сплошных сред, но в случае обсуждаемых здесь вычислительных моделей традиционно скалярные физические параметры жидкости должны быть представлены в качестве ортогональных тензоров, определяемых с помощью соответствующих диагональных матриц.

С тензором плотности связывается предыстория деформации элементарных расчетных объектов, что на Лагранжевом этапе вычислительного эксперимента позволяет считать крупные частицы жидкости участвующие в свободном (корпускулярном) движении. Кроме тензора плотности на всех расчетных этапах будут проявляться упругие и вязкостные свойства жидкости в виде связанных тензоров и , реальные физические величины которых зависят от текущего состояния тензора формы . Если же в вычислительном эксперименте используется модель идеальной жидкости на неподвижных узлах сетки _Ω, то суть вычислений будет сводиться к обеспечению изотропности среды, когда тензоры  и образуются только равными по величине диагональными компонентами.

Свойства вычислительных операций и элементарных объектов. Все операции вычислительной гидромеханики выполняются исключительно в размерной форме, что необходимо для контроля физической корректности моделируемых процессов для каждой отдельно взятой частицы жидкости и, при необходимости, использования гибридных схем в подобластях, где происходит изменение режимов течений и автоматическое изменение используемых законов механики и их математических моделей. С вычислительными объектами связываются три группы операций:

1) операции логические, или эмпирические. Это физические законы, такие как: если есть тензор конвективных скоростей , то он порождает тензор напряжений по закону (2), и т.п. Такие операции задают способ формирования и методы анализа состояния тензорных объектов, а также принятие решений об изменении математических моделей по ходу расчетов;

2) операции сложения, применяемые только к величинам с одинаковыми физическими размерностями и при условии, что они определены и построены в одном и том же базисе. К сложным объектам сложение может также применяться в условиях, когда определены операции автоматического согласования систем отсчета и физических размерностей;

3) операция "произведение" применима только к объектам, у которых связываемые компоненты лежат в дуальных системах отсчета. Запрещается повышение ранга тензора в операциях произведения. Для построения векторных или тензорных объектов должны определяться специальные алгоритмы, опирающиеся на физическую постановку задачи о проведении вычислительного эксперимента.

Вычислительные объекты наделяются следующими свойствами: 1) скалярные или инвариантные величины, например время t, могут участвовать только в операциях произведения, с любыми объектами; 2) векторные величины участвуют в операциях сложения с сопоставимыми векторами и в операциях произведения с тензорами, если это необходимо для перехода из одной системы отсчета в другую. Применение умножения векторов для построения диад недопустимо; 3) тензорные величины определяют физические свойства элементарных частиц жидкости, их геометрические деформации [1], кинематические свойства и другие динамические процессы в сплошной среде. Тензоры могут участвовать как в операциях сложения, с сопоставимыми тензорами, так и  в операциях произведения со связанными векторными или скалярными объектами. С использованием тензорных объектов формулируются основные законы гидромеханики, а автоматическое конструирование и анализ этих объектов образует собственно и саму логико-математическую задачу о проведении вычислительного эксперимента в гидромеханике.

Движение контрольной точки в пространстве, сместившейся за один такт первого (Эйлерового) этапа вычислительного эксперимента, может быть представлено как пересчет координат точки в локальной системе отсчета в абсолютную, с учетом интервала времени t :
⁺ = + ·t + ··t²/2 + ( + ·t + ··t²/2)· ... ( 4 )

Следующим этапом вычислений может быть определение вектора приращения локальной скорости для этой точки (ускорения), по полю которого затем восстанавливается новое поле внешних сил и внутренней энергии для исходных крупных частиц жидкости:
⁺ = + ( + ·t )·. _{( 5 )}

Таким образом выполняется первый шаг вычислительного эксперимента, за которым следует контроль корректности решения на поле { r }; согласование условий сохранения физических параметров внутри элементарных ячеек жидкости { M }; а в завершение – перерасчет новых для каждой крупной частицы жидкости «внешних сил»:
= · · = · · · .

Схемой расщепления решения для уравнений второго порядка по физическим процессам будем называть деление процесса вычислений на два этапа, на каждом из которых выполняется дифференцирование первого порядка. Если на одном из этапов выполняется пространственное дифференцирование поля скорости, то затем по нему вычисляется новое поле напряжений, что позволяет на втором этапе выполнить расчет смещений частиц жидкости по скорректированному ускорениями полю скорости и с использованием конвективных перераспределений потоков согласовать условия сохранения, и снова восстанавливается поле скорости, но уже на новый момент времени.

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

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

Процесс трансформации полей, в случае исходного требования сохранения координат узловых точек {}, может быть разделен на три этапа:

I. Кинематика. Новое поле узловых точек:
{⁺ = + ·t + ··t²/2 }_{( 6 )}

Поле конвективных скоростей образуется алгоритмическим построением тензора:
{ } = { _+i - _0i } _{( 7 )}

Расчетное состояние нового поля внутренних свойств:
{ ⁺ } = {·+ } = { ( + ·t)· };
{⁺ } = { ·( +·t) } ... ( 8 )

II Динамика. Здесь производится сопоставление реологии жидкости с текущим состоянием вычислительной модели. Пусть определен закон сохранения количества движения на разнесенном по этапам вычислений интервале времени:
⁺· ₌ ·( + ^D ), исходя из (8):
^D = ( ⁺ - ) · · = ··t,  ( 9 )
получается уравнения Ньютона в форме Эйлера, справедливое для крупной частицы на неподвижных узлах расчетной области:
= ·· = ···. ( 1* )

Полученное выражение содержит в себе напряженное состояние жидкости, которое можно объяснить реологическими свойствами вычислительной модели течения:

{ } = { _+i - _0i } или = ··· = ··.

По форме новое уравнение соответствует записи напряжений в уравнениях Навье-Стокса. Реология реальной жидкости оформляется в виде законов (2),(3), связанных с тензором конвективных скоростей . Тогда суть расчетов может быть сведена к согласованию свойств вычислительной и физической моделей жидкости.

В результате разностного дифференцирования в выражении (9) была утеряна объемная составляющая ускорения, которая может быть получена при рассмотрении движения частицы с переменной массой, без учета деформации:
M = det( ), r = det( ),
  ₍₁₀₎

Рассматриваемая вычислительная модель всегда наделяет жидкость свойствами сжимаемости, вязкости и упругости [1]. В полученных уравнениях пока не рассматривались деформации тензора формы крупной частицы , что позволяет считать жидкость внутри нее изотропной. Выделением диагонального тензора ₀, такого, что след от остатка будет равен * будет равен нулю: = ₀ + ^*   (tr* =0), получим тензор шарового сжатия:
₀ = ·₀·t, _{( 11 )}
где компоненты тензора ₀  определяют давление, – коэффициент динамического сжатия.

Выделением кососимметричной части матрицы *, которая задает вращение одной крупной частицы относительно смежных, будет получен тензор вязких напряжений:
_Н = ·_H = ·( ^{* - * T})/2. _{( 12 )}

Оставшийся симметричный тензор связывается с упругой деформацией (6):
_Г = ·_Гt = ·( ^*- ^{* T})·t/2. _{( 13 )}

Полный тензор внутренних напряжений: = (·₀+·_Г)·t + ·_H. _{( 14 )}

Рассматриваемые динамические коэффициенты  , и отличаются от кинематических, коэффициентом скалярной плотности r. Под действием тензора напряжений , частица получает приращение скорости внутреннего (замкнутого) движения:
^D = · t /r _{( 15 )}

Если течение установившееся, то, за расчетный интервал времени, тензор приращения скоростей ^D должен компенсировать сам тензор конвективных скоростей: · + · ^D = 0. Это выражение на данном этапе вычислений является точным, поскольку не учитывается смещение крупных частиц за время t.

Последнее выражение является точным, если не учитывается смещение крупных частиц за время t.

3. Статика. На завершающем этапе необходимо провести восстановление поля скорости по вычисленным на втором этапе приращениям, при этом рассматриваются деформационные движения вокруг статических центров тяжести частиц.

От тензора ^D , определенного в локальном базисе, необходимо перейти к смешанному тензору, где локальные характеристики опираются на абсолютную систему отсчета:
^D = ^D ·. _{( 16 )}
Этот тензор не связан с локальной геометрией базиса .

Для перехода к исходной сетке построим новый локальный базис, который опирается на неподвижный узел, а в качестве смежных использует смещенные во времени пространственные точки:
⁺ = ₊⁺_i - ₀⁰_i, ( 17 )

Раскрывая выражение (7), использованное при построении тензора локальных скоростей: :    = _+i - _0i, по скоростям, связанным с новыми базисными векторами , получим алгоритм вычисления нового поля скорости:
⁺ = + S_i·^D_i, _{(18 )}
который суммирует приращения скорости от окружающих крупных частиц.

Выражения (15) - (18) раскрывают основные алгоритмические построения, позволяющие применять обратный закон Ньютона: = · .

Ускорения, полученные в векторной форме: для частицы с переменной "массой" (10); и для распределенных массовых сил, должны быть интерполированы с центров ячеек на исходные узлы расчетной области.

Словесно этапы вычислений определяются следующим образом:
I - на неподвижной Эйлеровой сетке производятся расчеты распределенных характеристик течения;
II - рассматриваются внутренние свойства частиц жидкости с целью построения тензоров "массы", в которых сохраняется предыстория деформации. Здесь же организуются итерационные процессы установления, в которых происходит согласование реологии вычислительной и физической моделей течения жидкости;
III - характеристики течения  интерполируются со смещенных в Лагранжевом движении центров тяжести крупных частиц на исходные узлы расчетной области. На этом же этапе можно рассмотреть условия на свободных (жидких) границах, где вместо алгоритмов интерполяции в разорванных узлах нерегуляризованной сетки внутри расчетной области будет использоваться экстраполяция с помощью заграничных центров особых (граничных) частиц жидкости.

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

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

Контроль корректности выполняемого решения сводится к контролю допустимости интерполяционных формул. В качестве примера разрушения решения можно привести случай, когда одна грань частицы, в результате свободного движения в заданном поле скоростей, обгонит противоположную грань этой же частицы. Этого можно избежать уменьшением расчетного интервала времени. Выражение, показывающее такой предельно допустимый интервал, обычно называют критерием устойчивости решения во времени. Известный критерий Куранта в точности соответствует вышеуказанному определению.

С контролем пространственной аппроксимации дело обстоит существенно сложнее, так как соответствующие критерии связываются не только с начальными условиями, но и с локальными геометрическими характеристиками расчетной области. На отладочных этапах реализации вычислительных экспериментов можно предложить только метод визуального контроля формы расчетных полей, прорисовываемых на графическом экране ЭВМ после выполнения каждого такта вычислений.

В соответствии с выражением (8): ⁺ = ( +·t )· , масса частицы фактически является сумматором тензоров локальных скоростей. На первом этапе используется для определения внутренних свойств частицы, на третьем - с помощью = ^-1 описывается взаимодействие смежных частиц жидкости. Вычислительная модель всегда содержит в себе три основных свойства жидкости: сжимаемость; вязкость; упругость, при этом соотношение интенсивностей указанных напряжений, может привести к критическому режиму или изменить режим течения жидкости (имеется в виду образование струй, вихревых слоев и кавитационных разрывов).

= ·₀·t·+ ·_Н·+ ·_г·t·= ₀ + _Н + _Г. _{( 14}^*₎

Тензор напряжений (14) можно определить в качестве характеристического полинома для внутреннего состояния расчетной частицы, в котором реологические параметры жидкости встанут на место главных инвариантов тензора:
₀ : I ¹ 0 - сжимаемость; _Н : II ¹ 0 - поворот; _г : III ¹ 0 - чистая деформация, другие инварианты этих тензоров равны нулю.

Если внутреннее состояние частицы рассматривается без учета смещения расчетных узлов, то можно поставить требование компенсации тензора локальных скоростей за счет приращений скорости вызываемых тензором напряжений:
·t + ··t²/2 = 0, ( 19 )
из которого, в качестве сумматора напряжений можно вывести тензорную плотность:
⁺ = + ( ₀ + _Н + _Г.)·t²/2 @ ₀ + _Н + _Г. _{( 20 )}

Рассмотрим два варианта разрушения тензора ⁺ при выполнении одного такта вычислений,  когда r = det( ) обращается в ноль.

_Н - не рассматривается;
₀ такой, что  det( ₀ ) < 0 - соответствует полю с разрежением;
* такой, что det( ₀+ _г ) = 0 - тензор плотности обращается в диаду, которая может быть определена плоскостью, перпендикулярной к главной оси растяжения тензора упругости _г. Если жидкость не выдерживает отрицательного давления, то данная плоскость должна быть использована в качестве свободной границы, проходящей через крупную частицу.

2. Образование свободной струи или турбулентного вихря.

На завершающем этапе вычислений исходная вихревая частица может обратить в нуль детерминант внутреннего поля сопряженной частицы, где r_o > 0. В этом случае можно воспользоваться эмпирическим алгоритмом переноса вихря с первого этапа вычислений в тензор "массы" третьего этапа. Тогда будет снято напряжение со смежной ячейки, а частица в свободном движении будет проявлять особенности, свойственные турбулентному вихрю. Приведенный алгоритм является эмпирической моделью, которую можно использовать за пределами аппроксимационного разрешения расчетной области, для получения точного решения сетка в районе вихря может быть сгущена.

Тензорная арифметика, использованная для проектирования полных и прямых вычислительных экспериментов в гидромеханике, в общем случае всегда имеет дело с реальными жидкостями, обладающими полным комплексом реологических свойств, которые задаются линейными/обратными зависимостями в пределах конечноразностного масштаба элементарных расчетных ячеек.  Основные иллюстрации выполнялись на примере решения гиперболической задачи о распространении упругих продольных волн (длинных волн в океане), которые допускают явные пространственные решения и потому очень удобны для физической интерпретации основных вычислительных процессов.

1. Реализация вычислительных экспериментов о распространении длинных волн в океане, построенная на основе тензорной модели, при реальном моделировании длинноволновых процессов в океане (приливы, штормовые нагоны, цунами) показывает хорошее соответствие модели натурным наблюдениям. Тензорная модель обладает ярко выраженными дисперсионными свойствами (увеличение периода при движении волны над неровным дном, или так называемое красное смещение); показывает хорошее согласование с наблюдениями за временем затухания длинноволновых процессов вблизи побережья (декремент затухания характеризует корректность условий сохранения энергии длинноволновых процессов); в случае успешного согласования аппроксимационных критериев получается хорошее соответствие формы и порядка следования длинноволновых импульсов моделируемого и реального цунами вблизи мелководных участков побережья (корректное протекание интерференционных, дисперсионных и рефракционных процессов для волновых процессов). В тензорной модели всегда содержатся «конвективные» особенности течения жидкости, что для традиционных численных методов связывается с необходимостью моделирования течений с «нелинейными свойствами». Тензорная модель легко позволяет задействовать гибридные схемы и внедрять в решение различные эмпирические зависимости, типа приливного геопотенциала, придонной вязкости и поверхностных напряжений при прохождении штормов и глубоких циклонов. Так как модель основана на принципах разделения решений по физическим процессам, то все инородные зависимости могут включаться или исключаться в процессе проведения вычислительных экспериментов.

2. Важным подтверждением корректности тензорных уравнений служит вычислительный эксперимент по моделированию обтекания нескольких тел произвольной формы потоком идеальной жидкости (т.е. потоком, внутри которого отсутствуют вихреисточники). В этом случае не только показанные здесь тензорные зависимости справедливы по отношению к смежным расчетным ячейкам, но и все дальние расчетные ячейки имеют активные связи, осуществляемые строго по прямым линиям. Это означает, что моделирование обтекания тел идеальной жидкостью можно свести к простому интегрированию движущегося тела по его объему и смоченной поверхности. Для этого объем твердого тела заменяется эквивалентными диполями, а поверхность - вихреисточниками. Такая задача решается успешно, в том числе для нескольких движущихся тел одновременно.

3. В современной вычислительной математике важное место занимает визуализация исходных данных и наглядный анализ результатов вычислительных экспериментов, при котором многомерные зависимости представляются на двумерном поле графического экрана ЭВМ. Показанные в данной работе численные алгоритмы являются простыми и эффективными методами пересчета векторно-тензорных величин между различными системами координат.

Рис. 4. В управлении подвижными объектами всегда используются многократные преобразования систем координат, связанных как с реальными объектами, так и с их представлением в счетно-решающих устройствах и терминалах ЭВМ

Во многих случаях показанные в работе тензорные преобразования активно используются для быстрого и наглядного построения двух- трехмерных интерполяционных процедур, работающих как преобразование систем координат, и в том числе – неортогональных, заведомо обладающих наивысшей численной эффективностью.

Последний пример (рис. 4), в случае использования произвольных тензорных базисов, может служить своеобразной моделью дальнего взаимодействия между энергонапряженными корпускулами – квантами, что может быть с успехом использовано для формализации многих законов электромагнитных взаимодействий [2].

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

Литература
1. Астарита Дж., Маруччи Дж. Основы гидромеханики неньютоновских жидкостей, М.: Мир, 1978. 307 с.
2. Астахов А.В., Широков Ю.М. Курс физики. Т.2. Электромагнитное поле. М.: Наука, 1990, 384 с.
3. Белоцерковский О.М., Давыдов Ю.М. .Метод крупных частиц в газовой динамике. М.: Наука, 1982. 370с.
4. Кильчевский Н.А. Элементы тензорного исчисления и его приложения в механике. М.: Гос. изд-во техн.-теорет. лит.,1954. 167 с.
5. Мак-Коннел А.Дж. , Введение в тензорный анализ, М.: Наука, 1963, 411с.
6. Храмушин В.Н. Трехмерная тензорная математика вычислительных экспериментов в гидромеханике // Высокопроизводительные вычисления и их приложения: Материалы Всерос. науч. конф. М.: НИВЦ МГУ, 2000. С. 114-117.