авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 | 2 || 4 | 5 |   ...   | 11 |

«УДК 004 ББК 32.81 М34 Рекомендовано Редакционно-издательским советом ГрГУ им. Я. Купалы. Ред а к ц и он н а я ...»

-- [ Страница 3 ] --

Австралийская система прогнозирования пожаров [18] при реализации модели Ротермела в уравнении (1) вовсе использует константу B = 1.1. В экспериментах [20] по сжиганию слоя опа да сосновых иголок получено значение B = 0.8. Американская система FCCS [8] по умолчанию использует B = 1.2, но дает возможность эксперту ”поиграть” с этой константой при моделиро вании.

2. О теоретических моделях лесных пожаров Область применения эмпирических моделей ограничена, а точность прогноза не высока в ситуациях, отличных от тех, при которых были собраны опытные данные. Поэтому актуально построение и исследование теоретических моделей лесных пожаров, которые базируются на зако нах газовой динамики, тепломассопереноса, теории горения, других фундаментальных законах физики, химии. Такие модели записываются, как правило, в виде системы дифференциальных уравнений в частных производных. Верификация теоретических моделей весьма затруднитель на, однако только они, описывая развитие лесного пожара на основе общих законов и с учётом большого количества факторов, позволяют отвечать на широкий круг вопросов [1, 21].

Как было отмечено в [1] относительно теоретических моделей распространения лесных пожа ров, чаще других используется в качестве основы и считается наиболее обоснованной модель А.М.

Гришина [21] – модель многофазной пористой реагирующей среды, которая описывает процес сы нагрева, сушки, пиролиза и горения древесины. Эта модель базируется на фундаментальных законах физики (сохранения массы, количества движения и энергии), теоретически обоснова на, частично подтверждена натурными и лабораторными экспериментами. Именно эта модель была принята авторами за основу при создании компьютерного комплекса по моделированию распространения вершинных верховых пожаров. В работах [13, 14, 15] приведены результаты детализации модели, уточнены определяющие функции и параметры применительно к харак терным условиям возникновения и протекания процессов распространения лесных пожаров в Беларуси, изложены основы алгоритмической и компьютерной реализации. Эффективность мо дели демонстрируется в основном на примере расчета беглых верховых лесных пожаров, как наиболее быстро распространяющихся и наносящих значительный материальный ущерб.

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

3. Адаптированная математическая модель низовых лесных пожаров В приведенной ниже математической модели распространения лесных пожаров вычисляется эволюция распределений следующих величин: T – осредненная по высоте слоя ЛГМ температура в Кельвинах;

j, j = 1, 2, 3, 4 – объемные доли многофазной реагирующей среды, где 1 соответ ствует сухому органическому веществу лесных горючих материалов, 2 – связанной с ЛГМ воде в жидко-капельном состоянии, 3 – коксику (конденсированному продукту пиролиза), 4 – ми неральной части ЛГМ (золе);

c, = 1, 2, 3 – массовые концентрации компонентов газовой фазы, где c1 соответствует кислороду, c2 – горючим газам (горючим компонентам продуктов пиролиза), c3 – смеси остальных газов (инертных компонентов воздуха, водяного пара, инертных продуктов реакций пиролиза, горения коксика и окисления горючих газов).

При построении модели учитываются следующие физико-химические процессы: ”подвод” теплоты в результате конвекции, теплопроводности и излучения, прогрев ЛГМ, их сушка и пиро лиз, горение газообразных и твердых продуктов пиролиза. Вывод уравнений, обоснование модели, численная схема и особенности организации расчетов приведены в публикациях [13 - 15, 22] для беглых верховых лесных пожаров. В данной работе делается ”формальный” переход от верховых пожаров к низовым, а именно, используется та же система дифференциальных уравнений, но с новыми значениями физико-химических величин, характерными для низовых лесных пожаров.

Относительно неизвестных функций j (j = 1, 2, 3, 4), c ( = 1, 2, 3), T, зависящих от вре мени и пространственных координат, формулируется начально-краевая задача.

3.1. Уравнения модели:

1 2 3 = 1 (1, T ), = 2 (2, T ), = 3 (1, 3, c1, c2, T ), = 0, (2) t t t t c1 + (V, gradc1 ) div(5 DT gradc1 ) = c1 (1, 2, 3, c1, c2, T ), (3) t c2 + (V, gradc2 ) div(5 DT gradc2 ) = c2 (1, 2, 3, c1, c2, T ), (4) t 5 cp5 (V, gradT ) div(T gradT ) T + = T (1, 2, 3, 4, c1, c2, T ). (5) t 5 cp5 + j j cpj j= 3.2. Определяющие функции:

