Основы моделирования динамики технических объектов в программе ПА10

УДК 681.3

, канд. техн. наук, доцент,

, канд. техн. наук, доцент,

студент,

МГТУ имени

Основы моделирования динамики технических объектов в программе ПА10

Можно выделить два подхода к моделированию динамики технических объектов. Первый подход предполагает разработку специализированных программ моделирования для конкретных предметных областей: PSPICE в электронике; ADAMS, DADS в механике и гидравлике; МВТУ для моделирования систем управления, в первую очередь, ядерными реакторами и т. п. Второй подход предполагает разработку универсальных программ моделирования технических объектов разной физической природы на основе метода аналогий переменных и уравнений математических моделей таких объектов: программы ПA6, ПA7, ПA8, ПA9 [1].

Поскольку в настоящее время встают задачи комплексного моделирования сложных технических систем, а сложные технические системы, в свою очередь, состоят из физически разнородных технических объектов, то актуальность второго подхода значительно возрастает. Однако имеющиеся в настоящее время соответствующие программы имеют ряд принципиальных недостатков. В первую очередь это интеграция математического ядра программ с функциональной частью, формирующей дифференциальные уравнения моделируемых систем, а также ограниченность выбранного базиса переменных моделирования [2, 3], что не позволяет практически реализовать новые, неявные методы интегрирования дифференциальных уравнений с порядком точности выше второго [4-7]. В данной статье предлагается и обосновывается новый, наиболее общий подход к моделированию динамики сложных технических систем на основе метода аналогий и обобщенных формальных схем, а также новый обобщенный базис переменных моделирования таких схем. Разрабатываемая на основе этих принципов программа ПA10 будет базовым математическим ядром системы моделирования динамических процессов нового поколения, а ее разработка будет основываться на опыте создания упомянутых выше программ и будет иметь следующие принципиальные отличия от этих программ:

- использование нового обобщенного базиса переменных модели;

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

ВВЕДЕНИЕ

Системы инженерного анализа проектных решений – Computer Aided Engineering (CAE) Systems используются в CAD/CAE/CAM системах, в основном, для решения задач моделирования и анализа объектов производства. При этом традиционно наибольшее внимание уделялось программному обеспечению моделирования и анализа технических объектов на уровне дифференциальных уравнений в частных производных, т. е. таким системам как MSC. Nastran, ANSYS, COSMOS и т. п., которые выполняют моделирование с очень высокой степенью точности и достоверности, но требуют значительных затрат времени и вычислительных ресурсов для моделирования. В настоящее время непрерывно растут требования со стороны организаций, выпускающих конечную продукцию, по снижению стоимости развертывания производства изделий и сокращению сроков выпуска изделий на рынок при одновременном повышении качества и снижении стоимости продукции. Возрастает роль быстрой и точной оценки принимаемых проектных решений. При этом крайне важна именно требуемая точность оценки, поскольку стоимость ошибок в принятии проектных решений на ранних стадиях проектирования чрезвычайно высока. В этих условиях центр тяжести в CAE системах перемещается на системы моделирования и анализа технических объектов на уровне обыкновенных дифференциальных уравнений, типа MSC. ADAMS, обеспечивающих как быстроту, так и точность оценки принимаемых проектных решений. Рынок подобных систем непрерывно рос на 8 – 9 % в год даже во время экономического кризиса (при общем уменьшении рынка CAD/CAE/CAM систем) и предполагается, что этот показатель возрастет до 14 % к 2007 году (по данным Daratech).

Среди подобных систем моделирования наибольший интерес представляют системы моделирования динамических процессов в технических объектах, предназначенные для анализа процессов в различных технических системах во временной и частотной области. Кроме того, подобные системы вместе с системами конструкторского проектирования (CAD) позволяют быстро создавать виртуальный или электронный прототип изделия, используемый для достаточно точной верификации проектных решений. Системы моделирования динамических процессов в технических объектах прошли достаточно большой путь развития (работы в этой области начались еще в 60-х годах прошлого века) и использовались при решении ряда сложных технических задач. Однако в настоящее время подобные системы эффективно и широко используются только в ECAD (Электронных CAD) при проектировании изделий электроники, микроэлектроники, радиоэлектроники и вычислительной техники. При проектировании же в других областях, и, в первую очередь, в MCAD (Машиностроительных CAD) системах, они используются эпизодически и не вносят существенного вклада в решение проектных задач. Более того, рынок систем моделирования динамических процессов в технических объектах в MCAD практически монополизирован фирмой MSC. Software Corporation и ее группой программных продуктов под общим названием MSC. ADAMS, которые широко используются в большинстве CAD/CAE/CAM систем. Программные модули MSC. ADAMS часто используются автономно и их основные потребители – ведущие производители автомобилей. В связи с этим большинство модулей MSC. ADAMS ориентированы на решение задач автомобильной промышленности.