R1 R2 c R1 Mc R 1 (1, T ) = 2 (2, T ) =,, 3 (1, 3, c1, c2, T ) =, (6) 1 2 3 M 1 1 R51 c1 Q (c1 c1 ), c1 (1, 2, 3, c1, c2, T ) = (7) 5 cp5 h 1 R52 c2 Q (c2 c2 ), c2 (1, 2, 3, c1, c2, T ) = (8) 5 cp5 h T ) 4R T q5 R5 q2 R2 + q3 R3 h (T T (1, 2, 3, 4, c1, c2, T ) =, (9) 5 cp5 + j j cpj j= 3 T c MC Q = (1 c )R1 + R2 + c = 1, 5 =, R3, (10) M T M M =1 = E1 E2 E R2 = k02 T 1/2 2 2 exp( R1 = k01 1 1 exp( ), ), R3 = k03 s 3 5 c1 exp( ), (11) RT RT RT R5 M R51 = R3 R52 = (1 c ) R1 R5,, (12) 2M M2 ECO R5 = 5 min(c2, c1 )kCO exp( ). (13) 2M1 RT Здесь t – время;

V – вектор равновесной скорости ветра;

T – невозмущенная температура окру жающей среды в Кельвинах;

j, j = 1, 2, 3, 4 – истинные плотности j-й фазы;

5 – плотность газовой фазы (смеси газов);

– невозмущенная плотность смеси газов (плотность воздуха);

c1 и c2 – массовые концентрации кислорода и горючих газов в невозмущенной атмосфере;

M, = 1, 2, 3 – молекулярные массы компонентов газовой фазы;

MC – молекулярная масса углерода, M – молекулярная масса (невозмущенная) воздуха;

cpj, j = 1, 2, 3, 4 – теплоемкости j-й фазы;

cp5 – теплоемкость газовой фазы;

Q – массовая скорость образования газовой фазы;

T – коэффициент турбулентной теплопроводности;

DT – коэффициент диффузии;

q2, q3 и q – тепловые эффекты процессов испарения, горения конденсированного горючего и газообразно го горючего продукта пиролиза соответственно;

h – высота слоя ЛГМ (h = h2 h1, где h и h1 – высоты верхней и нижней границ слоя соответственно);

– коэффициент теплообмена, между окружающим воздухом и моделируемым слоем ЛГМ;

R – интегральный коэффициент поглощения;

– постоянная Стефана-Больцмана;

R1 – массовая скорость реакции пиролиза (раз ложения под действием высокой температуры с выделением горючих газов) сухого органического вещества ЛГМ;

R2 – массовая скорость испарения воды из ЛГМ (сушки ЛГМ);

R3 – массовая скорость реакции горения коксового остатка;

R51, R52 – массовые скорости образования кислоро да, горючих газов;

R5 – массовая скорость реакции горения (окисления) горючих газов;

k01, k02, k03 – предэкспоненты химических реакций, E1, E2, E3 – энергии активации химических реакций, R – универсальная газовая постоянная, s – удельная поверхность конденсированного продук та пиролиза (коксика), C – коксовое число ЛГМ, – доля газообразных горючих продуктов пиролиза ЛГМ.

Записанная выше система (2) - (5) – нелинейная система уравнений в частных производных, она включает теплофизические коэффициенты, скорости сушки и реакций пиролиза и горения, целый ряд других эмпирических постоянных. Совокупность коэффициентов уравнений пред ставляет оснастку принятой математической модели, является ее существенной и неотъемлемой частью. Качественная оснастка модели позволяет адекватно воспроизводить в вычислительных экспериментах наблюдаемые реальные процессы распространения лесных пожаров. Наполнение модели осуществлялось с учётом приведенных в литературных источниках сведений, причём предпочтение отдавалось тем, где упоминались экспериментальные исследования или отмеча лось, что они получены с помощью решения обратных задач механики реагирующих сред. В настоящей работе также несколько величин определялись в вычислительных экспериментах на этапе адаптации модели.

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

3.3. О начальных и граничных условиях. Начальные распределения объёмных долей, массовых концентраций и температуры задаются во всей области определения решения. Обо значим область моделирования через G, а ее границу –. Область G можно разбить на три подобласти G = Gf ire G G+. Эти подобласти не обязательно односвязные, их взаимная гео метрия может быть довольно сложной, например, для многоочаговых пожаров. Через G и G+ обозначим подобласти сгоревших и не затронутых пожаром участков области G удаленные от очагов горения на достаточное расстояние (эти подобласти характеризуются невозмущенными значениями c1, c2 и T ).

В зоне горения Gf ire начальные распределения величин j (j = 1, 2, 3, 4), c ( = 1, 2, 3) и T должны быть ”самосогласованными”, так как между всеми этими величинами есть определённые связи физического характера, которые необходимо учитывать. Это довольно сложные вопросы, требующие отдельного рассмотрения.

Начальные распределения в зонах ”выжженная” G и ”невыжженная” G+ задаются следу ющим образом:

T |t=0 (G G+ ) = T ;

(14) c1 |t=0 (G c2 |t=0 (G c3 |t=0 (G G+ ) = 1 c1 c2 ;

G+ ) = c1, G+ ) = c2, (15) 0 1 |t=0 (G+ ) = 2 |t=0 (G+ ) = (1 )W 4 0;

,, 3 (G+ )|t=0 = 0, (16) 1 1 |t=0 (G ) = 0, 2 |t=0 (G ) = 0, 3 |t=0 (G ) = C. (17) Здесь предполагается, что в G ЛГМ полностью сгорели;

0 – плотность типичного слоя ЛГМ, плотности j (j = 1, 2, 3, 4), W – влагосодержание слоя ЛГМ и – зольность лесных го рючих могут быть неоднородными, соответствующие распределения в рассматриваемой области моделирования считаются заданными функциями координат.

На границе области моделирования G удобно использовать так называемые мягкие гра ничные условия:

T c1 c = 0, = 0, = 0. (18) n n n 4. О корректности приведенной постановки Система дифференциальных уравнений принятой теоретической модели является нелиней ной, коэффициенты и правые части уравнений зависят от всех искомых переменных. Физико химические процессы имеют ”взрывной” характер протекания, описывающие их функции выра жаются экспоненциальными зависимостями от температуры. Температура в зоне горения (обыч но относительно узкой) меняется от невозмущенного до значения более тысячи градусов Кельви на. Также в зоне пожара по направлению его распространения распределения массовых концен траций горючих газов и окислителя, объемных долей компонент ЛГМ имеют большие градиенты.

Записанная выше система уравнений даже в одномерном случае очень сложна;

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

В каком ракурсе авторами изучались вопросы корректности приведенной постановки зада чи? Анализировались численные решения и их поведение при возмущении определяющих функ ций и начальных условий. Уравнения аппроксимировались по явным однородным схемам, опи сание дано в [13, 14].

Следует заметить, что для низового пожара, по аналогии с верховым [14], возможны несколь ко режимов развития чрезвычайной ситуации. В частности, режимы затухания и устойчивого распространения. Например, для режима устойчивого распространения характерны устанавле ние определенной геометрии и формы профилей температуры, физико-химических параметров.

Они без существенных изменений переносятся во всех направлениях развития процесса (по, про тив ветра и во флангах).

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

Рисунок 1 – Распределения температуры и объемной доли сухих ЛГМ (время t = 0, 1, 2 и 5 секунд) На рисунке 1 для четырех моментов времени показаны графики, иллюстрирующие эволюцию рас пределений температуры и объемной доли сухих ЛГМ в элементе распространяющейся по ветру кромки пожара.

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

Рисунок 2 – Распределения температуры и объемной доли сухих ЛГМ. Второй комплект формы возмущений начального профиля температуры Рисунок 3 – Распределения температуры и объемной доли сухих ЛГМ (время t = 0 и 5 секунд).

Сдвиг распределения температуры относительно объемной доли сухих ЛГМ Графики рисунка 2 иллюстрируют другую конфигурацию возмущений начальных распределе ний температуры. Но, как и в предыдущих вариантах, с течением времени все распределения принимают одинаковые фиксированные формы и переносятся по потоку с постоянной скоро стью. Однако в отличие от графиков на рисунке 1 для ”пунктирных” профилей наблюдается сдвиг решений по пространству на фиксированное расстояние по сравнению со ”сплошными” графиками. Иллюстрации на рисунке 3 показывают, что на сдвиг оказывает влияние взаимное расположение распределений температуры и объемной доли сухих ЛГМ в начальный момент.

5. Об инструментарии вычислительных экспериментов Изучение особенностей численных решений на этапах исследования их точности, анализа поведения базовых и асимптотики возмущённых решений, обнаружения качественных особен ностей моделируемых процессов предполагают многовариантные расчёты. Пришлось провести сотни продолжительных по времени вычислительных экспериментов. Моделирование пожаров авторами проводилось на ”обычных” компьютерах. Расчет большинства вариантов в двумерной постановке по времени занимает до нескольких суток. Для реализации возможностей прове дения многовариантных расчетов с разными параметрами сетки и характеристиками модели, без повтора признанных удовлетворительными фрагментов расчётов, разработаны специальные модули автоматизации проведения и протоколирования серий вычислительных экспериментов.

Запрограммированный инструментарий позволяет изменять параметры моделирования, сетки и сохранять любые наборы результатов в базе данных [22]. Конечно же ввиду вычислительной сложности задачи очевидным становится необходимость разработки и исследования методов распараллеливания вычислительного алгоритма задачи, проведение соответствующих расчетов на суперкомпьютере.

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

Приведем представительные результаты численного моделирования распространения пожа ра для случая равномерно распределенного сильно высушенного слоя опада и мхов, когда пло щадь возгорания – квадрат. При анализе получаемых численных решений основное внимание уделялось контролю распределений вблизи центральной линии по направлению ветра, где реше ния локально одномерны, а также скоростей распространения фронта и тыла. В расчётах исполь зовались следующие значения: V = 2 м/с, W = 10%, = 0, 0 = 10 кг/м3, h = 0.1 м, T = С, 1 = 500 кг/м3, 2 = 1000 кг/м3, 3 = 200 кг/м3, 4 = 200 кг/м3, = 1.15 кг/м3, c1 = 0.23, c2 = 0, M1 = 32, M2 = 28, M3 = 29, MC = 12, M = 29, cp1 = 2000 Дж/(кг·К), cp2 = Дж/(кг·К), cp3 = 900 Дж/(кг·К), cp4 = 1000 Дж/(кг·К), cp5 = 1000 Дж/(кг·К), q2 = 3 · 106 Дж/кг, q3 = 1.2 · 107 Дж/кг, q5 = 107 Дж/кг, T = 1000 Дж/(м·с·К), DT = 1 м2 /с, = 300 Вт/(м2 ·K), c = 0.1, = 0.8, R = 0.3, s = 1000 м1, E1 /R = 9400 К, E2 /R = 6000 К, E3 /R = 10000 К, k01 = 3.63 · 104 с1, k02 = 6 · 105 К0.5 с1, k03 = 1000м·с1.

Графики рисунка 4 иллюстрируют рассчитанные распределения температуры T и объемной доли сухих ЛГМ 1. Для приведенных моментов времени на начальной стадии форма профилей меняется – этап развития процесса. Далее профили температуры и объемной доли сухих ЛГМ устанавливаются, в последующем их формы меняются незначительно, возмущения переносятся по и против ветра с разными скоростями – режим устойчивого распространения. Для рассмат риваемого варианта профили устанавливаются после 15 с.

При величине скорости ветра V = 2 м/с рассчитанная скорость фронта пожара по направ лению ветра равна f ront = 0.24 м/с, против ветра (задняя кромка): back = 0.19 м/с.

На рисунке 5 приведена полученная по результатам расчётов серии вариантов зависимость скорости распространения фронта и тыла кромки пожара от равновесной скорости ветра V при остальных фиксированных переменных для случая режима устойчивого распространения. Из графика видно, что зависимость скорости распространения фронта пожара от скорости ветра близка к линейной;

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

Обобщая приведенные и другие уже полученные, но не обсуждаемые в настоящей работе результаты расчётов, следует отметить, что, используя разработанный программный комплекс, методику вычислительного эксперимента и интерпретации результатов, можно составить таб лицы зависимостей скорости распространения фронта и тыла кромки пожара от равновесной скорости ветра V для типичных по каждому региону наборов значений основных характеристик территорий. К таковым относятся коэффициенты влагосодержания, плотности сухого орга нического вещества лесных горючих материалов, связанной с ЛГМ воды в жидко-капельном состоянии, коксика, минеральной части ЛГМ. Реализация сказанного позволит на основе исполь зования табулированных данных для каждого конкретного пожароопасного региона значительно расширить возможности методики Ротермела, повысив точность прогнозов распространения пожаров.

Рисунок 4 – Графики распределения температуры и объемной доли сухих ЛГМ Рисунок 5 – Зависимости скоростей фронта и тыла кромки пожара от скорости ветра Список литературы 1. Mathematical models and calculation systems for the study of wildland re behaviour / E. Pastor [et al] // Progress in Energy and Combustion Science. – 2003. – Vol. 29. – P. 139–153.

2. McArthur, A. G. Weather and grassland re behaviour / A. G. McArthur. – Commonwealth of Australia:

Forestry Research Institute, Leaet No. 100, 1966. – 23 p.

3. Rothermel, R. C. A Mathematical model for predicting re spread in wildland fuels / R. C. Rothermel. – Ogden: USDA Forest Service, Research Paper INT–115, 1972. – 43 p.

4. Баровик, Д. В. Состояние проблемы и результаты компьютерного прогнозирования распространения лесных пожаров / Д. В. Баровик, В. Б. Таранчук // Вестник БГУ. Серия 1. – 2011. – № 3. – С. 78–84.

5. Баровик, Д. В. О развитии методики Ротермела и реализации двумерной компьютерной модели про гноза распространения лесных пожаров / Д. В. Баровик, В. Б. Таранчук // Веснiк Вiцебскага дзяр жаўнага ўнiверсiтэта. – 2011. – № 6 (66). – С. 5–11.

6. Heinsch, F. A. BehavePlus re modeling system, version 5.0: Design and Features. General Technical Report RMRS-GTR-249 / F. A. Heinsch, P. L. Andrews. – Fort Collins, CO: USDA Forest Service, Rocky Mountain Research Station, 2010. – 111 p.

7. Finney, M. A. FARSITE: Fire Area Simulator – model development and evaluation / M. A. Finney. – Ogden:

USDA Forest Sevice, Research Paper RMRS–RP–4, 2004. – 47 p.

8. An overview of the Fuel Characteristic Classication System – Qualifying, classifying, and creating fuelbeds for resource planning / R. D. Ottmar [et al.] // Can. J. For. Res. – 2007. – № 37. – Р. 2383–2393.

9. US Forest Service. The Fuel Characteristic Classication System. [Электронный ресурс]. – Режим доступа:

http://www.fs.fed.us/pnw/fera/fccs/index.shtml 10. Scott, J. H. Standard re behaviour fuel models: a comprehansive set for use with Rothermel’s surface re spred model / J. H. Scott, R. E. Burgan // USDA Forest Service, General Technical Report RMRS–GTR– 153, 2005. – 72 p.

11. Баровик, Д. В. Методические и алгоритмические основы программного комплекса ”Расчет и визуали зация динамики лесного пожара” / Д. В. Баровик, В. И. Корзюк, В. Б. Таранчук // Чрезвычайные ситуации: предупреждение и ликвидация. – 2011. –№ 2 (30). – С. 22–33.

12. Баровик, Д. В., Таранчук В. Б. Адаптация модели Ротермела для реализации в программном ком плексе прогноза распространения лесных пожаров / Д. В. Баровик, В. Б. Таранчук // Научный интернет журнал ”Технологии техносферной безопасности”. – 2011. – № 6 (40). – Режим доступа:

http://ipb.mos.ru/ttb/2011-6/2011-6.html 13. Barovik, D. V. Mathematical modelling of running crown forest res / D. V. Barovik, V. B. Taranchuk // Mathematical Modelling and Analysis. – 2010. – Vol. 15, N 2. – P. 161–174.

14. Баровик, Д. В. Об особенностях адаптации математических моде-лей вершинных верховых лесных по жаров / Д. В. Баровик, В. Б. Таранчук // Вестник БГУ. Серия 1, Физика, Математика, Информатика.

– 2010. – N 1. – С. 138–143.

15. Баровик, Д. В. Численная реализация математической модели вер-ховых лесных пожаров / Д. В. Ба ровик, В. Б. Таранчук // Весцi БДПУ. Серия 3, Физика, Математика, Информатика. – 2010. – N 2. – С. 40–44.

16. Reformulation of Rothermel’s wildland re behavior model for heterogeneous fuelbeds / D.V.Sandenberg [et al.] // Can. J. For. Res. – 2007. – № 37. – P. 2438–2455.

17. Nelson, R. M. Flame characteristics of wind-driven surface res / R. M. Nelson, C. W. Adkins // Can. J.

For. Res. – 1986. – № 16. – P. 1293–1300.

18. McCaw, W. L. Predicting re spread in Western Australia mallee–heath shrubland. PhD Thesis / W. L. McCaw;

University of New South Wales. – Canberra, Australia, 1998. – 256 p.

19. Rate of spread of free–burning res in woody fuels in a wind tunnel / W. R. Catchpole [et al] // Combustion Science and Technology. – 1998. – Vol. 131. – P. 1–37.

20. Pagni, J. Flame spread through porous fuels / J. Pagni, G. Peterson // Fourteenth Symposium (International) on Combustion, USDA Forest Service, Washington, DC. – Pittsburgh: The Combustion Institute, 1973. – P. 1099–1107.

21. Гришин, А. М. Математическое моделирование лесных пожаров и новые способы борьбы с ними / А. М. Гришин. – Новосибирск: Наука, 1992. – 408 с.

22. Баровик, Д. В. Базы данных результатов численного моделирования (на примере задачи распростра нения лесных пожаров) / Д. В. Баровик // Вестник БГУ. Серия 1, Физика, Математика, Информатика.

– 2010. – N 2. – С. 170–174.

Баровик Дмитрий Валентинович, ассистент кафедры кафедры компьютерных технологий и систем БГУ, кандидат физико-математических наук, BarovikD@gmail.com.

Корзюк Виктор Иванович, заведующий кафедрой математической физики БГУ, доктор физико математических наук, профессор, член-корреспондент НАН Беларуси, заведующий отделом математиче ской физики Института математики НАН Беларуси, Korzyuk@bsu.by.

Таранчук Валерий Борисович, заведующий кафедрой компьютерных технологий и систем БГУ, док тор физико-математических наук, профессор, Taranchuk@bsu.by.

УДК 621.372. П. В. БОЙКАЧЕВ, Г. А. ФИЛИППОВИЧ, В. Ф. БЕЛЕВИЧ МОДЕЛИРОВАНИЕ ПАРАМЕТРОВ СВЧ ТРАНЗИСТОРОВ В ШИРОКОМ ДИАПАЗОНЕ ЧАСТОТ Предложены модели для входного и выходного сопротивлений сверхвысокочастотного транзистора, составленные по результатам измерения S-параметров. Модели характери зуют свойства транзистора в широком диапазоне частот.

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

- расчет параметров на основе физической эквивалентной схемы транзистора;

- измерение параметров транзистора на дискретном ряде частот.

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

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

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

Подобный подход использован специалистами Санкт-Петербургского электротехнического института для моделирования сопротивлений антенн [1]. Однако их подход использует эквива лентные схемы, состоящие не более чем из трех элементов, для расчетов которых используются результаты измерений на двух частотах. Такой подход не обеспечит необходимую точность для моделирования параметров СВЧ транзисторов. Поэтому нами была разработана модель СВЧ транзистора, включающая до 6 элементов физической эквивалентной схемы.

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

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

Входное сопротивление транзистора Для иллюстрации подхода возьмем S-параметры современного малошумящего транзисто ра ATF-34143, измеренные в диапазоне частот 0.5... 18 ГГц при параметрах питания Vds = V, Ids = 20 mA. При определении входных и выходных сопротивлений транзистор полагается однонаправленным [2], а сопротивления рассчитываются по измеренным S-параметрам. Резуль таты расчета активной Re Z(f ) и реактивной Im Z(f ) составляющих входного сопротивления транзистора Z(f ) представлены на рис. 1.

Рисунок 1 – Функции входного сопротивления по результатам измерений Из анализа полученных зависимостей нетрудно установить наличие последовательного и параллельного резонансов во входной цепи транзистора и небольшого сопротивления в ее по следовательной ветви. Такому характеру частотных зависимостей Re Z(f ) и Im Z(f ) можно поставить в соответствие эквивалентную схему, изображенную на рис. 2.

Рисунок 2 – Эквивалентная схема входного сопротивления Действительная и мнимая части сопротивления схемы на рис. 2 определяются следующими выражениями:

2 RL2 Re Z(j) = r + (1) R2 (1 2 L2 C2 )2 + (L2 ) (1 2 L1 C1 )2 2 R2 L2 (1 2 L2 C2 ) Im Z(j) = +2 (2) R (1 2 L2 C2 )2 + (L2 ) C Определим параметры элементов эквивалентной схемы. Значения элементов схемы r и R можно установить непосредственно из частотной зависимости Re Z(f ) на рис.1: r = min Re Z(f ), R = max Re Z(f ) r. Для расчета остальных параметров потребуются значения и в четырех точках частотного диапазона. С этой целью выберем частоты последовательного и параллельного резонансов 01 и 02 соответственно. Частота 1 будет использоваться для расчета L1 и C1, поэтому ее целесообразно выбрать в той области частотного диапазона, где элементы L1 и C создают наибольший вклад в Im Z(f ), т.е. в области резонансной частоты 01. Это позволит повысить точность расчета L1 и C1. Из этих же соображений выбирается частота 2 в области резонанса 02.

Используя значения Re Z(f ) и Im Z(f ) на выбранных частотах, находим значения элемен тов:

2 1 R(2 )R ;

R(2 ) = Re Z(f2 ) r L2 = ;

C2 = 2 (3) 2 (R R( )) 2 02 C R2 L2 (1 2 L2 C2 ) [ImZ(f1 ) X(1 )] 1 L1 = ;

C1 = 2L ;

X(1 ) = (4) 2 2 R2 (1 2 L2 C2 )2 + (L2 ) 1 02 01 Полученные выражения (1) – (4) можно использовать для расчета значений элементов экви валентной схемы на рис. 2. Используя значения действительной и мнимой составляющих входного сопротивления транзистора, приведенные на рис. 1, находим: r = 8 Ом, R = 270 Ом, 1 = 2.4 пф, L1 = 0.658 нГ, L2 = 0.698 нГ, C2 = 0.27 пф.

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

Рисунок 3 – Функции входного сопротивления модели (кривые 1,3) в сравнении с функциями входного сопротивления построенной по результатам измерений (кривые 2,4).

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

Выходное сопротивление транзистора Результаты расчета активной Re Z(f ) и реактивной Im Z(f ) составляющих выходного сопро тивления транзистора Z(f ) представлены на рис.4. Поведение зависимостей Re Z(f ) и Im Z(f ) указывает на наличие слабо выраженного параллельного резонанса, емкостного характера со противления в последовательной ветви и заметного резистивного сопротивления этой ветви. Эк вивалентная схема с такими свойствами представлена на рис.5.

Рисунок 4 – Функции выходного сопротивления по результатам измерений.

Рисунок 5 – Эквивалентная схема выходного сопротивления.

2 RL2 r Re Z(j) = r1 + +2 (5) 1 + 2 R2 2 C1 2 R (1 2 L2 C2 )2 + (L2 ) R2 2 C1 2 R2 L2 (1 2 L2 C2 ) Im Z(j) = +2 (6) 1 + 2 R2 2 C1 2 R (1 2 L2 C2 )2 + (L2 ) Расчет значений элементов r2, R, L2 и C2 аналогичен расчету элементов предыдущей схемы с учетом того, что элементу r2 соответствует элемент r. Выбор частоты производится в обла сти частот, где существенно влияние элементов r2 и C1, т.е. в низкочастотной части диапазона.

Совместное решение выражений (5) и (6) дает следующий результат:

(X(1 ) + Im (Z(1 ))2 + (Re Z(1 ) r1 R(1 )) r2 = (7) Re Z(1 ) r1 R(1 )) X(1 ) + Im (Z(1 ) C1 = (8) 1 [(X(1 ) + Im (Z(1 ))2 + (Re Z(1 ) r1 R(1 ))2 ] где R( + 1) и X(1 ) – последние составляющие выражений (5) и (6) соответственно. Они определяются на частоте для рассчитанных ранее R, L2 и C2.

Результаты расчета элементов эквивалентной схемы: r1 = 30 Ом, r2 = 99.5 Ом, R = Ом, 1 = 3.2 пф, L2 = 0.756 нГ, C2 = 0.25 пф. Частотные зависимости действительной и мнимой частей сопротивления эквивалентной схемы в сравнении с функцией выходного сопротивления по результатам измерений приведены на рис.6.

Рисунок 6 – Функции выходного сопротивления модели (кривая 1,3) в сравнении с функцией выходного сопротивления по результатам измерений (кривая 2,4).

Сопоставление зависимостей на рис. 6 свидетельствует о хорошем приближении эквива лентной схемы.

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

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

Список литературы 1. Проектирование радиопередающих устройств с применением ЭВМ/ Под ред. О. В. Алексеева // М.:

Радио и связь, 1987. – 392 с.

2. Шварц, Н. З. Линейные транзисторные усилители СВЧ / Н. З. Шварц // М.: Радио и связь, 1980. – 368 с.

Бабаев Петр Александрович, профессор кафедры прикладной электродинамики Воронежского го сударственного университета, доктор медицинских наук, профессор, babaev.pa@mail.ru.

Дешко Михаил Сергеевич, аспирант кафедры прикладной электродинамики Воронежского государ ственного университета, deshko@yahoo.com.

Ледков Егор Владимирович, студент 5 курса факультета математики и информатики Воронежского государственного университета, egoretz@uk.yahoo.com.

УДК 616.125.4-008.314-073. В. К. БОЙКО, А. М. КАДАН ПРОГРАММНАЯ РЕАЛИЗАЦИЯ АЛГОРИТМА РЕШЕНИЯ ЗАДАЧИ О МИНИМАЛЬНОМ ЧИСЛЕ ВХОДОВ ЛИНЕЙНЫХ СТАЦИОНАРНЫХ СИСТЕМ Одним из алгоритмов решения задачи о минимальном числе входов является проверка условий управляемости на специально построенном множестве матриц управления. В работе приводится программная реализация такого алгоритма, оцениваются его вычис лительные возможности и обсуждаются возможности его применения для различных классов систем Введение В ряде практических задач для уже заданного объекта нужно подобрать оператор управле ния, при котором объект будет управляем, причем, иногда требуется, чтобы выбранный оператор был в некотором смысле самым простым. Одним из наиболее распространенных требований такого рода, является выбор оператора управления, обеспечивающего свойство управляемости при минимальном числе управляющих воздействий, то есть при минимальной размерности вектора управления. Задача в такой постановке называется задачей о минимальном числе входов. Впервые она была поставлена и решена в работах [1, 2] для линейных стационарных дифференциальных систем, где было показано что минимальное число входов такой системы равно числу нетривиальных инвариантных многочленов матрицы системы. Там же указан способ построения матрицы управления, при которой система управляема таким числом входов. Для стационарных систем с запаздыванием аналогичный подход был использован В.М.Марченко.

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

Постановка задачи Рассмотрим систему dx(t) = Ax(t) + Bu(t), t [0, t1 ], t1 0 (1) dt Будем считать, что матрица A системы (1) задана, а матрицу B разрешается выбирать так, чтобы полученная система (1) стала управляемой, причем число столбцов матрицы B должно быть минимальным. Эту задачу будем называть задачей о минимальном числе входов (управля ющих воздействий) системы.

Как известно, система (1) управляема при заданной матрице управления, тогда и только тогда, когда rank B, AB, A2 B,..., An1 B = n (2) Это значит, что среди миноров матрицы условия (2) есть отличный от нуля минор порядка n.

Элементы матрицы условия (2) являются линейными формами относительно элементов матрицы, поэтому любой минор n-го порядка будет представлять собой однородную форму степени n от n r переменных bpq матрицы.

Приведем один вспомогательный результат. Пусть n (y1,..., ym ) – однородная форма сте пени n от m переменных y1,..., ym m j j jm aj1 j2 jm y11 y22 ym n (y1,..., ym ) = (3) i=1ji =n,ji Известно, что число коэффициентов этой формы равно числу размещений с повторениями из m элементов по n (n + m 1)!

n Cm = (m 1)!n!

Выберем систему значений переменных y1l... yml, l = 1,..., Cm так, чтобы определитель n m m j yili ;

ji = n;

ji 0;

det (4) n i=1 i=1;

l=1,...,Cm был отличен от нуля.

Справедлива Лемма. Если для системы значений переменных y1l,..., yml, l = 1,..., Cm справедливы ра n венства (y1k,..., yml ) = 0 (5) то n (y1,..., ym ) 0, т.е. форма (3) есть тождественный нуль.

Справедливость леммы следует из того, что если для формы (3) равенство (5) рассматривать как линейную однородную систему относительно коэффициентов aj1 j2...jm, то условие (4) говорит, что система (5) имеет только тривиальное решение, т.е. все коэффициенты aj1 j2...jm = 0.

Укажем способ построения системы значений y1l,..., yml удовлетворяющих условию (4).

Пусть mi n yil = ln, = 0, = ±1, i = 1,..., m, i = 1,..., Cm (6) Тогда m m1 +j m2 +...+j ) j yil = l(j1 n 2n m, i= причем различным наборам чисел j1, j2,..., m таких, что m ji = n, ji 0, i = 1,m i= соответствуют различные числа j1 nm1 + j2 nm2 +... + jm, как числа в системе счисления с основанием n. Определитель (4) в этом случае является опре делителем Вандермонда и отличен от нуля.

Вернемся к задаче вычисления минимального числа входов. Запишем элементы матрицы в виде строки b11 b21... bn1 b12... bn2... b1r... bnr и, как в случае формы (3), построим систему значений вида (6) nr(p1)rq bl = ln, = 0, = ±1, (7) pq которые образуют семейство {B(n, r)} матриц управления размерности n r, для которых вы полняется условие аналогичное (4). Если на семействе {B(n, r)} критерий (2) не выполняется, то это означает, что любой минор n-го порядка, как однородная форма степени n от n · r пере менных bpq, тождественно равен нулю, и система (1) не управляема ни при какой n r матрице управления.

Таким образом, справедлива Теорема. Минимальное число входов системы (1), равно наименьшему r, при котором в семействе {B(n, r)} существует матрица, на которой выполняется критерий (2).

Построив семейства {B(n, 1)}, {B(n, 2)},..., {B(n, r)},..., {B(n, n 1)} и, проверяя на них критерий (2), мы решаем задачу о минимальном числе входов.

Для матрицы управления с минимальным числом r входов выполняется условие rankB = r поэтому, удалив из семейств {B(n, 1)}, {B(n, 2)},..., {B(n, r)},..., {B(n, n 1)} матрицы, неудо влетворяющие этому условию, мы упростим проверку критерия (2).

Показано, что условие (4) можно заменить рядом независимых условий меньшей размерно сти. Этот результат и положен в основу программируемого алгоритма.

Алгоритм построения семейства матриц управления {B(n, r)} 1. Задаем размерность системы n.

2. Задаем число входов r, 1 r n 1.

3. Вводим запись n r матрицы B в виде строки b11 b21... bn1 b12 b22... bn2... b1r b2r... bnr (8).

4. Из однородной формы степени n с n r неизвестными строки (8) выбираем произведения, образованные одинаковыми элементами строки (8), и строим определитель, удовлетворяющий условию (4). Матрицу такого определителя будем называть клеткой.

1. Учитывая, что матрица управления B, соответствующая минимальному числу входов, должна иметь полный ранг rankB = r, начинаем построение с клетки, образованной = r элементами строки (8). Обозначим набор этих элементов = (1... ).

2. Сосотавляем всевозможные разбиения числа n на положительных слагаемых. Записываем эти разбиения в виде i = (1,..., ), где i – номер разбиения. Число таких разбиений равно i i (n 1)!

( 1)!(n )!

Это и есть размер клетки. Наборы запоминаем в массиве Набор показателей степеней.

3. Строим строку клетки.

1. Из чисел 1, n выбираем набор из чисел, значения которых присваиваем 1.... Число таких наборов равно (n 1)!

n!

Cn = !(n )! ( 1)!(n )!

Наборы запоминаем в массиве Набор оОснований степеней.

2. Для каждого набора степеней (1,..., ) строим произведение i i i i i = 1 1...

где i = (1... ), i – набор из массива Набор показателей степеней.

i i 3. Записываем эти произведения в строку (n1)!

1 2 (1)!(n)!

...

Формирование одной строки завершено.

Значения 1... записываем в массив Окончательный набор оснований степеней 4. Выбираем из массива Набор оснований степеней следующую запись и присваиваем ее значения переменным 1...

Повторяем шаги 4.3.ii, 4.3.iii.

В результате построена следующая строка клетки. Назовем ее строка k+1.

Проверяем ее на независимость от предыдущих строк согласно условию строка строка rank =k+1 (9)...

строка k+ Если (9) выполнено, то значение переменных сохраняется в массив Окончательный набор оснований степеней и переходим к построению новой строки.

Если (9) не выполнено, то значение переменных не сохраняется и переходим к построению новой строки.

(n1)!

5. Процесс останавливается, когда сформировано (1)!(n)! строк, удовлетворяющих условию (9).

Если массив Набор оснований степеней исчерпан, а остановка не произошла, формируем сооб щение На базе набора оснований степеней формироание клетки не возможно.

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

6. Из строки b11... bnr строим всевозможные сочетания по элементов. Их количество равно (nr)!

Cnr =.

()!(nr )!

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

1. Из массива Вспомогательный набор матриц для -блока берем запись Bl = (bl 1 j1,..., bl j ), l = 1, Cnr i i s и присваиваем значения B l = s, где s - s-я запись в Набор оснований степеней. s = (n1)!

1, (1)!(n)!

2. Для каждого s строим матрицу в которой выбранные элементы определим в 5.1, а остальные Bls элементы равны нулю.

3. Проверим условие rankBls = r (10) Если (3) не выполнено, берем следующее Bls.

4. Первую матрицу, для которой выполнено (10), записываем в массив Матрицы управления для -блока. Обозначим ее B 5. Следующую матрицу Bls, для которой выполнено (10), проверяем на условие rank(Bt, Bls ) r, t = 1, p (11) где p - количество матриц Bt в массиве Матрицы управления для -блока в момент проверки.

Если для всех t = 1, p условие (11) выполняется, то Bls добавляется в массив Матрицы управ ления для -блока.

Если условие (11) не выполнено при каком-то t, то дальнейшая проверка прекращается и пере ходим к следующей записи Bls из Вспомогательный набор матриц для -блока.

В результате получим массив Матрицы управления для -блока.

6. В 4.1 полагаем = r + 1,..., n и повторяем всё до 5.5. Если условие (4) выполнено для всех t = 1, p, то далее проверим условие rank(Bqt, Bls ) r q = r, 1 : (12) по всем матрицам каждого массива Матрицы управления для q-блока Если условие (12) где-то не выполнилось, то переходим к следующей записи Bls из массива Вспо могательный массив матриц для -блока.

Таким образом получаем массивы Матрицы управления для r-блока, Матрицы управления для r + 1-блока,..., Матрицы управления для n-блока.

7. Соединяем все полученные массивы и получаем массив n r матрицы управления Особенности реализации алгоритма Для реализации вычислений согласно предложенной схеме был разработан компекс про грамм на языке Python.

Выбор языка Python обусловлен BSD-подобной лицензией на свободное ПО [4], совмести мой с GNU General Public License (GPL). Это позволяет в дальнейшем избежать ограничений, связанных с коммерческим использование разработанного комплекса программ.

При расчетах использована математическая библиотека Python – numpy, раздел линейной алгебры.

Комплекс включает модули - для подготовки (генерации и отбора) матриц управления;

- для тестирования построенных массивов n r матриц управления.

В настоящее время получены мартрицы размерности n r, 3 n 8, 2 r n.

Информация о количестве полученных матриц управления и времени, затраченном на их получение, приведена в таблицах 1-4. Вычисления выполнялись с использованием процессора i5, RAM 4GB.

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

Первоначально был использован метод отбора с вычислительной сложностью O(n2 ), где n – количество сгенеритрованных матриц. Метод отбора предполагал последовательное папарное сравнение всех матриц A1 и A2 размерности n r и проверки, имеет ли образованная из них матрица (A1 |A2 ) размерности n 2r ранг, равный r. В случае положительного ответа матрица A2 исключалась из даленейшего рассмотрения. Время выполнения операции отбора оказалось весьма большим.

Для сокращения времени выполнения операции отбора все сгенерированные матрицы пред вартельно были разбиты на классы эквивалентности по следующему критерию. Матрица A отно силась к классу Ki1 i2...ip, если все строки матрицы A с номерами i1 i2... ip содержали ненулевые элементы.

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

Таблица 1 – Количество сгенирированных матриц n/r 2 3 4 5 6 7 6 11392 58030 73920 19800 720 – – 7 74004 628761 1609580 2150800 125160 99567 – 8 477332 2386660 7159980 16467954 19761544 2823077 Таблица 2 – Время генерации матриц (сек.) n/r 2 3 4 5 6 7 6 2 13 46 125 191 – – 7 11 122 624 21547 30256 44287 – 8 113 1099 6674 20022 37516 52673 Таблица 3 – Количество матриц после очистки n/r 2 3 4 5 6 7 6 6557 6188 734 27 1 – – 7 49070 32658 26754 1533 776 102 – 8 349547 389012 307654 210633 178654 109876 Таблица 4 – Время на очистку (сек.) n/r 2 3 4 5 6 7 6 10 21 15 4 1 – – 7 201 1674 2487 509 197 89 – 8 2934 4098 7893 9046 6300 3897 Web-сервис для расчета минимального числа входов Данный алгоритм построения семейства матриц управления {B(n, r)} положен в основу раз рабатываемого авторами web-сервиса, который позволял бы научным работникам, инженерам, а также студентам использовать его для решения задачи нахождения минимального числа входов.

Интерактивные сервисы подобного типа в настоящее время достаточно популярны и активно используются как компоненты современных образовательных платформ [5, 6], но для решения не только учебных задач.

Компонентами Web-сервиса являются:

База матриц управления {B(n, r)} для различных n, r;

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

Программный модуль для расчета минимального числа входов;

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


Список литературы 1. Vogt, W. G. The minimum number of inputs required for the complete controllability of linear stationary dynamical systems / W. G. Vogt, G. A. Gullen / IEEE Trans. automat. Contr. – 1967, 12, № 3.

2. Габасов, Р. К теории управляемости динамических систем / Р. Габасов // Дифференцю уравнения. – 1968. – Т. 4, № 9.

3. Марченко, В. М. К задаче вычисления минимального числа входов / В. М. Марченко // Известия АН БССР. Сер. физ.-мат. наук. – 1974, № 5.

4. History and License [Electronic resource] / Python Programming Language – Ocial Website.

Documentation. – Mode of access: http://docs.python.org/2/license.html. – Date of access: 05.04.2013.

5. About edX [Electronic resource] / The Future of Online Educationhttps – Mode of access:

http://www.edx.org/. – Date of access: 05.04.2013.

6. Class2Go [Electronic resource] / Stanford Online Classes. – Mode of access: http://class.stanford.edu/. – Date of access: 05.04.2013.

Бойко Валерий Константинович, заведующий кафедрой высшей математики Гродненского государ ственного университета им.Янки Купалы, доцент, кандидат физико-математических наук, boiko@grsu.by.

Кадан Александр Михайлович, заведующий кафедрой системного программирования и компью терной безопасности Гродненского государственного университета им.Янки Купалы, доцент, кандидат технических наук, kadan@mf.grsu.by.

УДК 539.3+(616.314-089.23) С. М. БОСЯКОВ, А. В. ВИНОКУРОВА, А. Н. ДОСТА КОНЕЧНО-ЭЛЕМЕНТНЫЙ АНАЛИЗ ПЕРЕМЕЩЕНИЙ КОСТЕЙ ЧЕРЕПА ЧЕЛОВЕКА ПОСЛЕ АКТИВАЦИИ ОРТОДОНТИЧЕСКОГО АППАРАТА ДЛЯ РАСШИРЕНИЯ ВЕРХНЕЙ ЧЕЛЮСТИ В настоящей работе представлены результаты конечно-элементного моделирования верх нечелюстного расширения черепа человека. Определены суммарные перемещения костей верхнечелюстного комплекса, а также в зубов верхней челюсти при различной высоте расположения пластинок аппарата по отношению к верхнему небу. Моделирование черепа и зубов верхней челюсти выполнено на основании данных спиральной томографии сухого черепа взрослого человека. Показано, что при уменьшении расстояния между пластинками и верхним небом перемещения точек черепа существенно изменяются как по величине, так и по направлению. Полученные результаты могут быть использованы при проектировании конструкций ортодонтических аппаратов Hyrax с учетом индивидуальных особенностей пациентов.

Введение Для лечения сужения верхней челюсти у взрослых пациентов с врожденными расщели нами губы и неба используют комплексный подход, который предусматривает хирургическую подготовку, ортодонтическое лечение, а также протетические мероприятия. Ортодонтическое лечение осуществляется с применением различных ортодонтических аппаратов, развивающих значительные усилия, которые могут приводить к разрыву срединного небного и других швов черепа, а также к смещению и повороту зубов верхнего зубного ряда [1, 2]. Клинические наблюдения показывают, что конструктивные особенности аппаратов для верхнечелюстного расширения влияют на интенсивность и характер перемещения зубов и челюстных костей. В связи с этим представляется актуальным и практически востребованным исследование влияния конструкции ортодонтического аппарата на напряженно-деформированное состояние черепа человека. Настоящая работа посвящена конечно-элементному моделированию перемещений костей верхнечелюстного комплекса после активации ортодонтического аппарата Hyrax.

Разработка модели черепа и ортодонтического аппарата Стереолитографическая модель черепа получена с применением программы для обработки медицинских изображений MIMICS 14.12 (Materialise’s Interactive Medical Image Control Systems, Materialise BV, Leuven, Belgium) на основании 210 томогрофических изображений высушенно го трупного интактного черепа взрослого человека с хорошо сохранившимися альвеолярными отростками и зубами. Шаг томографических срезов составляет 1 мм.

При генерации стереолитографической (STL) модели удалены первые о вторые премоляры верхней челюсти (14, 15 и 24, 25 зубы), а также постоянные моляры (16 и 26 зубы), на которые устанавливается ортодонтический аппарат. Конечно-элементная модель получена после обработ ки STL-модели в модуле 3-matic 6.1 пакета MIMICS. Импортирование дискретной модели черепа в программную среду пакета ANSYS Worbench 13 (ANSYS Inc., USA) выполнено с помощью его компонента Finite Element Modeler. Полученная модель содержит 26 445 узлов и 91 731 эле мент типа Solid72. Максимальный размер элемента равен 1 мм, конечно-элементное разбиение – свободное. Результат импортирования и исходная STL-модель показаны на рисунке 1.

Рисунок 1 – Стереолитографическая (A) и конечно-элементная (B) модели интактного черепа человека Твердотельные модели премоляров и первого моляра черепа человека также получены на основании томографических данных черепа человека с применением CAD-пакета SolidWorks (SolidWorks Corporation, USA). С помощью графических примитивов этого пакета построена мо дель ортодонтического аппарата Hyrax с коронками, устанавливаемыми на премоляры (14 и зубы) и постоянные моляры (16 и 26 зубы), а также стержнями, соединяющими коронки и дей ствующими на премоляры (15 и 25 зубы). Отметим, что в ортодонтическом аппарате Hyrax, используемом для раскрытия небного шва, активная часть (винт) фиксируется к пассивной (ко ронке или кольцу) посредством проволочных элементов, припаянных к небной поверхности колец или коронок, укрепляемых цементом на моляры и премоляры [3]. Реальная конструкция орто донтического аппарата и его модель в пакете SolidWorks представлены на рисунке 2.

Рисунок 2 – Ортодонтический аппарат (A) и его геометрическая модель (B), установленная на четвертые премоляры и постоянные моляры Длина и ширина пластинок модели ортодонтического аппарата, представленного на рисунке 2 (B), составляют 10 мм и 4 мм соответственно;

радиус поперечного сечения стержней равен мм, длина стержней составляет 27 мм, толщина коронок 0,2 мм.

Генерация конечно-элементной модели ортодонтического аппарата, премоляров и первого моляра выполнена с применением компонента Mechanical Model пакета Ansys Worbench 13. Об щее количество элементов и узлов составляет 13320 и 26375 соответственно;

максимальный раз мер полученных элементов для аппарата равен 1 мм, тип элемента – Mesh200. Модель аппарата и зубов после импортирования в модуль Finite Element Modeler добавлена в конечно-элементную модель черепа. При этом контакт между коронками ортодонтического аппарата и зубами верх ней челюсти задан с помощью контактных элементов CONTA174 и TARGE170. Контакт между черепом и зубами задан с помощью контактных элементов CONTA173 и TARGE170. Конечно элементная модель черепа с установленным ортодонтическим аппаратом показана на рисунке 3.

Рисунок 3 – Дискретная модель черепа человека с установленным на премоляры и моляры ортодонтическим аппаратом: A – вид спереди;

B – вид снизу Граничные условия, накладываемые на череп, соответствовали жесткой заделке узлов, находящихся в окрестности большого затылочного отверстия [2, 4, 5]. Перемещение каждой пластинки составляет 0,4 мм (соответствует активации винта ортодонтического аппарата на полоборота). Перемещения задавались только в горизонтальном направлении (в направлении оси 0x).

Определение перемещений и эквивалентных напряжений Конечно-элементный расчет напряженно-деформированного состояния черепа с установлен ным ортодонтическим аппаратом проводился для случая, когда модуль упругости материала, из которого изготовлены пластинка и стержни аппарата, составлял E = 200 ГПа, коэффици ент Пуассона = 0, 3. Модуль упругости компактной (кортикальной) костной ткани и зубов равен E = 15 ГПа и E = 20 ГПа соответственно, коэффициент Пуассона как для кортикальной костной ткани, так и для зубов = 0, 3 [6]. В результате получены диаграммы распределения эквивалентных напряжений и суммарных перемещений для черепа, опорных зубов и ортодонти ческого аппарата. На рисунке 4 приведены диаграммы суммарных перемещений для точек черепа в случаях A и B установки ортодонтического аппарата. В случае A стержни и пластинки орто донтического аппарата расположены в одной плоскости (в плоскости x0y), в случае B пластинки ортодонтического аппарата расположены на 8 мм выше (ближе к верхнему небу) по отношению к горизонтальному положению.

Из рисунка 4 видно, что направление перемещений существенно изменяется при установке пластинок ортодонтического аппарата у верхнего неба (рисунок 4, B) по сравнению с направле ниями перемещений при горизонтальном положении стержней и пластинок аппарата (рисунок 4, A). В случае A вектора суммарных перемещений верхней челюсти расположены над плоскостью x0y. В случае B вектора суммарных перемещений направлены практически вдоль отрицатель ного направления оси 0z.

В случае A наибольшие суммарные перемещения наблюдаются для точек верхней челюсти и составляют 0,2212 мм. Суммарные перемещения других костей черепа существенно меньше, в частности, наибольшие перемещения скуловой кости 0, 074 мм;

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


Рисунок 4 – Суммарные перемещения черепа (значения перемещений приведены в миллиметрах) после активации ортодонтического аппарата: A – стержни и пластинки отодонтического аппарата расположены в плоскости x0y;

B – пластинки ортодонтического аппарата смещены к верхнему небу на 8 мм Существенное отличие максимальных суммарных перемещений в случаях A и B согласует ся с распределением напряжений в костях черепа. Соответствующие диаграммы эквивалентных напряжений (напряжения по фон Мизесу) для двух рассматриваемых случаев закрепления ор тодонтического аппарата показаны на рисунке 5.

Рисунок 5 – Распределение эквивалентных напряжений в лицевой части черепа после активации ортодонтического аппарата (значения напряжений приведены в МПа): A – стержни и пластинки отодонтического аппарата расположены в плоскости x0y;

B – пластинки ортодонтического аппарата смещены к верхнему небу на 8 мм Из рисунка 5 видно, что в случае A значительные напряжения возникают только в области верхней челюсти;

максимальные эквивалентные напряжения составляют 127,15 МПа. В случае B напряжения возникают в области альвеолярных и лобных отростков верхней челюсти, а так же в области лобно-носового шва;

наибольшие эквивалентные напряжения равны 35, 5 МПа.

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

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

Список литературы 1. Chasonas, S. J. Observation of orthopedic force distribution produced by maxillary orthodontic appliances / S. J. Chaconas, A. A. Caputo // Am. J. Orthod. – 1982. – Vol. 82, № 16. – P. 492–501.

2. Provatidis, C. On the FEM modeling of craniofacial changes during rapid maxillary expansion / C. Provatidis, B. Georgiopoulos, A. Kotinas, J. P. McDonald // Med. Eng. Phys. – 2007. – Vol. 29. – P. 566–579.

3. Аболмасов, Н. Г. Ортодонтия / Н. Г. Аболмасов, Н. Н. Аболмасов. – М.: МЕДпресс-информ, 2008. – 424 с.

4. Gautam, P. Biomechanical response of the maxillofacial skeleton to transpalatal orthopedic force in a unilateral palatal cleft / P. Gautam, L. Zhao, P. Patel // Angl. Orthod. – 2011. – Vol. 81, № 3. – P. 503–509.

5. Jafari, A. Study of stress distribution and displacement of various craniofacial structures following application of transverse orthopedic forces A three-dimensional FEM study / A. Jafari, K. S. Shetty, M. Kumar// Angl. Orthod. – 2003. – Vol. 73, № 1. – P. 12–20.

6. Чуйко, А. Н. Особенности биомеханики в стоматологии / А. Н. Чуйко, В. Е. Вовк. – Харьков: Прапор, 2006. – 304 с.

Босяков Сергей Михайлоович, доцент кафедры теоретической и прикладной механики Белорусского государственного университета, кандидат физико-математичских наук, доцент, bosiakov@bsu.by.

Винокурова Анастасия Владимировна, магистрант кафедры теоретической и прикладной механики Белорусского государственного университета, janeraven@mail.ru.

Доста Андрей Николавевич, доцент кафедры ортопедической стоматологии Белорусского государ ственного медицинского университета, кандидат медицинских наук, доцент, dostastom75@mail.ru.

УДК 519. Ю. М. ВУВУНИКЯН ВЫЧИСЛЕНИЕ СПЕКТРАЛЬНЫХ ХАРАКТЕРИСТИК ПОЛИНОМИАЛЬНЫХ ЭВОЛЮЦИОННЫХ ОПЕРАТОРОВ, ПОРОЖДЁННЫХ ДИФФЕРЕНЦИАЛЬНЫМИ УРАВНЕНИЯМИ Методы математического моделирования являются основными средствами исследования сложных нелинейных динамических систем, для описания которых в последнее время используются полиномиальные эволюционные операторы с обобщенными спектральными характеристиками. При этом нелинейные и динамические свойства системы полностью характеризуются последовательностью многомерных импульсных характеристик. В качестве нелинейных систем в данной работе рассматриваются системы, определяемые нелинейными дифференциальными уравнениями первого порядка. Для аналитического и численного исследования систем используются комплексные коэффициенты передачи в виде спектральных характеристик. Заметим, что для описания систем, определяемых дифференциальными уравнениями, необходимо в качестве характеристик рассматривать обобщенные функции.

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

Для достижения этой цели поставлены следующие задачи:

* получение замкнутых формул для обобщенных импульсных характеристик для конкрет ных полиномиальных эволюционных операторов;

* получение формул для обобщенных спектральных характеристик для квазиобратных опе раторов к данным эволюционным операторам;

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

* построение и оптимизация алгоритма, позволяющего вычислять спектральные характери стики квазиобратного оператора.

Импульсные характеристики эволюционных операторов Определение. Полилинейным эволюционным оператором степени m называется оператор m A, определяемый следующим образом: Ax = Sn (an xn ), где Sn – оператор сокращения n= переменных порядка n, действующий по формуле:

Sn : ((t1, t2,..., tn ) f (t1, t2,..., tn ) (t f (t, t,..., t))), an -– обобщенная функция, называемая импульсной характеристикой оператора A, xn – тензор ная степень, – n-мерная свертка.

В данной работе будем рассматривать полилинейные эволюционные операторы, порождён ные дифференциальными уравнениями первого порядка следующего вида:

1. x + ax2 = f (t).

Уравнение 1) определяет полиномиальный эволюционный оператор A второй степени со следующими импульсными характеристиками: a1 =, a2 = a = a 2, где – дельта-функция, первая обобщенная производная от дельта-функции.

2. x + ax2 + bx = f (t).

Уравнение 2) определяет полиномиальный эволюционный оператор A второй степени со следующими импульсными характеристиками: a1 = + b, a2 = a = a 2.

3. x + ax3 + bx2 + cx + d = f (t).

Уравнение 3) определяет полиномиальный эволюционный оператор A третьй степени со сле дующими импульсными характеристиками: a1 = + c, a2 = b = b 2, a3 = a = a 3.