Широкое использование систем моделирования динамических процессов в практике проектирования ставит перед разработчиками таких систем новые задачи:

-  комплексное моделирование разнородных технических объектов - электрических, механических, гидравлических, тепловых и т. д., а также возникновение задач, которые практически невозможно решить с помощью существующих программ;

-  надежность и высокая точность математических ядер таких программ;

-  обеспечение связи с другими подсистемами моделирования, в первую очередь с подсистемами геометрического моделирования;

-  разработка удобных средств для создания индивидуальных библиотек моделей разнообразных элементов технических систем с различной степенью точности и адекватности.

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

ОБОБЩЕННЫЕ ФОРМАЛЬНЫЕ СХЕМЫ

При математическом моделировании состояние любого технического объекта (ТО) в конкретный момент времени можно описать множеством некоторых физических величин, имеющих как постоянные значения, так и переменные значения в течение заданного времени моделирования. Следует отметить, что одна и та же физическая величина при моделировании может быть как постоянной, так и переменной. Например, при моделировании старта ракеты массу ракеты можно считать постоянной и в математической модели старта ракеты эта физическая величина будет параметром модели, наоборот, при моделировании полета ракеты масса ракеты будет изменяться и в математической модели полета ракеты та же самая физическая величина будет переменной модели. Обозначим множество переменных модели ТО вектором Vто(t). Конечной математической моделью динамических процессов в ТО в общем случае будет некоторая система дифференциально-алгебраических уравнений (ДАУ), замкнутая относительно этих переменных:

F(Vто, dVто/dt, t)=0 (1)

Вместе с тем любой ТО можно представить как множество элементов, связанных между собой множеством связей, поэтому его можно отобразить схемной моделью в виде некоторой обобщенной формальной схемы (ФС), состоящей из набора элементов, связанных между собой в некоторых узлах этой схемы. Такая схемная математическая модель технического объекта должна отображаться аналогичным множеством переменных Vфс(t) -> Vто(t), а конечной математической моделью такой ФС должна быть система ДАУ, аналогичная системе (1). Тогда можно разработать программную систему моделирования обобщенных ФС (ПА10), с помощью которой можно будет моделировать разнообразные технические объекты, описываемые системой ДАУ вида (1).

Элементный состав обобщенной ФС и типы связей между элементами должны давать возможность моделирования как любых дифференциально-алгебраических уравнений вида (1), так и разнообразных элементов технических объектов любой физической природы (электрических, магнитных, механических, гидравлических, тепловых, оптических, акустических и др.) и разнообразных связей между такими элементами. Предлагаемый подход является обобщением и развитием методов, предложенных в работах [1, 2, 7].

Состояние каждого узла связи в ФС в некоторый момент времени t опишем двумя типами безразмерных переменных – рис.1.

Рис.1

1) Переменная типа потока, втекающего в узел связи или вытекающего из узла связи Ik(t). Эти переменные удовлетворяют уравнению непрерывности потоков для каждого узла ФС, которое будет базовым для получения модели ФС на основе моделей отдельных элементов ФС:

I1+I2+…+Ik+… = 0 (2)

2) Переменная типа потенциала Pj(t), которая определяется по отношению к некоторому базовому узлу, потенциал которого принимается равным нулю. Разновидностью переменной этого типа будет переменная типа напряжения (разности потенциалов) Uji(t)=Pj(t)-Pi(t), которая определяется как разность потенциалов между узлами j и i.

Назовем узлы связи со стороны элементов ФС полюсами, тогда простейшими базовыми элементами ФС будут двухполюсники.

БАЗОВЫЕ ДВУХПОЛЮСНИКИ