4. ax + bx + cxn = f (t).

Уравнение 4) определяет полиномиальный эволюционный оператор A степени m со следую щими импульсными характеристиками: a1 = a + b, am = c m, an = 0 n = m.

Обобщённое преобразование Лапласа an импульсной характеристики an называется спек тральной характеристикой порядка n эволюционного оператора A, а семейство (n ) - системой a = 1, и то, что обоб спектральных характеристик эволюционного оператора A. Учитывая, что щенное преобразование Лапласа тензорного произведения равно тензорному произведению обоб щенных преобразований Лапласа, для эволюционных операторов, порожденных уравнениями 1)-4), получим систему спектральных характеристик.

1. Для эволюционного оператора A, порожденным уравнением 1), имеем:

a1 (1 ) = 1, a2 (1, 2 ) = a.

2. Для эволюционного оператора A, порожденным уравнением 2), имеем:

a1 (1 ) = 1 + b, a2 (1, 2 ) = a.

3. Для эволюционного оператора A, порожденным уравнением 3), имеем:

a1 (1 ) = 1 + c, a2 (1, 2 ) = b, a3 (1, 2, 3 ) = a.

4. Для эволюционного оператора A, порожденным уравнением 4), имеем:

a1 (1 ) = a1 + b, am (1, 2,..., m ) = c, an = 0 n = m.

Квазиобратные операторы к рассматриваемым эволюционным операторам Определение. Пусть A - полиномиальный эволюционный оператор степени m:

m Sn (an xn ), (x X), Ax = n= B - полиномиальный оператор Вольтера-Винера степени r:

r Sk (bk xk ), (x X).

Bx = k= И пусть C – оператор Вольтера-Винера, являющийся композицией операторов B и A, т.е. C = BA, а F – оператор Вольтера-Винера, являющийся композицией операторов B и A, т.е. F = AB.

Оператор B называется левым квазиобратным степени r к оператору A, если C1 = I и Cn = при 2 n r.

Оператор B называется правым квазиобратным степени r к оператору A, если F1 = I и Fn = 0 при 2 n r.

Определение. Оператор B называется квазиобратным степени r к оператору A, если он од новременно является левым и правым квазиобратным степени r к оператору A.

В соответствии с теоремой композиции имеем:

n m (1 + 2 +... + n, n +1 +... + n +n, n +...+n cn () = b m1 +1 +... + nm ) 1 1 1 2 m=1 n1 +n2 +...+nm =n n1 (1, 2,...n1 )n2 (n1 +1,...n1 +n2 )...nm (n1 +n2, n1 +...+nm1 +1,...nm ).

a a a Оператор B является квазиобратным степени r к оператору A, если выполняются условия:

( c ), c1 () = т.е.

1 ()1 () = 1,, b a и для любого n = 2;

3;

...;

r выполнено равенство n m (1 + 2 +... + n, n +1 +... + n +n, n +...+n b m1 +1 +... + nm ) 1 1 1 2 m=1 n1 +n2 +...+nm =n n1 (1, 2,...n1 )n2 (n1 +1,...n1 +n2 )...nm (n1 +n2, n1 +...+nm1 +1,...nm ) = 0.

a a a Симметричные спектральные характеристики квазиобратных операторов Рассмотрим теперь случай симметричности спектральных характеристик квазиобратных операторов. Отметим, что если спектральные характеристики не симметричны, то нужно приме нить к каждой спектральной характеристике оператор симметризации sym:

symf = f, m!

Gm где суммирование ведется по группе Gm всех перестановок степени m. При этом будем предпола гать, что спектральные характеристики a1, a2,..., an исходного эволюционного оператора A уже симметричны.

Обозначим симметричные характеристики через s, s,..., s. Получаем:

b2 b3 br s = b2 (1, 2 ) + b2 (2, 1 ) = 1 · a2 (1, 2 ) a2 (2, 1 ) b2 + 2 2 a1 (1 )1 (2 )1 (1 + 2 ) a1 (1 )1 (2 )1 (1 + 2 ) a a a a a2 (1, 2 ) a1 (1 )1 (2 )1 (1 + 2 ).