Двухполюсники опишем двумя типами переменных, которые будут составлять исходный базис переменных ФС:

Iд – поток, протекающий через двухполюсник;

Uд – напряжение между полюсами двухполюсника.

Предположим, что ФС состоит только из n взаимосвязанных двухполюсников, тогда состояние ФС в любой момент времени в выбранном базисе переменных однозначно определяется множеством потоков и напряжений для всех двухполюсников и будет иметь размерность 2n:

Vфс(t) = (I1,U1,I2,U2,…In, Un).

Покажем, что для моделирования динамики любых технических систем, описываемых обыкновенными дифференциальными уравнениями вида (1), достаточно определить всего три типа базовых двухполюсных элементов, вместо 5, предложенных в [1], но при этом расширим возможности этих элементов:

1. Элемент типа зависимый источник потока (рис.2).:

Рис.2.

Элементы этого типа описываются двумя переменными – iI и uI. Математическая модель элементов данного типа представляет собой одно алгебраическое уравнение вида:

iI = f (Vфс, t), где

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

Обозначим через Jf – вектор частных производных базисной функции по переменным ФС - Vфэс:

Jf = ¶f/¶Vфс,

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

Вторая переменная в элементах этого типа не определяется:

uI = var.

2. Элемент типа емкость (рис.3).

Рис.3.

Элементы этого типа описываются двумя переменными – iC и uC. Математическая модель элементов данного типа представляет собой одно дифференциальное уравнение вида:

iC = CduC/dt, где C=f(Vфс, t).

Обозначим производную duC/dt через xpC (производная дифференцируемой переменной uC(t) по времени), тогда для целей использования современных неявных методов решения алгебро-дифференциальных уравнений ФС необходимо расширить базис выбранных переменных и включить эту переменную в вектор переменных ФС:

xpC=duC/dt ÎVфс.

В моделях элементов данного типа необходимо по тем же причинам, что и для вышеописанного элемента, вычислять элементы вектора частных производных:

¶iC/¶Vфс=xpCJf,

причем необходимо будет вычислять

¶iC/¶xpC=xpC¶C/¶xpC +C.

3. Элемент типа сопротивление (рис.4)

Рис.4

Элементы этого типа описываются двумя переменными – iR и uR. Математическая модель элементов данного типа представляет собой одно алгебраическое уравнение вида:

uR = R iR, где R=f(Vфс, t).

В моделях элементов данного типа необходимо по тем же причинам, что и для вышеописанных элементов, вычислять элементы вектора частных производных:

¶uR/¶Vфс=uR Jf, причем ¶uR/¶iR=iR ¶R/¶iR + R.

Прежде, чем рассматривать расширение состава базовых моделей ФС и расширение базиса переменных ФС покажем, что для любой системы ДАУ вида (1) можно всегда получить ФС, состоящую только из единичных емкостей, единичных сопротивлений и источников потока управляемых напряжениями на этих единичных элементах.

Систему ДАУ (1) всегда можно преобразовать к подсистемам:

dXp/dt = Fd (Xp,Ya, t) (3),

Ya = Fa(Xp,Ya, t) (4),

где размерность векторов Xp и Ya в общем случае может оказаться больше, чем размерность векторов X и Y. Систему (1) можно преобразовать к подсистемам (3) и (4), используя следующие простые эквивалентные математические преобразования:

1)  Выделение очередной переменной или ее производной в очередном уравнении системы (1) и перенос этой переменной в левую часть очередного преобразуемого уравнения.

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

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

4)  Одновременное прибавление и вычитание нужной переменной в очередном преобразуемом уравнении системы (1) и перенос этой переменной в левую часть преобразуемого уравнения.

Рассмотрим пример преобразования системы ДАУ 3-го порядка с тремя неизвестными:

d2x1/dt2 + (dx2/dt)2 + f1(y1,x1,x2) = 0

y13 + f2(y1,x2) =0

f3(dx1/dt, y1) = 0

Введем две дополнительные переменные – одну дифференцируемую x3 = dx1/dt и одну алгебраическую y2 = dx2/dt и преобразуем исходную систему ДАУ к подсистемам вида (3) и (4):

dx1/dt = x3

dx2/dt = y2