a a Аналогично получим формулы для s,..., s.

b3 br Алгоритм построения разбиений Для реализации формулы вычисления спектральных характеристик на языке высокого уров ня необходимо составить алгоритм, состоящий из простых действий, выполняемых функциями и процедурами. Это требует разложение формулы на более простые составляющие. Реализация суммы n1 +n2 +...+nm =n требует построения разбиений числа n на m слагаемых. При построении характеристик исполь зовались функции для нахождения разбиений и для нахождения количества разбиений. Нахож дение количества разбиений не составляет особенных трудностей. Однако, построение разбиений – уже не тривиальная задача. Проанализировав данную задачу, можно сделать следующие выво ды: Первое разбиение строится следующим образом: n-m+1,11,...,1m-1. Максимальный элемент в разбиениях равен n-m+1. Первый элемент в разбиениях принимает все значения от 1 до n-m+1.

Сгруппируем разбиения по первому элементу и обозначим его как i. Тогда можно заметить, что каждая группа за вычетом первого элемента представляет собой разбиения числа n-i на m- слагаемое. Это наблюдение позволяет построить разбиения рекурсивно со следующим услови ем выхода из рекурсии: при m=1 возвращать входящее n. Данный алгоритм позволяет строить разбиения для любых n и m, n=m. По такому же принципу можно строить и обычные разби ения. В этом случае дополнительно осуществляется выбор максимального элемента разбиения как min(n,m), а условием выхода из рекурсии будет n=1. Данный принцип построения разбиений не производит перебор в лексико-графическом порядке, что позволяет повысить производитель ность.

Программный модуль на основе данного алгоритма позволяет значительно ускорить вре мя нахождения спектральных характеристик, независим других коммерческих продуктов. Про грамма работает с использованием фреймворка.NET 4.0 в современных операционных системах семейства Microsoft Windows.

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

Список литературы 1. Вувуникян, Ю. М. Эволюционные операторы с обобщёнными импульсными и спектральными харак теристиками: монография/ Ю. М. Вувуникян. – Гродно: ГрГу, 2007. – 224 с.

2. Рейнгольд, Э. Комбинаторные алгоритмы. Теория и практика/ Э. Рейнгольд, Ю. Нивергельт, Н. Део.

– Москва: Мир, 1980. – 476 с.

3. Эндрюс, Г. Теория разбиений/ Г. Эндрюс. – М.: Наука, 1982. – 256 с.

Вувуникян Юрий Микиртычевич, профессор кафедры теории функций, функционального анализа и прикладной математики Гродненского государственного университета им. Янки Купалы, доктор физико математических наук, профессор, vuv@grsu.by.

УДК 621. С. В. ГАРКУША, Д. В. АНДРУШКО МАТЕМАТИЧЕСКАЯ МОДЕЛЬ РАСПРЕДЕЛЕНИЯ ЧАСТОТНО-ВРЕМЕННОГО РЕСУРСА В НИСХОДЯЩЕМ КАНАЛЕ СВЯЗИ ТЕХНОЛОГИИ WIMAX Предлагается математическая модель распределения частотно-временного ресурса в нисходящем канале связи стандарта IEEE 802.16. Предложенная модель направлена на формирование одного пакета данных нисходящего канала для каждой пользовательской станции, что позволяет минимизировать количество служебных сообщений передаваемых по используемому частотному каналу связи.

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