y2 = (-dx3/dt – y1 –f1(y1,x1,x2))/y2

y1 = -f2(y1,x2)/y12

dx3/dt = dx3/dt + f3(dx1/dt, y1)

Каждому i-му уравнению подсистемы (3) можно сопоставить ФС, состоящую из двух зависимых источников потока, единичной емкости и единичного сопротивления, показанных на рис. 5

Рис. 5

В этой подсхеме напряжение на единичной емкости UCdi соответствует переменной xdi, а напряжение на единичном сопротивлении URdi соответствует производной этой переменной dxdi/dt. Источник потока Idi соответствует правой части i-го уравнения.

Каждому j-му уравнению подсистемы (7) можно сопоставить ФС, состоящую из зависимого источника потока и единичного сопротивления, показанных на рис. 6

Рис. 6

В этой подсхеме напряжение на единичном сопротивлении URaj соответствует j-ой алгебраической переменной подсистемы (4). Источник потока Iai соответствует правой части j-го уравнения. Следовательно, любая система ДАУ вида (3) и (4) может быть представлена множеством однотипных подсхем, показанных на рис. 5 и рис. 6.

Таким образом, для любой системы ДАУ вида (1) можно всегда сопоставить ФС на основе введенных трех базовых двухполюсных элементов, поэтому данный набор базовых элементов можно считать необходимым и достаточным для моделирования любых технических объектов, динамика поведения которых описывается системами ДАУ вида (1). Однако для сохранения преемственности с предыдущими программами продолжим расширение базовых двухполюсных элементов и базиса переменных ФС.

4. Элемент типа проводимость (рис.7).

Рис.7.

Элементы этого типа описываются двумя переменными – iG и uG. Математическая модель элементов данного типа представляет собой одно алгебраическое уравнение вида:

iG = G uG, где G=f(Vфс, t).

В моделях элементов данного типа необходимо по тем же причинам, что и для вышеописанных элементов, вычислять элементы вектора частных производных:

¶iG/¶Vфс=uG Jf, причем ¶iG/¶uG=uG ¶G/¶uG +G.

5. Элемент типа источник напряжения (рис. 8).

Рис. 8.

Элементы этого типа описываются двумя переменными – iE и uE. Математическая модель элементов данного типа представляет собой одно алгебраическое уравнение вида:

uE = f(Vфс, t).

В моделях элементов данного типа необходимо по тем же причинам, что и для вышеописанных элементов, вычислять элементы вектора частных производных:

¶uE/¶Vфс=Jf.

Вторая переменная в элементах этого типа не определяется:

iE = var.

6. Элемент типа индуктивность (рис.9)

Рис.9.

Элементы этого типа описываются двумя переменными – iL и uL. Математическая модель элементов данного типа представляет собой одно дифференциальное уравнение вида:

uL = L diL/dt, где L=f(Vфс, t).

Обозначим производную diL/dt через xpL (производная дифференцируемой переменной ФС iL(t) по времени), тогда для целей использования современных неявных методов решения алгебро-дифференциальных уравнений ФС необходимо расширить базис переменных моделирования и включить эту переменную в вектор переменных ФС:

xpL=diL/dt ÎVфс.

В моделях элементов данного типа необходимо по тем же причинам, что и для вышеописанного элемента, необходимо вычислять элементы вектора частных производных:

¶uL/¶Vфс=xpL Jf, а также ¶uL/¶xpL=xpL ¶L/¶xpL +L.

Обозначим через m общее количество двухполюсников типа C и L, тогда вектор базисных переменных ФС, состоящей из n двухполюсников, будет иметь размерность 2n+m, поскольку в этот вектор добавлены производные xpC и xpL для каждой емкости и индуктивности. Обозначим через k общее количество узлов, связывающих между собой двухполюсники, за исключением базового узла, обозначим также через Ф(t) – вектор узловых потенциалов ФС, размерностью k и расширим базис переменных, включив этот вектор в вектор базисных переменных ФС: Ф(t)ÎVфс(t), который будет иметь в результате размерность 2n+m+k. Назовем этот базис обобщенным базисом ФС, в отличие от базиса переменных состояния, узловых потенциалов и расширенного, предложенных в работах [1,2].