Исходя из технологических особенностей технологии WiMAX, задача распределения частотного и временного ресурса должна быть сформулирована как задача распределения слотов между пользовательскими станциями (ПС) сети и соответствующего объединения их в пакеты данных в зависимости от заявленной скорости передачи и параметров QoS.

Анализ известных решений Анализ известных решений [1–5] показал, что повышение производительности технологии WiMAX и обеспечение QoS может быть обеспечено путем как раздельного, так и согласованного решения задач распределения частотного и временного ресурса. Так варианты решений зада чи распределения частотного ресурса приведены в работах [1, 2]. Подход, предложенный в [3], направлен на решение задачи распределения временного ресурса. Кроме того, в работах [4, 5] предложены подходы, направленные на совместное решение задачи распределения частотного и временного ресурсов, сформулированные как задачи распределения слотов и формирования па кетов данных нисходящего канала связи. Однако подходы, предложенные в работах [4, 5] носят эвристический характер.

На основе недостатков известных решений [1–5], сформулированы требования к перспек тивным решениям задачи распределения слотов и формирования пакетов данных в нисходящем канале технологии WiMAX: ориентация на эффективное использование частотного и временно го ресурсов;

учет требований по скорости передачи ПС и качеству обслуживания;

минимизация количества служебных данных, передаваемых по каналу связи;

ориентация на преимущественно динамический характер решения задачи распределения слотов и формирования пакетов данных;

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

учет технологических особенностей сети (режима работы, ширины канала, количества подканалов, длительности кадра);

учет территориальной удаленности стан ций (определяет выбор схемы модуляции и кодирования (Modulation and Coding Scheme, MCS) для передачи сигнала пользовательской станции).

На основании проведенного анализа, а также сформулированных требований, предлагается единая модель, в которой описывается совместная процедура распределения как частотного, так и временного ресурса. Предложенная модель представлена в виде решения задачи распределения слотов и формирования их в пакеты данных нисходящего канала связи технологии WiMAX с учетом показателей качества обслуживания (требуемой скорости передачи ПС) и территориаль ной удаленности пользовательских станций сети (MCS для передачи сигнала пользовательским станциям).

Математическая модель распределения частотно-временного ресурса В предлагаемой модели предполагаются известными следующие исходные данные: N – коли чество пользовательских станций в беспроводной сети;

K – количество подканалов в частотном канале, определяемое используемым подрежимом DL FUSC или DL PUSC;

L – количество симво лов в кадре;

Rтрб – требуемая скорость передачи данных для обслуживания n-й ПС (Мбит/c);

S n – количество символов, формирующих один слот;

виды MCS, в зависимости от территориальной удаленности ПС;

M – количество слотов на одном подканале нисходящего канала для передачи полезной информации;

Q – количество слотов, предназначенных для передачи служебной инфор мации.

В ходе решения задачи распределения слотов и формирования пакетов данных, для переда чи полезной информации каждой ПС сети в нисходящем канале связи, в рамках предлагаемой модели необходимо обеспечить расчет булевой управляющей переменной (xn ), определяющей k,m закрепление подканалов и распределение слотов за ПС, на которых будут передаваться данные в нисходящем канале:

1, если m-й слот на k-м подканале выделен n-й ПС;

xn = (1) k,m 0, в противном случае.

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

1) Условие закрепления k-го подканала на протяжении передачи m-го слота не более чем за одной ПС N xn 1(k = 1, K;

m = 1, M );

(2) k,m n= 2) Условие закрепления за n-й ПС количества слотов, обеспечивающего необходимую ско рость передачи при используемом MCS M K n xn Rтрб (n = 1, N ), n RS (3) k,m m=1 k= SRn kn K где RS = (Tb +Tg )L+TRT G +TT T G – пропускная способность слота, закрепленного за n-й ПС, которая n s cb зависит от используемой MCS и представляет собой количество переданных бит за время равное длительности слота;

Rc – скорость кода, используемого при кодировании сигнала n-й ПС;

kb – n n битовая загрузка символа n-й ПС;

TT T G =105,7 мкс – длительность интервала переключения с приема на передачу;



Pages:     | 1 | 2 || 4 | 5 |   ...   | 11 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.