Обозначим через фi и фj узлы подключения некоторого двухполюсника, тогда для этого двухполюсника можно добавить второе линейно-независимое уравнение, помимо основных уравнений (алгебраического для двухполюсников типа I, G,E, R или дифференциального для двухполюсников типа C и L), приведенных выше:

Uд = фi – фj.

Таким образом, для n двухполюсников мы будем иметь 2n линейно-независимых алгебро-дифференциальных уравнений. Для каждого узла ФС в свою очередь можно добавить линейно-независимое уравнение для суммы потоков, втекающих в этот узел:

I1+I2+…+Ik+… = 0 (смотри уравнение (2)).

В результате мы получим 2n+k линейно-независимых уравнений, замкнутых относительно базисных переменных ФС. m уравнений для m производных базисных переменных по времени дадут формулы интегрирования, в результате на каждом шаге интегрирования необходимо будет решать замкнутую систему нелинейных алгебраических уравнений (НАУ) размерностью 2n+m+k относительно неизвестного вектора Vфс. Это весьма высокая размерность, как будет видно из примеров, но система НАУ будет чрезвычайно разреженной (для линейных систем как правило 2-3 ненулевых элемента в строке соответствующей матрицы Якоби). Предложенный обобщенный базис снимает многие ограничения при разработке моделей элементов технических объектов, присущие базисам переменных состояния, узловых потенциалов и расширенному базису, которые так или иначе связаны с методами интегрирования ДАУ. Введение производных дифференцируемых переменных (xpC и xpL) в базис позволяет разрабатывать абсолютно независимо и автономно программы интегрирования систем ДАУ, основанные, как на явных, так и неявных формулах интегрирования.

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

ФУНКЦИОНАЛЬНЫЙ МНОГОПОЛЮСНИК

Функциональным многополюсником (ФМ) назовем элемент ФС, который представляет собой схему соединения двухполюсников и внутренних ФМ, в которой выделены внешние полюса Pi для связи с другими элементами ФС – рис.10.

Рис.10

Каждый i-ый полюс описывается двумя переменными:

1.  Поток iPi, втекающий в i–ый полюс, который назовем полюсным потоком.

2.  Потенциал i-го полюса uPi, который назовем полюсным потенциалом.

Обозначим через N общее количество внешних полюсов ФМ. Предположим, что ФМ состоит только из n двухполюсников, которые описываются вектором внутренних переменных ФМ – Vвнфм, размерностью 2n+k+m, где k – число внутренних узлов ФМ, а m – число элементов типа C и L внутри ФМ. Математической моделью такого ФМ будут внутренние уравнения для внутренних двухполюсников ФМ, размерностью 2n+k, а также N полюсных уравнений для суммы потоков, втекающих в каждый i-ый внешний полюс внутри ФМ. Таким образом, общее количество переменных ФМ будет на N больше, чем общее число уравнений, составляющих математическую модель ФМ. При включении ФМ в общую схему ФС недостающие N линейно-независимых уравнений будут получены для суммы потоков, втекающих в каждый i-ый внешний полюс снаружи ФМ.

В качестве математического ядра ПA10 предполагается использовать стандартную программу решения систем ДАУ общего вида – DMAN [4], которая была разработана, не привязываясь к выбранному базису переменных моделирования.

ПРОГРАММА ИНТЕГРИРОВАНИЯ СИСТЕМЫ ДАУ

Программа DMAN предназначена для интегрирования систем ДАУ общего вида, не разрешенных относительно производных:

F(XP,X,Y,t) = 0 (5)

где X - вектор дифференцируемых переменных системы (4) размерностью m; XP - вектор производных дифференцируемых переменных по времени размерностью m, т. е. XP = dX/dt; Y - вектор алгебраических переменных системы (4) размерностью k; t - независимая переменная (например, время); F - вектор-функция системы (4) размерностью n, где n = (m + k). Заданы начальные условия для вектора дифференцируемых переменных: X0 = X(t0), где t0 - начальный момент времени интегрирования. Начальные значения остальных переменных, т. е. XP0 и Y0 рассчитываются в программе DMAN перед началом интегрирования автоматически.

В качестве примера рассмотрим систему обыкновенных дифференциальных уравнений (ОДУ) второго порядка в нормальной форме Коши (уравнения Дуффинга):

dx1/dt = x2

dx2/dt = xx2 - x13 + sin(OMt)

где x1,x2 - дифференцируемые переменные системы ОДУ; OM - параметр.

Эта система может быть представлена как система ДАУ в общем виде по разному, например:

1-ый вариант:

f1(X) = xp1 - x2 = 0

f2(X) = xp2 - x1 + 0.2x2 + x13 - sin(OMt) = 0

где F = (f1,f2); X = (x1,x2); XP = (xp1,xp2); xp1 = dx1/dt;

xp2 = dx2/dt; m = 2; k = 0; n = 2.

2-ой вариант:

f1(X) = xp1 - x2 = 0

f2(X) = xp2 - x1 + 0.2x2 + y1 - sin(OMt) = 0

f3(X) = x13 - y1 = 0

где F = (f1,f2,f3); X = (x1,x2); XP = (xp1,xp2); xp1 = dx1/dt;

xp2 = dx2/dt; Y = y1; m = 2; k = 1; n = 3.

В процессе интегрирования система (5) на каждом шаге интегрирования дополняется системой:

G(XP,X,h) = 0, (6)

где h - шаг интегрирования; G - вектор-функция размерностью m, которая соответствует выбранному методу интегрирования. Например для неявного метода Эйлера G имеет вид:

G = XP - (1/h) * (X - Xn) = 0,

где Xn - значение вектора дифференцируемых переменных в момент времени tn; tn - момент времени в начале шага интегрирования h;

X - значение вектора дифференцируемых переменных в текущий момент времени t; h = (t - tn).

В программе DMAN реализованы 4 метода интегрирования систем ДАУ:

М1 - А-устойчивый неявный метод Эйлера первого порядка точности;

М2 - А-устойчивый неявный метод второго порядка точности;

М3 - новый А-устойчивый неявный метод четвертого порядка точности;

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

Во всех методах на каждом шаге интегрирования решается система нелинейных алгебраических уравнений относительно переменных X, XP, Y методом Ньютона, для которого необходимо вычислять матрицу Якоби для системы (5) и системы (6). Матрица Якоби для системы (6) вычисляется в программе DMAN автоматически. Матрица Якоби для системы (5) должна вычисляться в специальной подпрограмме пользователя FCT, наряду с вычислением вектор-функции F. Эта матрица состоит из двух подматриц:

RJ = ( RJ1, RJ2 )

где RJ - матрица Якоби размером (nn1); n1 = m + n. RJ1 – матрица частных производных вектор-функции F по переменным XP размером (nm), т. е. матрица RJ1 = [ dF/dXP ]. RJ2 - матрица частных производных вектор-функции F по переменным Z размером (nn), т. е. матрица RJ2 = [ dF/dZ ], где Z - вектор переменных системы ОДУ размерностью n, который включает дифференцируемые переменные и алгебраические переменные системы (5), т. е. Z = ( X,Y ).

Для приведенного выше примера матрица Якоби и соответствующие подматрицы имеют вид:

1-ый вариант:

RJ

1	0	0	-1
0	1	(-1+3x12)	0.2

RJ1

1	0
0	1

RJ2

0	-1
(-1+3x12)	0.2

2-ой вариант:

RJ

1	0	0	-1	0
0	1	-1	0.2	1
0	0	(3x12)	0	-1

RJ1

1	0
0	1
0	0

RJ2

0	-1	0
-1	0.2	1
(3x12)	0	-1

Начальным приближением для метода Ньютона являются значения переменных в начале шага интегрирования, поэтому, если на некотором шаге нет сходимости итераций, то текущий шаг h уменьшается вплоть до минимального HMIN, после чего выдается сообщение об ошибке (обычно из-за ошибок в формулах для определения элементов матрицы Якоби в подпрограмме пользователя или из-за расходимости самой системы ДАУ).

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

В конечном счете, вычисление любой вектор-функции F системы ДАУ можно свести к вычислению различных базисных функций одного аргумента типа f(x). Рассмотрим, например, функцию f(x) с ограниченной областью определения х, например, функция f(x)= определена только для x>=0. Поэтому при каждом обращении к такой функции необходимо проверить значение аргумента х и, если аргумент выходит из области определения функции, то следует установить параметр NBAD=1, не вычисляя значения функции. Некоторые базисные функции, теоретически определенные для любых значений аргументов, практически могут принимать значения, выходящие за пределы разрядной сетки ЭВМ, или накладываемыми реально возможными значениями аргумента и функции. Например, функция f(x)=ex при abs(x)>40 обычно выходит за допустимые значения в реальных физических системах. В этом случае также следует использовать параметр NBAD=1.

Важнейшей практической базисной функцией f(x) является кусочно-нелинейная функция, определяемая на множестве заданных точек. В точках излома эта функция имеет разрыв всех производных и для интегрирования таких функций без потери точности следует использовать параметр NBAD=2, фактически, начиная интегрирование после каждой точки излома с минимальным шагом интегрирования и с новыми начальными условиями. Таким образом, в ПА10 будет решена проблема интегрирования функций с разрывами производных без потери точности интегрирования.

ПРИМЕРЫ ПОЛУЧЕНИЯ МОДЕЛЕЙ В НОВОМ ОБОБЩЕННОМ БАЗИСЕ

Перейдем к рассмотрению примеров получения моделей ФС для реальных технических объектов - двух электрических схем, одной линейной и одной нелинейной. Сначала рассмотрим пример получения модели ФС для линейной электрической схемы, приведенной на рис. 11, в новом обобщенном базисе.

Рис.11

Вектор базисных переменных для этой ФС будет иметь размерность 12 и включает в себя следующие переменные:

Vфс = (uE, iE, uR, iR, uC, iC, uL, iL, ф1,ф2,xpC, xpL).

Система ДАУ для этой ФС будет иметь размерность n=10, m=2 и включает в себя следующие уравнения:

uE=f(t) или uE-f(t)=0

uE=ф1 или uE-ф1=0

uR=RiR или uR-RiR=0

uR=ф1-ф2 или uR-ф1+ф2=0

iC=CxpC или iC-CxpC=0

uC=ф2 или uC-ф2=0

uL=LxpL или uL=LxpL=0

uL=ф2 или uL-ф2=0

iE-iR=0

iR-iC-iL=0.

Отметим, что в базисе узловых потенциалов размерность системы ДАУ будет равна 2 (в пять раз меньше!), но при этом возможности моделирования отдельных элементов ФС будут ограничены (например, нельзя моделировать идеальные источники напряжения).

Для этой задачи вектор X=(uC, iL), вектор XP=(xpC, xpL), а вектор Y= (uE, iE, uR, iR, iC, uL, ф1,ф2). Вектор-функция F включает в себя 10 вышеприведенных уравнений. Полученная система ДАУ была решена с помощью программы DMAN. Полученное решение сравнивалось с аналитическим, при этом были получены численно более устойчивые решения при изменении коэффициентов этой системы ДАУ, по сравнению с системой ДАУ, полученной в базисе узловых потенциалов.

Рассмотрим пример получения модели ФМ для нелинейной модели диода.

На рис. 12 показаны эквивалентная схема диода из двухполюсных элементов и его представление как ФМ.

Рис.12.

Вектор переменных для ФМ диода будет иметь размерность 9 (из них первые 5 переменных являются внутренними) и включает в себя следующие переменные:

Vфмдиод = (uCд, iCд, xpCд, uIд, iIд, uк, iк, uа, iа)

Система ДАУ для ФМ диода (математическая модель диода в обобщенном базисе) будет иметь размерность 6 и включает в себя следующие уравнения (первые 4 уравнения – для 2 внутренних двухполюсников ФМ диода):

iCд-CдxpСд=0

uCд-uа+uк=0

iIд-(iтe(uIд/фт)+1)=0

uIд-uа+uк=0

iа-iCд-iIд=0

iк+iCд+iIд=0

В этих уравнениях iт, Cд, и фт – это параметры данной модели диода.

Такие уравнения будем включать в библиотеку моделей ПA10 и использовать при получении моделей сложных ФС, содержащих диоды.

Рассмотрим пример получения в обобщенном базисе с использованием ФМ диода модели ФС для нелинейной электрической схемы, приведенной на рис.13.

Рис.13

Вектор переменных Vфс включает в себя первую очередь переменные диода, для которых узловые потенциалы заменяются на потенциалы узлов подключения диода:

Vфмдиод = (uCд, iCд, xpCд, uIд, iIд, ф3, iк, ф2, iа)ÎVфэс

Затем в вектор переменных ФС добавляются остальные переменные, в результате вектор ФС будет иметь размерность 19 и включает в себя следующие переменные:

Vфс=( uCд, iCд, uIд, iIд, ф3, iк, ф2, iа uE, iE, uR1, iR1, uC, iC, uR2, iR2, ф1, xpCд, xpC).

Система ДАУ для этой ФС будет иметь размерность 17 и включает в себя следующие уравнения:

Сначала в эту систему добавляются 6 уравнений диода, в которых узловые потенциалы заменяются на узлы подключения диода:

iCд-CдxpСд=0

uCд-ф2+ф3=0

iIд-(iтe(uIд/фт)+1)=0

uIд-ф2+ф3=0

iа-iCд-iIд=0

iк+iCд+iIд=0

Затем добавляются 2 уравнения для суммы токов в узлах подключения диода:

iR1-iа=0

iк+iC1+iR2=0

Затем добавляются 8 уравнений для 4 двухполюсников:

uE-f(t)=0

uE-ф1=0

uR1-R1iR1=0

uR1-ф1+ф2=0

iC-CxpC=0

uC1-ф3=0

uR2-R2iR2=0

uR2-ф3=0

И в завершении формирования замкнутой системы ДАУ добавляется 1 уравнение для суммы токов в узле ф1:

iE-iR1=0

Для этой задачи вектор X=(uCд, uC), вектор XP=(xpCд, xpC), а вектор Y= (iCд, uIд, iIд, ф3, iк, ф2, iа uE, iE, uR1, iR1, iC, uR2, iR2, ф1). Вектор-функция F включает в себя 17 вышеприведенных уравнений. Полученная система ДАУ была решена с помощью программы решения систем ДАУ общего вида DMAN. Полученное решение сравнивалось с решением, полученным для метода узловых потенциалов в программе ПA7, при этом были получены численно более устойчивые решения при изменении коэффициентов этой системы ДАУ, по сравнению с системой ДАУ, полученной в базисе узловых потенциалов.

ЗАКЛЮЧЕНИЕ, ВЫВОДЫ

В результате проведенных исследований выполнен анализ и обоснование выбора нового обобщенного набора базовых элементов ФС и нового обобщенного базиса переменных ФС, в качестве основы для разработки новой, мобильной, открытой программы моделирования динамики технических систем - ПА10. Программа ПA10 будет базовым математическим ядром системы моделирования динамических процессов нового поколения и будет иметь следующие принципиальные особенности:

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

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

Дальнейшие исследования будут включать в себя:

-  разработку структуры всей системы моделирования динамики технических объектов и организации межмодульных связей в этой системе;

-  разработку новых методов решения систем линейных алгебраических уравнений с учетом очень сильной разреженности этих систем для обобщенного базиса при сохранении численной устойчивости вычислений;

-  разработку новой методики тестирования программы, основанной на всестороннем тестировании математического ядра на основе методов тестирования математических программ типа Mathlab, Mathematica и т. п

-  использование универсальных СУБД для хранения библиотек моделей элементов моделируемых технических объектов;

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

Результаты работы могут быть использованы при разработке модулей анализа динамических процессов для САПР разнообразных технических изделий.

ЛИТЕРАТУРА

1. Системы автоматизированного проектирования: Учебное пособие для вузов: В 9 кн./Под ред. .- М., Высшая школа, 1986.

2. , , -Радиотехника,1986, №9.

3. Маничев методы интегрирования обыкновенных дифференциальных уравнений для программ анализа радиоэлектронных схем. - Изв. вузов. Радиоэлектроника, том 32, N6, 1989.-С.45-49.

4. Маничев алгоритмы для программ анализа динамики технических систем. Вестник МГТУ, сер. Приборостроение, вып. 1, 1996. Стр. 48-56.

5. Устойчивость методов Рунге-Кутты для жестких нелинейных дифференциальных уравнений.- М.: Мир, 1988.-334 с.

6. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи: Пер. с англ.-М.:Мир, 199с.

7. Решение обыкновенных дифференциальных уравнений. Жесткие задачи: Пер. с англ.-М.:Мир, 200с.