WWW.LIBRUS.DOBROTA.BIZ
БЕСПЛАТНАЯ  ИНТЕРНЕТ  БИБЛИОТЕКА - собрание публикаций
 

Pages:     | 1 ||

«образовательное учреждение высшего профессионального образования «Пермский национальный исследовательский политехнический университет» И.Ю. Зубко, Н.Д. Няшина ...»

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

Для исключения больших «выбросов» при определении удельной потенциальной энергии U для некоторой амплитуды возмущений A 0 атомы отклонялись не от своих начальных положений, соответствующих значению межатомного расстояния a0 и амплитуде A = 0, а от положений в решетке с параметром a AA, полученным при амплитуде A A, где A – шаг увеличения амплитуды. Например, сначала исследуется амплитуда 0,1a на исходной решетке со стандартным межатомным расстоянием a0. При этом определяется равновесное расстояние a0,1 = k a0, k 1, соответствующее этой амплитуде. В следующем опыте для большего значения амплитуды 0,2a берется уже не исходная решетка с параметром a0, а решетка с межатомным расстоянием k a0, полученным ранее. Такая процедура позволяет исключить излишнее сближение отдельных атомов с ростом случайных отклонений их положений в решетке, то есть с ростом температуры (амплитуды теплового движения) и, следовательно, уменьшить неизбежное во многих случайных процессах число «выбросов» значений удельной потенциальной энергии взаимодействия .

Зависимость безразмерных величин U / и A / для ГЦК-решетки при разных значениях n числа атомов на ребре образца (от 2 до 12) приведена на рис. 2.30. Для реализации алгоритма, отработанного в пакете «Mathematica», далее использовались языки программирования Intel Fortran Composer XE и PGI Accelerator Fortran v.12. Вычислительные эксперименты проводились на рабочей станции с процессором AMD Phenom X6 1055T и графическим ускорителем Nvidia GeForce GTS460 .

Расчеты показали более чем двукратное преимущество вычислений на графическом ускорителе (с использованием технологии PGI Accelerator) по сравнению с наиболее эффективными реализациями OpenMP-распараллеливания как для переменных одинарной, так и для переменных двойной точности .

Рис. 2.30. Зависимость безразмерной потенциальной энергии образца с ГЦК-решеткой от безразмерной амплитуды теплового движения атомов при различных размерах образца Результаты показывают, что все кривые удельной потенциальной энергии пересекаются в одной точке с критическим значением амплитуды A* 0,46. Кривые для разных значений n довольно быстро сходятся к некоторой кривой, которую можно отождествить с макроскопической. При условии A A* наименьшей удельной потенциальной энергией обладает самый большой из рассмотренных образцов. После точки пересечения кривых картина изменяется, и меньшая удельная энергия соответствует наименьшему из возможных образцов с 2 атомами на ребре. Это значение амплитуды можно сопоставить с температурой плавления. Удельная кинетическая энергия системы атомов, совершающих колебания с одинаковой частотой и амплитудой A, осредненная по периоду колебаний, равна A22 / 2. Принимая, что частота колебаний не зависит от амплитуды и является параметром материала, ее можно найти из условия

A*22 = 2k T*,

где T* – температура плавления выбранного материала .

Задания:

1. Провести исследование изменения равновесного межатомного расстояния и удельной потенциальной энергии для образцов различного размера с ОЦК-решеткой .

2. Найти такие значения n1 и n2 числа атомов на ребре кубического образцов с ГЦК- и ОЦК-решеткой соответственно, что общее количества атомов в кубических образцах с такими решетками будут мало отличаться .

3. Изобразить на одном рисунке графики зависимостей удельной потенциальной энергии образцов с ГЦК- и ОЦК-решетками атомов одного сорта от амплитуды теплового движения при найденных значениях n1 и n2 числа атомов на ребре образца. Исследовать возможность появления значений амплитуды, соответствующих фазовым переходам .





4. Провести исследование изменения равновесного межатомного расстояния и удельной потенциальной энергии для образцов различного размера с ГПУ-решеткой .

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

Список литературы к главе 2

1. Атомно-дискретное описание влияния анизотропных межатомных взаимодействий на упругие свойства ГПУ металлов / М.А. Баранов [и др.] // Физика твердого тела. – 2004. – Т. 46. – Вып. 2. – С. 212–217 .

2. Бертяев Б.И., Реут И.И. Об уравнении состояния, сжимаемости и внутреннем давлении в металлах с ОЦК-, ГЦК- и ГПУ-решeтками // Вестн. Самар. гос. техн. ун-та. Сер. Физ.-мат. науки. – 2008. – № 2 (17). – С. 215–223 .

3. Васильев Д.М. Физическая кристаллография. – М.: Металлургия, 1981. – 248 с .

4. Зубко И.Ю., Трусов П.В. Определение упругих постоянных ГЦКмонокристаллов с помощью потенциала межатомного взаимодействия // Вестник ПНИПУ. Механика. – Пермь: Изд-во Перм. нац. исслед. политехн. ун-та, 2011. – № 1. – С. 147–169 .

5. Вывод упругого закона монокристаллов металлов из потенциала межатомного взаимодействия / И.Ю. Зубко [и др.] // Вестник Нижегородского университета им. Н.И. Лобачевского. – Н. Новгород: Изд-во ННГУ им. Н.И. Лобачевского, 2011. – № 4. Ч. 5. – С. 2181–2183 .

6. Кривцов А.М. Деформирование и разрушение твердых тел с микроструктурой. – М.: Физматлит, 2007. – 304 с .

7. Кривцов А.М. Упругие свойства одноатомных и двухатомных кристаллов: учебное пособие. – СПб.: Изд-во С.Петерб. гос. политехн. ун-та, 2010. – 144 с .

8. Кривцов А.М., Подольская Е.А. Моделирование упругих свойств кристаллов с гексагональной плотноупакованной решеткой // Механика твердого тела. – 2010. – № 3. – С. 77–86 .

9. Поздеев А.А., Трусов П.В., Няшин Ю.И. Большие упругопластические деформации: теория, алгоритмы, приложения. – М.: Наука. 1986. – 232 с .

10. Симонов М.В., Зубко И.Ю. Определение равновесных параметров решетки различных ГПУ-монокристаллов с помощью потенциала межатомного взаимодействия Ми. // Вестник ПНИПУ. Механика. – Пермь:

Изд-во Перм. нац. исслед. политехн. ун-та, 2012. – № 3. – С. 205–218 .

11. Метод молекулярной динамики в физической химии / под ред .

Ю.К. Товбина. – М.: Наука, 1996. – 334 с .

12. Черных К.Ф. Введение в анизотропную упругость. – М.: Наука, 1988. – 190 с .

ГЛАВА 3. МЕТОД КЛЕТОЧНЫХ АВТОМАТОВ

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

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

3.1. МОДЕЛИРОВАНИЕ С ИСПОЛЬЗОВАНИЕМ ИМИТАЦИОННОГО ПОДХОДА

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

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

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

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

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

2. Если аналитические методы имеются, но математические процедуры трудно реализуемы, сложны и трудоемки .

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

4. Когда имитационный подход оказывается единственным способом исследования сложной системы из-за невозможности наблюдения явлений в реальной обстановке .

5. Когда необходимо контролировать протекание процессов в сложной системе путем замедления или ускорения процессов в ходе имитации .

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

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

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

Имитационный подход, в частности, оправдан, если вопросы, на которые должна ответить модель, относятся не к выяснению фундаментальных законов и причин, определяющих динамику реальной системы, а к анализу поведения системы, как правило, выполняемому исключительно в практических целях [3] .

Суть подхода, используемого при разработке имитаторов, состоит в том, что процесс функционирования сложной системы представляется в виде определенного алгоритма, реализуемого на компьютере. Для систем, состоящих из множества элементов, приходится строить модель не только всей системы, но и модели отдельных элементов. Как и для аналитического подхода, разработка имитатора ведется с использованием некоторой совокупности гипотез. Изменение даже одной гипотезы для одного из элементов системы может привести к необходимости пересмотра всей модели системы и поиску новых методов исследования (именно поэтому аналитические и численные подходы к моделированию сложных систем применяют после длительного и всестороннего изучения поведения как всей системы, так и отдельных ее элементов). Имитационный подход позволяет «максимально использовать всю имеющуюся в распоряжении исследователя информацию о системе» [3] .

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

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

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

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

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

При моделировании обычно используется три представления времени:

– реальное время моделируемой системы;

– модельное время, по которому организуется синхронизация событий в системе;

– машинное время имитации, отражающее затраты ресурсов времени ЭВМ на организацию имитации .

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

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

Процесс построения имитатора можно представить как технологический процесс, в котором можно выделить восемь этапов его построения [6]:

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

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

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

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

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

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

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

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

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

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

Успех или неудача проведения имитационных экспериментов с моделями сложных систем существенным образом зависят от инструментальных средств, используемых для моделирования, то есть от набора аппаратно-программных средств, предоставляемых пользователю-разработчику или пользователю-исследователю имитатора. В настоящее время существует большое количество специальных языков создания имитаторов на ЭВМ, которые называют языками моделирования. Перед разработчиком возникает проблема выбора языка, наиболее эффективного для целей моделирования конкретной системы. Языки моделирования заслуживают пристального внимания [8], так как, во-первых, число существующих языков и систем моделирования превышает несколько сотен и необходимо научиться ориентироваться в них. Во-вторых, почти каждый новый язык моделирования является не только средством, облегчающим доведение концептуальной модели до готовой машинной моделирующей программы, но и представляет собой новый способ «видения мира», то есть построения моделей реальных систем. К наиболее известным языкам моделирования систем с дискретными событиями относят SIMULA, SIMSCRIPT, GPSS, SOL, CSL. Достаточно полный обзор разработанных к середине 90-х годов языков моделирования, их классификация, обсуждение достоинств и недостатков, сравнение между собой представлены в [8] .

Отметим еще раз достоинства и недостатки имитационного подхода .

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

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

3.2. МЕТОД КЛЕТОЧНЫХ АВТОМАТОВ. ОСНОВНЫЕ ОПРЕДЕЛЕНИЯ

Одним из простых, но эффективных с вычислительной точки зрения имитационных подходов к моделированию физических процессов является подход клеточных автоматов .

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

Q внутренних состояний и V выходных сигналов, а также две функции:

функция переходов и функция выходов. Функция переходов определяет, в какое состояние q’ (из множества возможных состояний Q) перейдет автомат, если он находится в состоянии q и на вход поступил сигнал a (q,a) = q’, а функция выходов показывает, какой при этом образуется выходной сигнал v (из множества выходных сигналов V): (q,a) = v. Множество сигналов и состояний дискретны, и, кроме того, дискретно множество моментов, в которые поступают входные сигналы, выдаются выходные сигналы и меняются состояния. В этом случае число возможных значений аргументов функций переходов и выходов также конечно, и эти функции можно задавать таблично. Такие автоматы называются конечными автоматами .

–  –  –

Клеточные автоматы (КА) являются частным случаем [9] конечных автоматов, используемых для моделирования динамического поведения однородных сред. При этом пространство и время считаются дискретными, а физические величины в каждой точке моделируемой среды могут принимать конечное множество дискретных значений. Для КА существует достаточно развитая теория. В ее основу легли работы Дж. фон Неймана, который в 1948 году ввел в науку само понятие «клеточный автомат»

при разработке первой компьютерной модели биологического самовоспроизводства. Клеточный автомат можно представить как регулярную решетку «клеток», каждая из которых может находиться в конечном числе возможных состояний, например 0 или 1. Состояние системы полностью определяется значениями переменных в каждой клетке. Особенностями классических клеточных автоматов являются:

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

– переменные в каждой ячейке изменяются одновременно («синхронно»), исходя из значений переменных на предыдущем шаге;

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

Широко распространены [9] автоматы для моделирования динамического поведения двумерных сред. Каждой частице такой среды, занимающей некоторое пространство, ставится в соответствие элементарный автомат, имеющий форму квадрата (треугольника или шестиугольника) и именуемый клеткой. Совокупность клеток образует клеточное пространство, в котором функционирует КА. Как отмечалось выше, отдельная клетка имеет конечный набор состояний Q, а выходные сигналы клеток есть номера их состояний. Функция переходов определяется текущим состоянием q клетки, а также состоянием ее окружения. В зависимости от свойств моделируемого объекта выбирают различные типы окрестностей (рис. 3.2). В простейших случаях в качестве окрестности принимают четыре или восемь ближайших клеток. Однако возможны ситуации, когда состояние клетки зависит от состояния более удаленных клеток. Чтобы исключить особый вид окрестности для клеток, лежащих в крайних рядах, область замыкают, то есть для клеток крайнего левого столбца соседями слева считают клетки крайнего правого столбца, и наоборот. Аналогичное замыкание выполняют для клеток верхнего и нижнего ряда области .

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

Рис. 3.2. Типы плоских решеток и окрестности

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

Одним из подходов к построению таких правил является введение правил клеточного автомата при разделении поля автомата на подобласти.

При этом:

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

2. Формулируются правила для блока, задающие обновление содержимого всего блока (не отдельной клетки, как в обычном клеточном автомате) .

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

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

Рассмотрим простейшую схему разбиения:

массив клеток разбит на блоки размером 22;

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

Такую схему разбиения называют окрестностью Марголуса .

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

а б

Рис. 3.3. Блоки 2x2 окрестности Марголуса:

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

При описании диффузии с помощью подхода клеточных автоматов отдельные клетки поля моделирования считаются либо пустыми, либо содержащими частицу примеси (могут рассматриваться частицы разных типов). На каждом из двух проходов каждый блок независимо от соседних случайным образом «поворачивается» на четверть оборота по или против часовой стрелки. В результате каждая частица совершает случайное блуждание по полю моделирования. Как показали многочисленные исследования такого автомата [9], поведение большого (достаточного для статистической обработки) числа частиц на большом количестве клеток демонстрирует очень высокую точность соответствия такой модели решению классического уравнения диффузии примеси с концентрацией c (выведено в разд.

1.2 пособия):

c = Dc, t

где – оператор Лапласа, а D – коэффициент диффузии. Переходя к безразмерному виду уравнения и вводя масштаб времени t0, получим t = t0, где – безразмерная положительная величина, характеризующая время .

Для расстояний ведем масштаб L, тогда пространственные координаты будут выражаться через безразмерные параметры i как xi = i L. При этом оператор Лапласа выражается через безразмерный оператор как = L2. Тогда уравнение диффузии принимает вид

–  –  –

Полученный безразмерный комплекс D t0 D L2 определяется типом решетки автомата. Параметры t0 и L могут быть связаны с характерным временем и длиной процесса диффузии. В работе [7] для автомата, реализующего диффузию с помощью правил Марголуса, получены значения безразмерного коэффициента диффузии D для двумерного ( D = 3 / 2 ) и трехмерного ( D = 23 / 18 ) случаев .

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

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

3.3. ПРИМЕНЕНИЕ КЛЕТОЧНЫХ АВТОМАТОВ К ОПИСАНИЮ

ФОРМИРОВАНИЯ ДИСЛОКАЦИОННОЙ МИКРОСТРУКТУРЫ В МЕТАЛЛЕ

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

Дислокации могут быть краевыми, винтовыми и смешанного типа .

Дислокацию характеризуют два вектора: единичный вектор, направленный вдоль линии дислокации (на рис. 3.4, б – перпендикулярно плоскости рисунка) и особый (ненулевой только для линейных дефектов типа дислокаций) вектор – вектор Бюргерса b. Краевые дислокации представляют собой края лишних атомных полуплоскостей и для них b (рис. 3.4, б), для винтовых дислокаций b, а для смешанных взаимная ориентация векторов и b меняется вдоль линии дислокации. В зависимости от того, ниже плоскости залегания расположена лишняя атомная полуплоскость или выше, то есть в зависимости от знака вектора Бюргерса, краевые дислокации условно разделяют на положительные и отрицательные и обозначают символами и T. Линии дислокаций из геометрических соображений не могут обрываться внутри кристалла, а только на поверхности кристалла или на других дислокациях, образуя узлы дислокационной сетки. При низких гомологических температурах основным механизмом движения краевых дислокаций является их скольжение в плоскости залегания (вдоль горизонтали на рис. 3.4). При встрече краевых дислокаций противоположных знаков выстраивается полная атомная плоскость, в этом случае говорят, что дислокации аннигилируют .

а б Рис. 3.4. Контуры Бюргерса: а – в идеальном кристалле;

б – в кристалле с краевой дислокацией; мера несовпадения – вектор Бюргерса b Теоретически рассчитано, что для сдвига одной части кристалла относительно другой в бездефектном кристалле необходимо приложить касательное усилие, большее µ / 2, где µ – модуль сдвига. При наличии же краевой дислокации достаточно приложения существенно меньшего напряжения (Пайерлса) * = 10–5…10–4 µ для незначительного смещения атомов в ядре дислокации, в результате которого «перекидываются»

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

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

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

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

В реальных металлах с гранецентрированной кубической решеткой краевые дислокации, как правило, наблюдаются в сочетании с поверхностными дефектами упаковки. Длина вектора Бюргерса краевой дислокации не может быть меньше межатомного расстояния, иначе после скольжения дислокации не будет восстановлена правильная структура кристалла. Если же вдоль части АА плоскости скольжения совершить сдвиг на вектор b, меньший вектора b решетки, то линия А (рис. 3.6, а) будет границей области, по которой прошел сдвиг b, а линия А – границей сдвига b. На линии А получается скачок вектора сдвига b = b – b и дислокация с вектором Бюргерса b расщепляется на две частичные дислокации с векторами Бюргерса b и b. В кристалле на полосе АА нарушается правильная укладка структуры. Укладка атомов при этом может и отличаться от равновесной, но существуют величины сдвигов b и b, обеспечивающие минимум энергии в кристалле с локальным нарушением укладки. Такие поверхностные дефекты называются дефектами упаковки .

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

Рис. 3.5. Геометрическая модель структуры идеальной устойчивой полосы скольжения (а) и микрофотография (просвечивающая электронная микроскопия) структуры с мультипольными стенками (б) Дислокационными реакциями называются перестройки дефектной структуры кристалла, в результате которых одна или несколько дислокаций с векторами Бюргерса b1, b2, … превращаются в другие дислокации с векторами Бюргерса b1, b2, …. Полный вектор Бюргерса в результате реакций сохраняется. Примером реакции является аннигиляция дислокаций противоположных знаков. Другой простейшей реакцией является расщепление полной дислокации на частичные. В кристалле с гранецентрированной кубической решеткой происходят простые реакции дислокаций, лежащих в плоскости (111), с векторами Бюргерса b1 = [110], b2 = [101], b3 = [011] (рис. 3.6, б), образующими равносторонний треугольник: b1 + b2 b3. Если же в реакции участвуют расщепленные дислокации, то в результате реакции возникает прочный барьер ЛомераКоттрелла, блокирующий обе пересекающиеся плоскости скольжения (рис. 3.6, в). Эти барьеры играют большую роль в образовании прочных структур при множественном скольжении .

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

Время идет дискретными шагами. Рассматриваются прямолинейные краевые дислокации, бесконечно протяженные в третьем направлении. Все дислокации могут быть либо расщепленными, либо полными. Если дислокации расщеплены, то в результате дислокационных реакций появляются барьеры, в противном случае – другие нерасщепленные дислокации. Переползание расщепленных дислокаций невозможно. На плоскости моделирования дислокации представлены своими следами. Клетка может быть занята дислокацией любой из трех систем скольжения любого знака, источником дислокаций или барьером (рис. 3.7); клетка обладает симметрией 3-го порядка. Для этого идеально подходят равносторонние шестиугольники (рис. 3.8, а) .

Рис. 3.7. Возможные состояния клетки: а – пустая клетка; б, в – положительная и отрицательная полные краевые дислокации; г, д – положительная и отрицательная расщепленные краевые дислокации; е – дислокационный барьер, перекрывающий соответствующие плоскости скольжения; ж – источник дислокаций из некоторой плоскости скольжения

–  –  –

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

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

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

–  –  –

где F – сила, вызывающая скольжение; Fn – переползание дислокации с координатами {x, y} из наклонной плоскости за счет взаимодействия с дислокацией из начала координат горизонтальной системы скольжения;

– угол наклона нормали наклонной системы к вертикали. Сила, действующая на выбранную дислокацию, определяется с помощью окрестности дальнодействия .

Рис. 3.9. Окрестность дальнодействия выбранной дислокации, находящейся в центре окрестности, показана толстой черной линией .

Выделена дислокация, попадающая в окрестность

–  –  –

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

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

–  –  –

г д е Рис. 3.12. Основные правила для скольжения: а – свободное скольжение дислокации; б – отложенное скольжение; в – аннигиляция; г и д – реакция дислокаций (в случае расщепленных дислокаций образуется барьер соответствующей системы скольжения); е – остановка дислокации при встрече барьера, источника или дислокации, реакция с которой невозможна. Правила инвариантны относительно одновременной смены знаков у объектов в клетках, отражения относительно вертикальной оси (движение в противоположном направлении) и поворота на ±120° (другие системы скольжения) Второй проход отличается от первого смещением окрестностей на одну клетку по горизонтали, то есть ситуация рис. 3.12, б первого прохода на втором будет соответствовать ситуации рис. 3.12, а. Третий и четвертый проходы отличаются от первого и второго соответственно смещением на одну клетку по вертикали. Окрестности в течение каждого прохода не перекрываются, что обеспечивает консервативность субстанции при движении. Полагается, что реакция двух дислокаций разных систем происходит только в том случае, когда сумма векторов Бюргерса реагирующих дислокаций принадлежит третьей системе. При взаимодействии нерасщепленных дислокаций в результате реакции появляется новая нерасщепленная дислокация, расщепленных – дислокационный барьер .

Реализация переползания проводится по очереди для разных систем скольжения с помощью двухклеточной окрестности в два (четный и нечетный) прохода (рис. 3.13, а–б), обеспечивающих консервативность движения. Темным цветом показан способ выбора окрестности для одной и той же дислокации для проходов разной четности. В расчетах перебираются все дислокации и в зависимости от прохода и положения текущей дислокации, выбирается окрестность. Правила для переползания дислокаций приведены на рис. 3.14. Аннигиляции и дислокационных реакций при переползании не происходит .

–  –  –

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

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

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

Пусть L – длина большей стороны шестиугольника, тогда длина и высота прямоугольника (см. рис. 3.15) суть 3 (L – 1) и 3 (L – 1) .

Сдвиговая деформация от скольжения на одну клетку одной дислокации bl b некоторой системы вычисляется как = 3 3 (L –1)2 = 3 3 (N –1)2 l,

–  –  –

где k – номер системы скольжения; mk – единичный вектор, направленный вдоль линии скольжения; nk – единичный вектор нормали к линии скольжения; k – суммарный пластический сдвиг по k-й системе. Изменение объема в результате переползания дислокаций не учитывается .

Принимается гипотеза аддитивности однородных упругих и осредненных по области пластических деформаций = el + pl. Это соответствует структурной схеме материала в виде последовательного соединения двумерного упругого и двумерного пластического элементов. В плоском случае это можно представить заключением пластического элемента в упругую матрицу (при параллельном соединении пластический элемент помещается на упругую подложку). Упругий элемент моделирует кристаллическую решетку. Пластические деформации определяются модулем клеточных автоматов .

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

В соответствии со значениями напряжений на текущей итерации изменяется конфигурация дислокаций и пересчитывается суммарная (упругопластическая) деформация образца. При жестком нагружении накладываемое приращение деформаций представляется как = el + pl, то есть изменение деформаций кристаллической решетки el = – pl .

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

= C: el = C: ( – pl ()), (3.3) где тензор четвертого ранга C определяет линейно-упругие свойства решетки. Решая итерационно уравнение (3.3), можно построить диаграмму деформирования при нагружении по любой двумерной траектории в пространстве деформаций. Приращение пластической составляющей полной деформации определяется по текущей итерации в модели клеточных автоматов .

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

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

–  –  –

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

–  –  –

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

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

В построенной модели удается учесть физическое время, используя экспериментальный факт о линейной зависимости скорости движения дислокации от действующей на эту дислокацию силы при не слишком больших скоростях движения дислокаций. Зная силы, действующие на каждую дислокацию, можно определить их скорости. Согласно правилам клеточного автомата за одну итерацию каждая дислокация смещается не более чем на одну клетку, поэтому, принимая, что «самая быстрая» дислокация должна переместиться не дальше чем на одну клетку, можно определить длительность текущей итерации. Дислокации, которые за это время не успевают покинуть свою клетку, не перемещаются. Таким образом, разные итерации длятся различное физическое время. Отметим, что физическое время перемещения дислокации на расстояние, равное размеру клетки, очень мало. Поэтому исследование медленного нагружения требует рассмотрения очень большого числа итераций и размеров области моделирования. Результаты, получаемые за разумное вычислительное время, соответствуют очень быстрому изменению приложенных воздействий. Совместное движение дислокаций и примесных атомов реализуется как последовательный пересчет их конфигураций. На каждом новом шаге в качестве начальных условий используются конфигурации дислокаций и примесных атомов, полученные на предыдущем шаге. Учет физического времени позволяет исследовать процессы нагружения при любой скорости. Конфигурации дислокаций и примесных атомов, приведенные на рис. 3.18, получены при «быстром» нагружении .

–  –  –

Рис. 3.18. Конфигурации дислокаций и примесных атомов: а – начальное случайное распределение дислокаций; б – формирование первых дислокационных субструктур, миграция к ним атомов; в – сформировавшиеся дислокационные субструктуры и атмосферы В начале нагружения дислокации образуют диполи, по мере роста напряжений начинают образовываться дислокационные стенки. При этом атомы примеси диффундируют к отдельным дислокациям и простейшим субструктурам, образуя вокруг них скопления в виде облаков. Система приходит в устойчивую конфигурацию, дислокации совершают лишь незначительные колебания, так что суммарное приращение деформации остается нулевым или меняется крайне незначительно. Концентрация примесных атомов в это время становится все более неоднородной – повышается их концентрация в местах скопления дислокаций. Время нахождения системы в этой конфигурации соответствует вертикальному (или почти вертикальному) участку кривой на диаграмме – (рис. 3.19) .

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

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

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

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

Вопросы для самопроверки

1. Что такое периодические граничные условия?

2. Каким образом реализуется случайное блуждание частицы в клеточном автомате с правилами Марголуса?

3. Что такое «мягкое» и «жесткое» нагружение?

4. Как определяется приращение тензора пластических деформаций, если известно приращение сдвигов по различным системам скольжения?

5. Каков механизм формирования атмосфер примесных атомов вблизи дислокационных субструктур?

6. Как влияют атмосферы примесных атомов на подвижность дислокаций?

7. Объяснить механизм движения диполей краевых дислокаций за счет градиентов напряжений .

Список литературы к главе 3

1. Бусленко В.Н. Автоматизация имитационного моделирования сложных систем. – М.: Наука, 1977. – 240 с .

2. Введение в математическое моделирование / В.Н. Ашихмин [и др.];

под ред. П.В. Трусова. – М.: Логос, 2005. – 440 с .

3. Горстко А.Б. Познакомьтесь с математическим моделированием. – М.: Знание, 1991. – 160 с .

4. Зубко И.Ю. Модель клеточных автоматов для монокристалла с дислокациями трех систем скольжения // Математическое моделирование: сб. науч. тр. Вып. 12. / Перм. гос. техн. ун-т. – Пермь, 2004. – С. 45–53 .

5. Зубко И.Ю., Трусов П.В. Численное моделирование формирования дислокационных субструктур в поле примесных атомов // Вестник

УГТУ-УПИ. Механика микронеоднородных материалов и разрушение:

сб. науч. тр. / Урал. гос. техн. ун-т. – Екатеринбург; 2006. – № 11 (82). – С. 70–76 .

6. Максимей И.В. Имитационное моделирование на ЭВМ. – М.: Радио и связь, 1988. – 232 с .

7. Малинецкий Г.Г., Степанцов М.Е. Моделирование диффузионных процессов с помощью клеточных автоматов с окрестностью Марголуса // Журнал вычислительной математики и математической физики. – 1998. – Т. 38, № 6. – С. 1017–1020 .

8. Советов Б.Я., Яковлев С.А. Моделирование систем. – М.: Высшая школа, 1998. – 319 с .

9. Тоффоли Т., Марголус Н. Машины клеточных автоматов. – М.: Мир, 1991. – 280 с .

ГЛАВА 4. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ

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

Метод конечных элементов впервые был применен в инженерной практике в начале 50-х годов XX века. Первоначально он развивался по двум независимым один от другого направлениям – инженерному и математическому. На раннем этапе формулировки МКЭ основывались на принципах строительной механики, что ограничивало сферу его применения. И только когда были сформулированы основы метода в вариационной форме, стало возможным распространение его на многие другие задачи. Быстрое развитие МКЭ шло параллельно с прогрессом современной компьютерной техники и ее применением в различных областях науки и инженерной практики .

Значительный вклад в разработку МКЭ был сделан Дж. Аргирисом .

Им впервые дана общая матричная формулировка расчета стержневых систем на базе фундаментальных энергетических принципов, определена матрица податливости, а также введено понятие матрицы жесткости (как обратной матрице податливости). Работы Дж. Аргириса и его сотрудников, опубликованные в 1954–1960 годах, стала отправной точкой для матричной формулировки известных численных методов и применения ЭВМ в расчетах конструкций .

Для развития МКЭ особое значение имели вариационные принципы механики и математические методы, основанные на этих принципах .

Дискретизацию задачи на основе вариационного метода Ритца впервые в 1943 г. применил Р. Курант. Лишь в 50-е годы появились аналогичные работы Ж. Поли, Ж. Герша и др .

Первая работа, в которой была изложена современная концепция МКЭ, относится к 1956 году. Американские ученые М. Тэрнер, Р. Клафф, Г. Мартин и Л. Топп, решая плоскую задачу теории упругости, ввели элемент треугольного вида, для которого сформировали матрицу жесткости и вектор узловых сил. Название м е т о д к о н е ч н ы х э л е м е н т о в ввел в 1960 году Р. Клафф. В период 1960–1965 годов опубликованы работы, в которых на основе вариационных принципов получены конечные элементы для решения задач изгиба плит, тонких оболочек, массивов. Среди них можно отметить работы Р. Мак-Лейа, Р. Мелоша, Дж. Бесселина, Ф. де Веубеке, М. Джонса, Т. Пиана. В 1967 году издана первая монография о МКЭ О. Зенкевича и И. Чанга, в которой изложены основы метода и области его применения .

К семидесятым годам относится появление математической теории конечных элементов. Здесь можно выделить труды И. Бабушки, Р. Галлагера, Ж. Деклу, Дж. Одена, Г. Стренга, Дж. Фикса. Значительный вклад в разработку теоретических основ МКЭ внесли и российские ученые .

В.Г. Корнеев указал на совпадение математической сущности МКЭ и ВРМ .

Сопоставление МКЭ с рядом вариационных методов приведено в трудах Л.А. Розина. Под руководством А.С. Сахарова разработана моментная схема конечных элементов .

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

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

4.1. ПРОЕКЦИОННЫЙ И ВАРИАЦИОННЫЙ ПОДХОДЫ

К ПОСТРОЕНИЮ РАЗРЕШАЮЩИХ СООТНОШЕНИЙ МКЭ

Под термином метод конечных элементов (МКЭ) подразумевается семейство методов, основанных на проекционных методах решения уравнений или вариационных методах минимизации функционалов .

–  –  –

по нормам пространств E или H частичных сумм этих рядов). Введем последовательности конечномерных подпространств En и Hn, натянутых на первые n векторов, и обозначим через Pn и Qn операторы проектирования на эти пространства (En = Pn E, Hn = Qn H) .

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

–  –  –

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

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

Определение 1 .

Оператор L: H H (H – гильбертово пространство) называется симметричным, если (Lu, v) = (u, Lv), u, v H .

Определение 2 .

Симметричный оператор L называется положительно определенным, если существует константа C 0, такая, что ( Lu, u ) C u .

–  –  –

то эта производная называется слабой вариацией или просто вариацией функционала F в точке y0 и обозначается V(y0, h) .

Теорема 1 (необходимые условия экстремума). Пусть y0 Y – точка экстремума функционала F(y) и существует слабая вариация F(y0) .

Тогда F(y0) = 0 .

Определение 4 .

Квадратичный функционал F(u) = (Lu, u) – 2 (u, f) для положительно определенного оператора L называется функционалом энергии .

Теорема 2. Чтобы элемент u0 D(L) из области определения оператора L доставлял минимум функционалу энергии, необходимо и достаточно, чтобы этот элемент удовлетворял уравнению Эйлера Lu0 = f (причем такой элемент будет единственным) .

Таким образом, для определенного класса операторов задача минимизации функционала и задача решения уравнения Lu = f эквивалентны .

Для положительно определенного оператора задача о минимуме функционала энергии всегда имеет решение .

Для каждого положительно-определенного оператора можно ввести энергетическое пространство .

Определение 5. Пусть L: H H положительно определенный оператор, действующий в полном гильбертовом пространстве H. Построим новое пространство с элементами из области определения оператора D (L) и скалярным произведением [u, v]L = (Lu, v). Если оно окажется неполным, то пополним его. Пополнение пространства HL называется энергетическим .

Минимум функционала будем искать в энергетическом пространстве .

Определение 6. Элемент из энергетического пространства u0 HL, доставляющий минимум функционала энергии, называется обобщенным решением уравнения Lu = f .

Метод Ритца минимизации функционала .

Метод Ритца (Релея – Ритца) один из методов приближенного решения задачи минимизации, который состоит в отыскании минимума функционала на конечномерном пространстве Hn. Элементы v Hn называются пробными функциями. Аппроксимацией Ритца называется функция

uh Hn, минимизирующая F на пространстве Hn:

–  –  –

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

Матрицы СЛАУ, получаемые в методах Галеркина и Ритца, оказываются заполненными, т.е. большинство элементов этих матриц отличны от нуля. Решение систем с заполненными матрицами делает оба метода затратными. Поэтому необходимо выбрать базисные функции таким образом, чтобы коэффициенты матрицы и правая часть системы легко вычислялись, а сама система была разреженной. Такой выбор базисных функций приводит к методу конечных элементов (МКЭ) .

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

Процесс реализации МКЭ включает следующую последовательность шагов:

1. Дискретизация области: построение сетки, задание свойств элементов. Область, на которой решается задача, покрывается непересекающимися подобластями простого типа, которые называются конечными элементами (КЭ). Вершины КЭ называются узлами. Узлы предназначены для описания геометрии элементов и для задания компонент решения (неизвестная величина определяется в узлах). Узлы могут быть внешними и внутренними. Внешние узлы лежат на границе КЭ и используются для соединения элементов друг с другом. Узлы могут располагаться между угловыми узлами. КЭ может иметь внутренние узлы, такие элементы, являясь элементами высокого порядка точности, обеспечивают более точное описание искомых функций. Компоненты неизвестной функции в узле называются степенями свободы. В зависимости от рассматриваемой задачи число степеней свободы в узле различно. Например, в задаче теплопроводности в каждой точке определяется скалярная величина температуры – одна степень свободы. В двумерной задаче теории упругости относительно перемещений в каждом узле определяется двумерный вектор перемещений, тогда в каждом узле КЭ определяется две степени свободы. В качестве степеней свободы могут фигурировать не только значения искомой функции, но и ее производные по пространственным координатам, как в задачах расчета мембран или оболочек. Кроме того необходимо задать свойства материала, из которого изготовлена конструкция. Например, для задач теории упругости для изотропного материала необходимо задать модуль Юнга и коэффициент Пуассона, в стационарных задачах теплопроводности – коэффициент теплопроводности .

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

3. Формирование СЛАУ с учетом вкладов от элементов и узлов, введение граничных условий в систему уравнений. Например, если задача решается с помощью метода Галеркина, формируются интегралы из условия ортогональности невязки и базисных функций. Если решается задача в вариационной постановке с помощью метода Ритца минимизации функционала, то СЛАУ получается после приравнивания к нулю производной функционала. Интегралы по области разбиваются на интегралы по элементам M

–  –  –

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

3. Показать, что решение уравнения двумерной стационарной теплопроводности с краевыми условиями = 1 на и n = h на q, h – функция от х, у минимизирует функционал

–  –  –

найти уравнение Эйлера (уравнение Эйлера описывает малые отклонения нагруженного троса, покоящегося на упругом основании жесткостью k) .

Показать, что этот функционал равен потенциальной энергии системы .

5. Для функционала

–  –  –

6. Какой оператор называется симметричным, положительно определенным? Приведите примеры .

7. Сформулируйте необходимые условия экстремума функционала .

8. Что называется функционалом энергии?

9. Чем является решение уравнения Эйлера для функционала энергии?

10. В чем суть метод Ритца (Релея–Ритца) минимизации функционала?

11. В чем суть метода конечных элементов?

12. Какова последовательность шагов при реализации метода конечных элементов?

4.2. ВИДЫ КОНЕЧНЫХ ЭЛЕМЕНТОВ И ИХ КЛАССИФИКАЦИЯ

–  –  –

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

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

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

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

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

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

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

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

Рис. 4.1. Типы конечных элементов: а – стрежневой (простой фермовый);

б – плоские (основные); в – трехмерные (основные); г – осесимметричный;

д – изгибаемый пластинчатый; е – осесимметричный тонкостенный оболочковый [4] Осесимметричные оболочковые элементы, изображенные на рис. 4.1, е, важны на практике так же, как и осесимметричные сплошные элементы .

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

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

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

Приведем их вид [3, 6]:

Линейные (два узла для каждого элемента)

–  –  –

женного значения искомой функции в узлах элемента .

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

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

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

–  –  –

2 = 2 1, 3 = 2( 3 ), (4.9) 4 = 0,25(15 4 18 2 + 3), 5 = 7 5 10 3 + 3 .

Функции 0 и 1 имеют вид (4.4). Производные этих функций ортогональны, что приводит к матрице разрешающих соотношений диагональной формы .

Еще один вид элементов – эрмитовы элементы. Базисные функции для них могут быть получены аналогичным образом как для лагранжевых элементов, но с использованием полиномов Эрмита вместо полиномов Лагранжа. При этом в каждом узле элемента будут определяться не только значения искомой функции, но и ее производные [5] .

Рассмотрим одномерный элемент с s узлами, причем узлы не обязательно расположены равномерно. Пусть у каждого узла две степени свободы – функция и ее производная s/x. Аппроксимация искомой функции на элементе имеет вид

–  –  –

Первый индекс у базисной функции ji обозначает порядок дифференцирования соответствующей узловой переменной, а второй – номер узла. Для того чтобы выполнялись условия

–  –  –

Узловые производные по координате x заменяются на производные по : / x = h 1 /, при этом величина h = x2 – x1 в последних соотношениях сокращается .

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

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

–  –  –

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

Такой элемент называется лагранжевым элементом первой степени .

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

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

–  –  –

но единице только в узле (r, s) и тождественно равно нулю во всех других узлах .

Рис. 4.2. Три элемента лагранжева семейства и некоторые соответствующие базисные функции: а – элемент первого порядка;

б – элемент второго порядка; в – элемент третьего порядка Отметим, что поскольку базисные функции (4.12) определены произведением многочленов степени р отдельно от х и у, то число членов в этом произведении превышает число членов, необходимое для воспроизведения многочлена степени р от х и у .

На рис. 4.3, а построен треугольник Паскаля, по которому можно определить число членов в многочлене произвольной степени (например в кубическом многочлене – 10 членов). На рис. 4.3, б дополнительно выделены шесть членов, появляющихся в кубической лагранжевой базисной функции, определенной соотношением (4.12). Можно уменьшить число узлов, ассоциированных с элементами высокого порядка, обеспечив, чтобы базисные функции элемента содержали бы только члены, входящие в многочлен соответствующей степени. Это семейство так называемых серендиповых элементов (рис. 4.3, в). Здесь узлы располагаются (по возможности) на границах элементов, а базисные функции (для узлов, расположенных не в вершинах элемента) получаются умножением членов степени р по одной переменной на линейные члены по второй переменной (рис. 4.4, а, б) .

Функции формы для четырехугольных серендиповых элементов 2-го и 3-го порядка, содержащие 8 и 12 узлов (см. рис. 4.4) соответственно, представляют собой полиномы = 1 + 2 + 3 + 4 + 5 2 + 6 2 + 7 2 + 8 2, = 1 + 2 + 3 + 4 + 5 2 + 6 2 + 7 2 + 8 2 + + 9 3 + 10 3 + 11 3 + 12 3 .

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

Рис. 4.3. Треугольник Паскаля: а – заштрихованы элементы, образующие полный кубический полином; б – заштрихованы элементы, образующие кубический лагранжев полином; в – заштрихованы элементы, образующие кубический серендипов полином [3] Рис. 4.4. Элементы серендипова семейства и некоторые соответствующие базисные функции: а – элемент второго порядка; б – элемент третьего порядка [8] Итак, общая формула для четырехугольных серендиповых элементов 2-го и 3-го порядка имеет вид [9]

–  –  –

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

–  –  –

в случае квадратичного элемента .

Пример. Построить функцию формы для квадратичного четырехугольного элемента для 0-го узла (см. рис. 4.4, а) .

Начнем с определения произведения:

–  –  –

0 = (1 )(1 ) ( 1 + 2 + 3 ) .

Произведение (1 – ) (1 – ) тождественно обращается в нуль в узлах 3, 4, 5, 6, 7. Константы 1, 2 и 3 должны быть выбраны так, чтобы 0 была равна единице в нулевом узле и нулю в узлах 1 и 7, то есть 0 : 0 (1; 1) = 2 ( 1 2 3 ) = 1, 1: 0 (0; 1) = 2 ( 1 3 ) = 0, 7 : 0 (1; 0) = 2 ( 1 2 ) = 0 .

–  –  –

Окончательно функция формы для 0-го узла имеет вид 0 = 1 / 4(1 )(1 ) (1 + + ) .

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

Определяем произведение (так как узел принадлежит только одной – первой стороне):

–  –  –

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

Двумерные базисные функции для треугольников. L-координаты .

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

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

Рис. 4.6. Семейство треугольных элементов: а – треугольник Паскаля;

б – элемент первого порядка; в – элемент второго порядка;

г – элемент третьего порядка

Рассмотрим построение двумерного линейного треугольного элемента. Линейный треугольный элемент содержит три узла i, j, и k, расположенных в вершинах треугольника, перенумерованных против часовой стрелки (рис. 4.7). Координаты узлов в системе координат хОу:

(xi, yi), (xj, yj), (xk, yk) соответственно. Функции формы для каждого узла записываются в виде

–  –  –

Для определения неизвестных шести констант необходимо решить СЛАУ из шести уравнений, полученных из условия равенства базисной функции 1 в своем узле и равенства нулю в остальных узлах:

–  –  –

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

Запишем линейную зависимость между локальными координатами

L1, L2, L3 и декартовыми координатами произвольной точки внутри треугольника:

x = L1 xi + L2 x j + L3 xk, x = L1 yi + L2 y j + L3 yk, (4.21) 1 = L1 + L2 + L3 .

Здесь (xi, yi), (xj, yj), (xk, yk) – декартовы координаты вершин треугольника. Очевидно, что функция L1 должна принимать значение 1 в i-м узле и обращается в 0 в остальных вершинах – совершенно так же, как базисные функции линейного треугольного элемента .

Несложно показать, что L1 выражаются через координаты вершин треугольника в виде соотношений

–  –  –

где n – порядок треугольного элемента (число узлов на стороне –1);

F = Li c – функция L-координат, определяется из уравнений n линий уровня функций L1, L2, L3, которые проходят через все узлы, за исключением узла, для которого определяется функция формы; F ( L1, L2, L3 ) – значение этих функций в узле, для которого составляется функция формы. Соотношение (4.25) является полиномом Лагранжа, записанного в L-координатах .

Пример. Составить функцию формы для 1-го узла кубического элемента .

Порядок элемента равен числу узлов на стороне треугольника минус 1:

–  –  –

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

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

1.0 1.0 0.5 0.0 0.5 0.0 0.5

–  –  –

0.0 1.0 0.0 0.0 1.0 0.5 0.5 0.5 1.0 0.0 0.0

–  –  –

0.5 0.0 1.0 1.0 1 .

0.5 0.5 0.0 0.0 0.5 0.0 0.0 0.5 0.5 1.0

–  –  –

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

На рис. 4.12 показаны типичные иерархические квадратичная и кубическая базисные функции .

–  –  –

Функции формы степени p лагранжева семейства для прямоугольных параллелепипедов в локальной нормированной системе координат ( 1,, 1 ) могут быть записаны как обобщение на трехмерный случай соотношений (4.12):

–  –  –

Здесь (xi, yi), (xj, yj), (xk, yk), (xl, yl) – декартовы координаты вершин тетраэдра. Очевидно, что функция L1 должна принимать значение 1 в i-м узле и обращается в 0 в остальных вершинах – совершенно так же, как базисные функции линейного тетраэдрального элемента. Можно показать, что i = L1, j = L2, k = L3, l = L4 .

–  –  –

4.2.4. Разбиение области на элементы .

Требования к конечным элементам Дискретизация области. Процесс дискретизации можно разделить на два этапа: разбиение расчетной области на элементы и нумерация элементов и узлов. Сложности второго этапа связаны с требованиями повышения эффективности вычислений [9] .

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

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

Наиболее просто разбить треугольную подобласть на элементы, если выбрать определенное число узлов вдоль каждой стороны, соединить соответствующие узлы прямыми линиями и точки пересечения этих линий считать узлами (рис. 4.14). Узлы на сторонах подобласти необязательно располагать на равных расстояниях. Варьирование расстояния между ними позволяет изменять размеры элементов. Если треугольная подобласть криволинейная, то криволинейные границы элементов заменяются на прямые отрезки. Если на стороне треугольной подобласти выбрано п узлов, число треугольных элементов в результате разбиения равно (п – 1)2 .

Четырехугольные зоны обычно разбиваются на элементы соединением узлов на противоположных сторонах. Пересечения линий определяют внутренние узловые точки (рис. 4.15). Внутренние четырехугольники могут рассматриваться как элементы; но они могут быть разбиты и на треугольные Рис. 4.14. Деление области Рис. 4.15. Деление области в виде треугольного вида на линейные четырехугольника на линейные треугольные элементы [9] треугольные элементы [9] элементы проведением короткой диагонали в каждом внутреннем четырехугольнике. Разбиение с использованием короткой диагонали предпочтительно, потому что элементы, близкие по форме к равностороннему треугольнику, приводят к более точным результатам, чем длинные узкие треугольники .

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

Для оценки качества элементов вводится параметр – характерный размер сетки diam ((i)). Под этим параметром будем понимать некий характерный линейный размер, например максимальную диагональ элементов или их максимальную высоту, при стремлении которого к нулю размеры всех конечных элементов также стремятся к нулю. Кроме того, введем величину ((i)) – наибольший из радиусов, вписанных в конечные элементы окружностей. Качество конечно-элементной сетки из прямосторонних треугольных или четырехугольных элементов оценивается величиной ((i))/diam((i)). Причем при прочих равных условиях погрешность аппроксимации тем хуже, чем это соотношение меньше. Поэтому построение сеток следует осуществлять так, чтобы величина ((i))/diam ((i)) была по возможности больше или хотя бы ((i)) / diam ((i)) с0 0, i = 1, 2, …, Nэл, равномерно по всем элементам .

На практике удобно вычислять эквивалентные им величины 2mes((i)) / aibi или 2 mes((i)) / bici, где ai bi ci – длины сторон треугольника, а mes ((i)) – его площадь, подсчитанная с помощью определителя, составленного из координат узлов элемента. Идеальным элементом для такого критерия является правильный треугольник. Следует избегать сильно вытянутых треугольных элементов с углом менее 30° .

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

6mes((i)) / abc, где mes((i)) – объем тетраэдра, а величина abc – наименьшее из произведений длин ребер, исходящих из одной вершины. Идеальным элементом с точки зрания такого критерия является правильный тетраэдр .

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

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

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

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

Наиболее простой способ существенного изменения размеров элементов Рис. 4.16. Деление расчетной области на треугольные и четырехугольные зоны с последующим разбиением на треугольные элементы [9] заключается в применении четырехугольных подобластей с неравным числом узлов на противоположных сторонах. Хорошим вариантом является случай расположения двух узлов на одной стороне против каждых трех узлов на противоположной стороне (рис. 4.17) .

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

если узел может перемещаться только в одном направлении, используется

–  –  –

символ подвижного шарнира (рис. 4.18, б). Изображенные подвижные шарниры допускают перемещение в вертикальном направлении и не позволяют двигаться в горизонтальном направлении) .

Нумерация узлов. Нумерация узлов влияет на эффективность вычислений, необходимых для получения решения. Использование МКЭ приводит к СЛАУ, большое число коэффициентов равно нулю. Все ненулевые коэффициенты матрицы и некоторые нулевые находятся между двумя линиями, параллельными главной диагонали (рис. 4.19). Расстояние между этими линиями и главной диагональю называется шириной полосы матрицы. Все коэффициенты вне этой полосы равны нулю, и они не должны храниться в машинной памяти. Эффективная вычислительная программа использует только те коэффициенты матрицы, которые находятся внутри указанной полосы. Уменьшение ширины полосы приводит к сокращению размеров требуемой машинной памяти, а также к сокращению времени вычислений. Ширина полосы В вычисляется по формуле [9] B = (R + 1) Q, где R – максимальная по элементам величина наибольшей разности между номерами узлов в отдельном элементе; Q – число неизвестных (число степеней свободы) в каждом узле. Минимизация величины B связана с минимизацией величины R, что, в частности, может быть осуществлено последовательной нумерацией узлов при движении в направлении наименьшего размера тела .

Два разных способа нумерации узлов показаны на рис. 4.20 .

–  –  –

Для разбиения, представленного на рис. 4.20, а, величина R – максимальная по элементам величина наибольшей разности между номерами узлов в отдельном элементе, равна 3. Тогда ширина полосы для задачи теплопроводности (где в каждом узле отыскивается одна неизвестная величина – одна степень свободы) равна 4, а для задачи упругости (две степени свободы в каждом узле) – 8 .

Для разбиения, представленного на рис. 4.20, б, величина R равна 6 .

Тогда ширина полосы для задачи с одной степенью свободы в узле равна 7, а для задачи упругости (две степени свободы в каждом узле) – 14 .

Правильная нумерация узлов в этом примере сокращает машинную память почти на 50 % .

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

Построение глобальной матрицы из набора локальных называется ансамблированием (конденсацией, сборкой). Рассмотрим метод прямой жесткости [9] на примере двумерной задачи теории упругости. Предположим, что область разбита на 4 линейных треугольных конечных элемента, глобальные номера узлов обозначены арабскими цифрами, локальные номера – i, j, k, арабскими цифрами в кружках обозначены номера элементов (рис. 4.21). В каждом узле определяется вектор перемещений, состоящий из двух компонент, то есть по две степени свободы в каждом узле. Таким образом, необходимо найти 12 неизвестных узловых перемещений, и глобальная матрица будет содержать 12 строк и столбцов. Покажем, как включить локальные матрицы для элементов 1 и 2 в глобальную матрицу .

–  –  –

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

Задачи:

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

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

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

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

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

6. Разбить четырехугольник на 24 конечных элемента, используя пять узлов вдоль одной пары сторон и четыре узла вдоль другой пары. Пронумеровать узлы так, чтобы получить минимальное значение величины R .

7. Разбить прямоугольный треугольник примерно на 60 элементов, предварительно выделив две треугольные и одну четырехугольную подобласти. Поместить наименьшие по размеру элементы вблизи прямого угла .

Вопросы для самопроверки

1. К чему приводит кусочная определенность базисных функций в МКЭ?

2. Перечислите основные типы конечных элементов. В чем их различия, где они применяются?

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

4. Запишите вид функций формы для иерархических многочленов .

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

5. Для чего используется аппроксимация решения с помощью полиномов Эрмита?

6. Как построить функцию формы второго порядка для прямоугольного элемента лагранжева типа? Запишите функцию формы второго порядка для углового узла; для узла, расположенного в середине стороны;

для узла, расположенного в центре элемента .

7. В чем принципиальная разница между функциями формы для элементов лагранжева и серендипова типов?

8. Как построить функцию формы второго порядка для прямоугольного элемента серендипова типа? Запишите функцию формы второго порядка для углового узла; для узла, расположенного в середине стороны .

9. Что такое L-координаты для треугольного конечного элемента?

10. Как построить квадратичные и кубические функции формы на треугольнике с помощью L-координат?

11. Как построить иерархические функции формы для треугольного конечного элемента с помощью L-координат?

12. Как построить функции формы лагранжева типа для четырехугольного элемента?

13. Как построить функции формы серендипова типа для четырехугольного элемента?

14. Как вводятся трехмерные L-координаты? Как с их помощью записать функции формы на тетраэдре?

15. Каким образом разбить треугольную или четырехугольную форму на треугольные конечные элементы?

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

17. Как размер ширины ленты глобальной матрицы жесткости зависит от нумерации узлов конечных элементов в расчетной области?

18. Как на схеме обозначаются закрепленные узлы, узлы на границах симметрии?

19. Опишите процедуру построения глобальной матрицы жесткости с помощью метода прямой жесткости .

4.3. ОТОБРАЖЕНИЕ И ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ

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

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

–  –  –

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

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

Для этого элемент площади в глобальных координатах заменяется элементом площади в локальных координатах:

d x d y = det Jd d. (4.30) Тогда любой интеграл по четырехугольному элементу

–  –  –

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

–  –  –

Здесь (xl, yl) – глобальные координаты точки, в которую требуется отобразить узел l в плоскости (, ) .

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

Следует следить, чтобы параметрическое отображение не вырождалось, т.е. J 0, что происходит, если внутренний угол четырехугольника больше 180° или если при использовании квадратичного отображения расстояние между центральным и угловыми узлами меньше трети длины стороны четырехугольника .

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

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

При дифференцировании по L-координатам следует учесть, что только две из них являются независимыми (см. (4.213)). Пусть L1, L2 – независимые координаты, тогда i i x i y = +, L1 x L1 y L1 (4.33) i i x i y = +, L2 x L2 y L2

–  –  –

4.3.2. Субпараметрические, изопараметрические и суперпараметрические элементы С понятием параметрического отображения тесно связано деление элементов на субпараметрические, изопараметрические и суперпараметрические [9] .

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

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

Возможны три варианта соотношения между числами узлов в этих двух множествах:

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

4.5, а, б; 4.11, б, в; 4.12, а, б) .

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

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

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

Однако применение суперпараметрических элементов имеет ряд недостатков. В частности, локальные матрицы в этом случае определяются с помощью численного интегрирования (компоненты матрицы градиентов зависят от координат). Кроме того, необходимо очень аккуратно вводить обозначения узлов, так как такой элемент может содержать до 30 узлов, их локальные номера неудобно обозначать i, j, k, l,…и поэтому используются числовые значения. В случае двумерных элементов нумерация начинается в произвольной точке и соответствует обходу против часовой стрелки .

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

–  –  –

Эти соотношения можно рассматривать как систему (m + 1) нелинейных алгебраических уравнений относительно (2n + 2) неизвестных C0, C1, C2,...Cn, x0, x1, x2,...xn. Для разрешимости этой системы уравнений число уравнений должно быть равно числу неизвестных m + 1 = 2n + 2, откуда m = 2n + 1 .

Таким образом, нужно выбрать такую сетку n, содержащую n + 1 узел, чтобы точно численно проинтегрировать полином степени

m = 2n + 1:

–  –  –

Пример 1. На отрезке [– 1; 1] построить сетку, содержащую один узел x0 и определить коэффициент C0 для интегрирования полинома степени m = 2n + 1 = 1 .

Учитывая, что n = 0, а = –1, b = 1, из (4.38) получим систему

–  –  –

Пример 2. На отрезке [– 1; 1] построить сетку, содержащую два узла x0 и x1 определить коэффициенты C0, C1 для интегрирования полинома степени m = 2n + 1 = 3 .

Опять учтем, что n = 1, а = –1, b = 1, из (4.38) получим систему

–  –  –

Значения координат узлов (точек интегрирования) и соответствующих весовых коэффициентов в формулах Гаусса для достаточно больших п приводятся во многих книгах по численному анализу и МКЭ (табл. 4.1) .

–  –  –

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

Квадратурные формулы Гаусса в случае двух и трех переменных .

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

I = f (,)d d .

Так как интегрирование проводится по квадрату, то простейшим подходом будет осуществление численного интегрирования независимо от и .

Таким образом, начиная с вычисления внутреннего интеграла и используя формулу (4.37), получим

–  –  –

в которой координаты узлов ( k, j ) определяются аналогично табл. 1 .

Если формулы интегрирования точны отдельно по и для полинома степени m, то выражение (4.39) будет давать точные значения для всех выражений вида m1 m2, m1, m2 m. Стандартные квадратурные правила иллюстрируются на рис. 4.22 .

Рис. 4.22. Схема узлов квадратур Гаусса для четырехугольника [8]

–  –  –

Координаты узлов и весовые коэффициенты для треугольных элементов приведены в табл. 4.2 [3, 8, 9] .

Возможно обобщение на трехмерный случай. Правила для интегрирования по тетраэдрам приведены в табл. 4.3 [3, 8, 9] .

–  –  –

1, 0, 0 6 5 0, 1, 0 0, 0, 1 6 0 0,2250000000, 1, 1

–  –  –

Задачи:

1. В плоскости (х, у) координаты узлов четырехугольного элемента:

(0, 0), (1, 0), (0,4, 0,4), (0, 1). Построить изопараметрическое отображение на квадратный билинейный элемент, определенный условием 1, 1 .

Показать, что отображение является вырожденным .

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

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

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

Построить изопараметрическое отображение треугольного элемента, координаты узлов которого в плоскости (х, у) – (1, 1), (5, 1), (8, 2), (7,5, 4), (4, 5), (2, 3) на описанный выше равнобедренный прямоугольный треугольник .

5. При конечно-элементных вычислениях часто требуется находить e e интегралы вида I = k l m d x d y и I = k le e d x d y по элеменx x m e e ту е. Найти число узлов, требующееся для точного вычисления этих интегралов при использовании квадратур Гаусса, если е а) билинейный элемент с четырьмя узлами; б) серендипов элемент с восемью узлами;

в) лагранжев элемент с девятью узлами .

Вопросы для самопроверки

1. Что такое параметрическое отображение?

2. В чем особенность дифференцирования и интегрирования функций форм по локальным координатам? Как перейти к производным по глобальным пространственным переменным?

3. Как от дифференцирования по L-координатам в треугольном конечном элементе перейти к производным по глобальным пространственным переменным?

4. Что такое изопараметрические, суперпараметрические и субпараметрические элементы? Приведите примеры. В чем разница их применения?

5. В чем особенность применения квадратурных формул Гаусса для численного интегрирования полиномов?

6. Какая связь между степенью интегрируемого полинома и количеством точек интегрирования в квадратурных формулах Гаусса для функций одной переменной?

7. Как обобщить квадратурные формулы Гаусса для интегрирования функций двух и трех переменных?

4.4. РАЗРЕШАЮЩИЕ СООТНОШЕНИЯ МКЭ

ДЛЯ КВАЗИСТАТИЧЕСКИХ ЗАДАЧ ТЕОРИИ УПРУГОСТИ

Разрешающие соотношения МКЭ для квазистатической задачи теории упругости могут быть получены с помощью проекционных методов [6] или на основе вариационного подхода. В задачах теории упругости соотношения выводятся на основе принципа виртуальной работы .

Коротко рассмотрим применение вариационного подхода на примере пространственной задачи теории упругости [4] .

Пусть du – виртуальное перемещение элемента объема. Тогда виртензор туальные деформации в этом объеме d, где = 1 2 u + uT малых деформаций .

Согласно принципу виртуальной работы сумма потенциала внешних нагрузок dV и величины запасенной энергии деформации d U при виртуальных перемещениях du равны нулю:

dV + dU = 0 .

Поскольку величина запасенной энергии деформации

–  –  –

Теперь можно получить локальные разрешающие соотношения МКЭ (для одного элемента) в матричной форме .

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

Перемещения аппроксимируются на элементе с т узлами соотношением

–  –  –

где {u r } – вектор узловых перемещений, для пространственной задачи имеет три компоненты (для плоской – две) .

Аналогично для виртуальных перемещений

–  –  –

имеют размер 31 и встанут в I-ю строку глобального вектора нагрузок .

После ансамблирования, исключения поверхностных нагрузок, описывающих взаимодействие между элементами, и учета статических граничных условий получаем СЛАУ для всего рассматриваемого тела:

[ K ]{u} + {F}0 + {F}0 + {F} f + {F}t = 0, где [K], {F}f, {F}, {F}, {F} – соответственно глобальные матрица жестt кости и вектора нагрузок, обусловленных объемными силами, начальными напряжениями и начальными деформациями; вектор поверхностных нагрузок, учитывающий статические граничные условия t = n, x, y, z t .

Обозначим {F } = {F }0 + {F }0 + {F } f + {F }t, тогда глобальная система уравнений относительно узловых перемещений в окончательном виде [ K ]{u} + {F } = 0. (4.47) Полученная СЛАУ имеет симметричную ленточную вырожденную матрицу. Сумма элементов матрицы по строкам или столбцам равна нулю (выполнение этого может выступать некоторым критерием правильности построения глобальной матрицы). Учет заданных кинематических граничных условий сделает матрицу невырожденной .

Заметим, что статические граничные условия были учтены автоматически при построении системы: они вошли в вектор поверхностных нагрузок в виде {F }t = [i ]Т {F }d, где компоненты вектора {F} соответt ствуют заданному вектору t = n, x, y, z t .

Для того чтобы учесть кинематические граничные условия (условия на перемещения), необходимо преобразовать глобальную систему уравнений следующим образом:

Пусть ит – заданное перемещение по направлению т-й степени свободы (всего М степеней свободы), тогда

1) из вектора {F } вычитается т-й столбец матрицы [К], умноженный на ит:

<

–  –  –

Задачи:

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

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

3. Запишите вид матрицы градиентов для двумерной задачи с использованием линейных треугольных КЭ, квадратичных треугольных КЭ .

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

5. Пусть на границе расчетной области задана распределенная нагрузка t = {q, 0}, q = const. Получите соответствующий локальный вектор распределенных нагрузок для линейного треугольного КЭ, квадратичного треугольного КЭ .

6. Пусть на границе расчетной области задана распределенная нагрузка t = {0, q}, q( y ) = ay + b. Получите соответствующий локальный вектор распределенных нагрузок для линейного треугольного КЭ, квадратичного треугольного КЭ .

Вопросы для самопроверки

1. Получить разрешающие соотношения МКЭ для задачи теории упругости с помощью проекционного метода – метода Галеркина .

2. Сформулируйте принцип виртуальной работы для механической системы .

4.5. ПОСТАНОВКА НЕСТАЦИОНАРНЫХ И ДИНАМИЧЕСКИХ ЗАДАЧ

–  –  –

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

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

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

–  –  –

Граничные условия задаются в каждый момент времени так же, как в стационарной задаче .

Физических задач, описываемых уравнением (4.48), достаточно много. Например, при = 0 это уравнение становится параболическим и является обычным уравнением нестационарной теплопроводности или задачей консолидации грунта как разновидность задач нестационарной фильтрации и т.д .

При µ = 0 уравнение (4.48) превращается в волновое уравнение, описывающее, например, электромагнитные волны, поверхностные волны в жидкости, волны расширения-сжатия и т.д .

При 0 и µ 0 (4.48) является волновым уравнением с демпфированием и используется, например, для описания волновых явлений в механике жидкости и газа [8] .

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

Перемещения упругого тела во времени обусловлено наличием двух систем дополнительных сил. Первую из них составляют силы инерции, которые характеризуют ускорение и согласно принципу Даламбера могут быть записаны в виде {u} .

(4.59) t 2 Здесь {u} – обобщенные перемещения. Эти силы совпадают по направлению с перемещениями {u} и обычно отнесены к единице объема, а – плотность (масса единицы объема) .

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

в общем случае они связаны нелинейной зависимостью со скоростью .

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

µ (4.60) t Здесь µ – некоторый коэффициент .

Тогда получаем эквивалентную квазистатическую задачу, причем распределенные объемные нагрузки определяются соотношением

–  –  –

где [K] и { F } – глобальные матрица жесткости и вектор нагрузок, полученные ансамблированием соответствующих локальных величин. Матрицы [C] и [M] получаются ансамблированием локальных матриц, задаваемых в виде <

–  –  –

Здесь е – площадь треугольного элемента. Если толщина t элемента предполагается постоянной в пределах элемента, для локальной матрицы масс из (4.642) имеем

–  –  –

. .

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

Выше различные нестационарные и динамические задачи были сведены к системам обыкновенных дифференциальных уравнений (задача Коши) с использованием конечно-элементой аппроксимации только для пространственных переменных. Для однозначного решения таких систем необходимо задать соответствующие начальные условия (значения в начальный момент времени на саму функцию для параболических уравнений или значения самой функции и ее скорости в начальный момент времени для гиперболических уравнений). Для решения этих систем можно применять стандартные методы решения систем ОДУ (метод Эйлера, Рунге – Кутты и т.д.) или использовать конечно-разностные аппроксимации для производных по времени [6, 11, 12] .

4.5.2. Разрешающие соотношения МКЭ для пространственной и временной конечно-элементной аппроксимации Возможны и другие подходы к решению нестационарных и динамических задач [8] .

Во-первых, к дифференциальному уравнению, описывающему задачу, может быть непосредственно применена процедура Галеркина (или другой проекционный метод), но функции формы будут зависеть не только от пространственных переменных, но и от времени:

–  –  –

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

Во-вторых, можно применить вариационный принцип по пространственным переменным и времени. В частности, для простых динамических задач вариационный принцип непосредственно следует из принципа Лагранжа. Ищется стационарное значение интеграла t2 = L dt, L = U + W + T, t1 где U, W, T – энергия деформации, потенциальная и кинетическая энергии системы соответственно. Здесь также строятся пространственные и временные конечные элементы .

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

Задачи, описываемые дифференциальными уравнениями первого порядка по времени. Типичной для этого класса является задача, описываемая уравнением [8] [ K ]{} + [C ] {} + {F } = 0. (4.65) t

–  –  –

Можно рассмотреть квадратичную трехточечную аппроксимацию производной по времени на интервале времени t [0, 2t].

При этом аппроксимируем функции с помощью квадратичных лагранжевых элементов (4.2):

<

–  –  –

Для решения этого уравнения необходимо задать два начальных условия: на саму функцию и ее скорость. Поскольку в каждый момент необходимо определять в узлах саму функцию и ее производную по времени, то для интерполяции функции на интервале времени t [0, t] выбираются эрмитовы конечные элементы. Таким образом имеем

–  –  –

Однослойную временную схему можно получить, если взвесить невязку дважды по Н01 и по Н11 (необходимо найти значение функции и ее скорость на конец временного интервала, на начало интервала они должны быть известны):

–  –  –

Очевидно, можно применить и более сложные конечные временные элементы с дополнительными степенями свободы .

Одним из классов динамических задач являются задачи на колебания и собственные значения. Здесь мы не будем останавливаться на особенностях их решения. За подробностями можно обратиться к [8] .

Задачи:

1. Получить аппроксимацию уравнения с первой производной по времени с помощью линейного конечного элемента (соотношение (4.70)) .

2. Вывести разностную схему (аналогичную (4.70)) для аппроксимации уравнения с первой производной по времени с использованием «временных» КЭ 2-го порядка лагранжева типа .

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

Вопросы для самопроверки

1. Приведите примеры нестационарных процессов. Запишите постановки соответствующих краевых задач .

2. Показать что дифференциальное уравнение (4.49) является уравнением Эйлера для функционала (4.50) .

3. Получить локальные разрешающие соотношения (4.57) из условия минимума функционала (4.50) для линейного треугольного элемента .

4. Получить локальные разрешающие соотношения (4.57) из условия минимума функционала (4.50) для билинейного четырехугольного элемента .

5. Вывести локальные разрешающие соотношения (4.58) для нестационарной задачи в компонентной форме .

6. Объясните суть временной конечно-элементной аппроксимации (аппроксимации первой производной по времени с помощью конечных элементов) .

4.6. ФИЗИЧЕСКИ НЕЛИНЕЙНЫЕ ЗАДАЧИ. ПЛАСТИЧНОСТЬ, ПОЛЗУЧЕСТЬ

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

а) линейной связи между деформациями и перемещениями (геометрическая линейность);

б) линейной связи между напряжениями и деформациями (физическая линейность) .

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

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

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

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

Такие задачи будут рассмотрены в следующем разделе .

Следует отметить существенный факт. В нелинейных задачах, в отличие от линейных, решение часто является неединственным .

При выводе СЛАУ для линейной теории упругости (4.53) использовался линейный закон Гука в виде {} = [ D ]({} { 0 }) + { 0 }. (4.77) Кроме того, предполагалась линейная связь между деформациями и перемещениями, перемещения считались непрерывными, и уравнения равновесия удовлетворялись приближенно .

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

F ({},{}) = 0. (4.78)

Если удастся найти такое решение (4.47), что при соответствующем подборе входящих в (4.77) параметров [D], {0} или {0} это уравнение и соотношение (4.78) удовлетворяются при одинаковых значениях напряжений и деформаций, то полученное решение будет искомым .

При решении таких задач целесообразно использовать итерационные методы .

Если при итерациях подбирается матрица [D], то приходим к методу переменных параметров упругости (методу переменной жесткости) [6] .

Если же подбираются {0} или {0}, то это так называемые методы начальных деформаций или начальных напряжений .

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

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

4.6.1. Метод переменных параметров упругости

Метод переменных параметров упругости можно использовать в случае, если связь между напряжениями и деформациями (4.73) можно представить в форме (4.72), где матрица упругих констант зависит от достигнутого уровня деформации:

–  –  –

которое можно решить различными итерационными методами [11, 12] .

Например, метод простых итераций. Сначала предполагается, что {u} = 0, вычисляется [K({u}(0))] = [K(0)] и определяется {u}(1) = [ K (0) ]1{F } .

(0) Процесс повторяется в соответствии с формулой

–  –  –

до тех пор, пока перемещения перестанут изменяться .

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

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

Итак, деформационная теория пластичности предполагает наличие однозначной зависимости между суммарными деформациями и напряжениями в упругопластическом теле, которые для изотропно тела имеют вид

–  –  –

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

Определяют поле перемещений, по нему находят поле деформаций (и интенсивность деформации) и интенсивность соответствующих упругих * напряжений: i * E i. Затем вычисляют = i, переопределяi ют Е* и * и решают упругую задачу с новыми эффективными значениями параметров упругости (рис. 4.23). Расчет заканчивается при достаточной близости двух соседних приближений для напряжений .

Рис. 4.23. Схема расчета по методу переменных параметров упругости [15] Учет ползучести по теории течения или упрочнения проводят по этапам времени. На начальном этапе решается упругая задача с эффективными параметрами упругости с определяющими соотношениями (4.82). Для следующего момента времени t определяющие соотношения принимают вид

–  –  –

где ic* – накопленная к моменту времени t деформация ползучести .

Функция текучести определяется для значений параметров, соответствующих началу этапа нагружения (i(0), T(0), t0 = 0, ic* (0) = 0), величины ij(0) представляют собой напряжения в упругом материале при отсутствии деформаций ползучести .

Для следующего этапа нагружения t2 уравнения имеют вид

–  –  –

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

Наиболее просто изложенный метод реализуется при отсутствии пластических деформаций, и деформации ползучести развиваются в упругом теле. Тогда на всех этапах полагают Е* = Е и * =. Если уже на первом этапе нагружения t = 0 возникают пластические деформации (напряжения превосходят предел текучести материала), то для расчета используют метод переменных параметров упругости. Этот метод применяют и на втором этапе нагружения, причем расчет считается достоверным, если в точках, в которых имелась пластическая деформация в конце первого этапа нагружения, пластическая деформация возрастает (увеличивается интенсивность напряжений). Если это условие не выполняется, то расчет проводят снова, причем в точках разгрузки принимают Е* = Е и * =. Аналогично выполняется расчет последующих этапов нагружения .

Метод переменных параметров упругости при сложном нагружении. В [15] рассматривается алгоритм применения метода переменных параметров упругости для задач пластичности и ползучести при сложном нагружении. В этом случае нагружение разбивается на ряд этапов. В большинстве случаев это разбиение целесообразно проводить по этапам действительной истории нагружения во времени. На каждом этапе нагружения должны выполняться уравнения равновесия в приращениях k + k f = 0

–  –  –

Таким образом, упругопластическую задачу представили в матричной форме упругой задачи, где параметры упругости зависят от напряженного состояния .

Интегрируя по времени, получим приращение полной деформации для k-го этапа нагружения

–  –  –

т.е. принимают параметры упругости, соответствующие напряженному и деформированному состоянию в конце предыдущего этапа нагружения .

Принимают сначала, что характер нагружения (нагрузка или разгрузка) остается таким же, как и на предшествующем этапе. Аналогичное предположение принимают для вектора

–  –  –

Далее проводят расчет k-го этапа нагружения упругого тела с эффективными упругими коэффициентами и заданными дополнительными деформациями ползучести. После расчета проверяют условия возникновения пластической деформации:

–  –  –

где i(k) – интенсивность напряжений в конце k-го этапа нагружения;

т ( ip k ) ) – предел текучести, соответствующий накопленной интенсивности пластической деформации в конце k-го этапа нагружения .

Если характер нагружения на k-м этапе остается таким же, как и на предыдущем (k – 1)-м, то корректировка значения [a]k–1 не требуется. Если на k-м этапе предполагали нагружение, а после проверки критерия пластичности он оказался невыполненным, то расчет этапа следует провести заново, положив

–  –  –

Если на k-м этапе предполагали разгрузку (т.е. решалась чисто упругая задача без пластики), но в действительности было нагружение, то расчет следует провести заново, принимая

–  –  –

В результате рассмотренной процедуры в первом приближении устанавливают приращение напряжений и деформаций на k-м этапе и значение [a]k .

Если принятые в начале этапа величины [a]k–1 и полученные после расчета [a]k достаточно близки, то расчет этапа заканчивают и приступают к следующему шагу. Если расхождения параметров упругости велики, полагают ([a]k 1 + [a]k ) [a ] = и проводят расчет второго приближения .

При проведении расчетов подобным образом уточняют значения вектора деформаций текучести {c} .

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

4.6.2. Методы начальных (дополнительных) напряжений Если определяющие уравнения разрешены относительно напряжений = f (), то их можно свести к соотношениям типа (4.77) для упругого материала, задавая соответствующим образом {0}. Tак как {0} влияет на силы {F }, то приходим к уравнению

–  –  –

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

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

В этом методе на каждом итерационном шаге используется одна и та же матрица жесткости .

Для определения матрицы жесткости [K(0)] следует использовать начальные значения упругих констант, если поведение материала в основном описывается соотношениями линейной теории упругости и отклонения от него локализованы. Однако если нелинейность появляется для всех напряжений, то для ускорения сходимости можно скорректировать упругие постоянные после первой итерации .

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

= f (). (4.86) Приведение соотношения (4.86) к виду (4.77) может быть достигнуто соответствующим выбором {0}. Уравнение (4.78) опять решается итерационным методом, но теперь упругие деформации, получаемые на каждом шаге, сравниваются с деформациями, соответствующими определяющему соотношению (4.86), и их разность используется для оценки невязки силы {F ( n ) }. В остальном процесс идентичен описанному выше. Матрица жесткости остается постоянной на любом шаге .

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

Метод дополнительных деформаций при простом нагружении .

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

Определяющие соотношении запишем в виде

–  –  –

Далее находим дополнительные деформации p (1) = 1 e(1) .

Во втором приближении рассматривается та же упругая задача, но при наличии дополнительных деформаций (2) = I + 2G ( (2) p (1) ) .

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

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

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

Метод дополнительных деформаций при сложном нагружении .

В [15] рассматривается алгоритм применения метода переменных параметров упругости для задач пластичности и ползучести при сложном нагружении .

Вектор приращения полных деформаций представляют в форме

–  –  –

ложение изображающей точки Аk–1 на кривой деформирования. Если точка Аk–1 находится внутри отрезка Оk–1Аk–1 (например, А’k–1), то первое приближение завершает расчет. Рассмотрим случай, когда точка Аk–1 расположена на кривой деформирования и i(k–1) = т (k–1) (рис. 4.25) .

При расчете первого приближения значения [a]e, {c} принимают соответствующими точке Аk–1 .

После проведения расчета проверяют выполнение условия нагружения i ( k ) = т ( ip k ) ). Если имеет место нагружение, то определяют вели

–  –  –

где Eк( k 1) – касательный модуль в точке Аk–1, проводят расчет второго приближения, в котором дополнительную пластическую деформацию считают равной F (i ){S} (1)i .

k В расчете принимают среднее значение функции F (i ) по ее значениям для точек Аk–1 и Ak(1) .

Подобным образом проводят расчет следующего приближения. Расчет заканчивается при достаточной близости двух соседних приближений. В результате расчета находят приращения напряжений и деформаций на k-м этапе нагружения, векторы напряжений и деформаций в конце этапа нагружения, значения i(k), ip k ). Расчет следующего (k + 1)-го этапа *( проводят по описанной процедуре .

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

–  –  –

напряжения {0}(1), тогда как в методе начальных деформаций значения деформаций корректируются поправочным членом {0}(1). Отметим, что когда с ростом напряжений деформации быстро увеличиваются, предпочтительнее использовать первый метод, а когда с ростом деформации быстро увеличиваются напряжения – второй .

4.6.4. Метод Ньютона для решения нелинейной системы уравнений, полученной для физически нелинейной задачи Заметим, что систему уравнений (4.80) можно рассматривать как систему нелинейных уравнений (х) = 0, для которой применимы все известные методы решения нелинейных систем уравнений и, в частности, метод

Ньютона (метод Ньютона–Рафсона) [7]:

–  –  –

Очевидно, что методы переменных параметров упругости, начальных напряжений и деформаций относятся к этим двум категориям [8] .

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

Используя принцип виртуальной работы, вновь запишем равенство изменений внешней и внутренней работ:

–  –  –

где вектор {F } содержит все внешние силы, обусловленные приложенными нагрузками. Если для вариации деформации справедливо соотношение

–  –  –

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

–  –  –

{u ( n+1) } = [ K т ( n ) ]1 {}( n ), (4.87) где [Kт(n)] – матрица касательных упругих постоянных, определенная для перемещений и деформаций, соответствующих приближенному решению {u(n)}. Заметим, что формула справедлива, если для любого приближения {u(n)} матрица [Kт(n)] имеет обратную .

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

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

Вопросы для самопроверки

1. Что такое физически нелинейные задачи? Приведите примеры постановок в МДТТ таких задач .

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

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

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

5. В чем суть метода дополнительных нагрузок для решения физически нелинейных задач?

6. В чем суть метода дополнительных деформаций для решения физически нелинейных задач?

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

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

9. В чем принципиальное отличие применения метода дополнительных деформаций при простом и сложном нагружении?

10. Сформулируйте идею применения метода Ньютона к решению системы нелинейных алгебраических уравнений, полученную в МКЭ для физически нелинейной задачи. Что такое матрица касательных упругих постоянных?

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

4.7. ГЕОМЕТРИЧЕСКИ НЕЛИНЕЙНЫЕ ЗАДАЧИ

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

Независимо от величины перемещений и деформаций внутренние и внешние силы должны удовлетворять условиям равновесия, которые вытекают из принципа виртуальной работы (или преобразованной по формуле Грина взвешенной невязки, полученной подстановкой приближенного значения напряжения в уравнение равновесия) [8]:

–  –  –

Черта сверху означает, что при больших перемещениях деформации нелинейно зависят от перемещений и матрица [ B ] зависит от {u}. Эту матрицу удобно представить виде

–  –  –

где [B0] – матрица, определяющая бесконечно малые деформации, а матрица [BL] зависит от перемещений. В общем случае [BL] является линейной функцией перемещений .

Очевидно, что уравнение (4.88) следует решать итерационными методами .

При использовании метода Ньютона необходимо найти зависимость между d{u} и d{}. Вариация (4.88) по d{u} имеет вид

–  –  –

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

Первое слагаемое в (4.92) можно записать в виде

–  –  –

где [K] – симметричная матрица, зависящая от величины напряжения .

Эта матрица называется матрицей начальных напряжений или геометрической матрицей. Таким образом,

–  –  –

где [K] – полная матрица тангенциальных жесткостей .

Итерационная процедура метода Ньютона строится следующим образом:

а) в качестве первого приближения {u} строится решение линейной теории упругости;

б) с помощью соотношения (4.88) определяется {(1)} для заданной матрицы [ B ] и напряжений, определяемых равенством {} = [ D ]({} { 0 }) + +{ 0 } (или другим определяющим соотношением);

в) строится матрица [K];

г) определяется поправка d{u (1) } = [ K ]1{(1) } .

Процесс повторяется до тех пор, пока величина d {u( n ) } не станет достаточно малой .

Все решения можно находить за один шаг для полной действующей нагрузки. Однако, как и во всех нелинейных задачах, возникает возможность неединственности решения и при этом может быть найдено решение, не имеющее физического смысла. В таких случаях целесообразно задавать нагрузки отдельными приращениями и получать нелинейное решение для каждого приращения. С вычислительной точки зрения это часто экономичнее, поскольку эффекты нелинейности на каждом шаге становятся меньше. Если приращения нагрузки достаточно малы по величине, то каждое решение в приращениях с достаточной степенью точности может быть найдено за один шаг (то есть описанный метод эквивалентен методу Эйлера; его можно уточнить, применяя методы Рунге–Кутты). Однако необходимо периодически проверять выполнение условия равновесия с помощью нелинейного соотношений (4.88) .

4.7.2. Общий случай больших деформаций и напряжений Рассмотрим случай больших градиентов перемещений и деформаций, которые описываются тензором Грина[8] = ( u + u + u u ), компоненты которого в декартовой ортогональной системе координат

–  –  –

где [M] – матрица размерностью 99 из шести компонент вектора напряжений, расставленных как показано в (4.104). Очевидно, что матрица [K] – симметрична .

Заметим, что все предыдущие матричные соотношения получены для конечного элемента (локальные), индекс элемента опущен для простоты .

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

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

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

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

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

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

Задачи:

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

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

3. Используя свойства, полученные в предыдущих задачах 1, 2, вывести соотношение (4.109) .

Вопросы для самопроверки

1. Какие задачи в МДТТ называются геометрически нелинейными?

Приведите примеры постановок таких задач .

2. Как привести матричную постановку геометрически нелинейной задачи к виду итерационного метода Ньютона для решения систем нелинейных уравнений?

3. Что собой представляют матрица начальных перемещений (матрица больших перемещений), матрица начальных напряжений (геометрическая матрица), полная матрица тангенциальных жесткостей?

4. Опишите алгоритм метода Ньютона для геометрически нелинейной задачи .

ВЫВОДЫ ПО ГЛАВЕ 4

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

Популярность МКЭ связана с рядом преимуществ этого метода:

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

простота исследования неоднородных и анизотропных тел (в частности, когда направление анизотропии переменное);

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

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

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

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

МКЭ позволяет рассматривать граничные условия с разрывной поверхностной нагрузкой, а также смешанные граничные условия .

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

Хотя ресурсы совершенствования МКЭ практически исчерпаны, однако ведется разработка численных методов, а также реализующих их программных комплексов, позволяющих более экономично использовать вычислительные ресурсы и гарантировать эффективное решение многовариантных задач анализа и проектирования. Например, создан комбинированный метод конечных и граничных элементов (КМКиГЭ), реализующий достоинства МКЭ и не имеющий его недостатков [17] .

Список литературы к главе 4

1. Каримов И. Электронный учебный курс для студентов очной и заочной формы обучения [Электронный ресурс]. – URL: http: //www.stroitmeh .

ru/lect31.htm

2. Сагдеева Ю.А., Копысов С.П., Новиков А.К. Введение в метод конечных элементов: метод. пособие. – Ижевск: Изд-во Удмурт. ун-та, 2011. – 44 с .

3. Зенкевич О., Морган К. Конечные элементы и аппроксимация. – М.: Мир, 1986. – 318 с .

4. Галлагер Р. Метод конечных элементов. Основы: пер. с англ. – М.:

Мир, 1984. – 428 с .

5. Норри Д., де Фриз Ж. Введение в метод конечных элементов. – М.:

Мир, 1981. – 304 с .

6. Бояршинов М.Г. Численные методы: учеб. пособие для студентов направления «Прикладная математика». Ч. 3 / Перм. гос. техн. ун-т. – Пермь, 2002. – 134 с .

7. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. – М.: Наука, 1977. – 832 с .

8. Зенкевич О. Метод конечных элементов в технике: пер. с англ. – М.: Мир, 1975. – 543 c .

9. Сегерлинд Л. Применение метода конечных элементов. – М.: Мир, 1979. – 392 с .

10. Терпугов В.Н., Лалин В.В. Конечно-элементные технологии построения расчетных алгоритмов для решения задач механики сплошных сред: метод. пособие. – Пермь: Изд-во Перм. гос. ун-та, 2008. – 380 с .

11. Бояршинов М.Г. Численные методы: учеб. пособие для студентов направления «Прикладная математика». Ч. 2 / Перм. гос. техн. ун-т. – Пермь, 1999. – 200 с .

12. Самарский А.А., Гулин А.В. Численные методы: учеб. пособие для вузов. – М.: Наука, 1989. – 432 с .

13. Оден Дж. Конечные элементы в нелинейной механике сплошных сред: пер. с англ. – М.: Мир, 1976. – 464 с .

14. Клованич С.Ф. МКЭ в нелинейных задачах инженерной механики. – Запорожье: Изд-во журнала «Свiт геотехнiки». – 2009. – 400 с .

15. Биргер И.А., Шорр Б.Ф., Иосилевич Г.Б. Расчет на прочность деталей машин: справочник. – М.: Машиностроение, 1993. – 640 с .

16. Коннор Дж., Бребия К. Метод конечных элементов в механике жидкости. – Л.: Судостроение, 1979. – 264 с .

17. Борисов А.В. Численное моделирование физических процессов с применением метода конечных элементов на базе COMSOL Multiphysics[Электронный ресурс]. – URL: http: //edu2.tsu.ru/html/3024_new/

ГЛАВА 5. ОБРАТНЫЕ И НЕКОРРЕКТНЫЕ ЗАДАЧИ

5.1. ОСНОВНЫЕ ПОНЯТИЯ И ПРИМЕРЫ Работы по обратным и некорректным задачам появились в первой половине XX века. Они были связаны с физикой (обратные задачи квантовой теории рассеяния), геофизикой (обратные задачи электроразведки, сейсмики, теории потенциала), астрономией и другими областями естествознания. С появлением ЭВМ область приложения обратных и некорректных задач охватила практически все научные дисциплины, в которых используются математические методы. В прямых задачах математической физики необходимо определить поля физических величин (температуру, напряжения, скорости и т.д.), при этом свойства исследуемой среды (коэффициенты уравнений), а также начальное состояние и условия на границе предполагаются известными. Однако на практике именно свойства среды часто являются неизвестными, встает проблема надежной идентификации материалов, т.е. определения физических характеристик, используемых для инженерных расчетов. И тогда возникают обратные задачи, в которых по информации о решении прямой задачи требуется определить коэффициенты уравнений. Эти задачи в большинстве случаев некорректны (неустойчивы по отношению к погрешностям измерений или имеющие неединственное решение) [1]. Решение обратных задач может помочь определить местоположение, форму и структуру включений, дефектов, источников (тепла, колебаний, напряжения, загрязнения и т.д.). Теория обратных и некорректных задач в настоящее время является одной из наиболее стремительно развивающихся областей науки. Обратные задачи имеют постоянно расширяющиеся области применения в инженерной практике, среди которых выделим следующие:

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

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

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

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

восстановление формы трехмерного тела по набору его ортогональных проекций на различные плоскости (рентгеновская томография, ЯМР-томография, УЗИ и т.д.);

экология (диагностика состояния воздуха, воды, космический мониторинг и т.д.);

экономика (теория оптимального управления, финансовая математика и т.д.) .

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

1) область определения независимых параметров (координат, времени);

2) уравнение, описывающее данный процесс;

3) начальные условия (для нестационарного процесса);

4) условия на границе исследуемой области .

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

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

Обратная задача называется ретроспективной, если требуется определить начальные условия .

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

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

Обратная задача называется задачей об источнике, если требуется определить источник .

Обратная задача называется коэффициентной, если требуется восстановить коэффициенты, входящие в основные уравнения .

В 1902 году Ж. Адамар сформулировал понятие корректности поставленной задачи для дифференциальных уравнений. Задача поставлена корректно в классическом смысле (по Ж. Адамару), если [2]

1) решение задачи существует;

2) это решение единственно;

3) решение задачи непрерывно зависит от входных данных .

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

Каждую некорректную задачу можно сформулировать как обратную к некоторой корректной прямой задаче .

В [1] приведены примеры корректных и, соответствующих им, некорректных задач:

–  –  –

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

Пример Тихонова [3]. Линейная алгебраическая система с плохо обусловленной (вырожденной) матрицей

–  –  –

Матрица системы и расширенная матрица имеют ранг 1. Таким образом, система совместна и имеет множество решений, зависящих от одного параметра. Возьмем последовательность 100, 200, 300 десятичных знаков и так далее в приближении иррациональных чисел во втором уравнении. Тогда определитель системы не равен нулю и можно построить единственное решение хп, уп, соответствующее сохранению п значащих цифр в представлении иррациональных чисел. Оказывается, что построенная таким образом последовательность решений не имеет предела при п, хотя исходная информация становится все точнее и точнее. Таким образом, задача решения такой СЛАУ с вырожденной матрицей будет некорректной .

В более общем случае, если определитель равен нулю, то система имеет единственное (тривиальные) решение только при нулевой правой части. Иначе решение такой системы неединственно или вообще не существует, т.е. задача решения СЛАУ с вырожденной матрицей или близкой к ней является некорректной. Эта ситуации типична, если матрица находится в результате некоторых приближенных вычислений и если в результате расчетов получено, например, что det A = 10–30. В этом случае ситуации, когда определитель равен нулю и решение существует только при определенных ограничениях на правую часть либо не равен нулю и система имеет единственное решение, в принципе неразличимы .

Пример Адамара [3].

Рассмотрим следующую задачу Коши для уравнения Лапласа:

utt ( x, t ) = u xx ( x, t ), t 0, u ( x,0) = 0, ut ( x,0) = sin k x .

k Несложно показать, что решением этой задачи будет функция

–  –  –

при k. Таким образом, возмущения в «начальном» условии, сколь бы малыми они ни были, неограниченно возрастают при t T .

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

Будем говорить, что задача поставлена корректно по Тихонову, если [2]:

1) априори известно, что решение задачи существует в некотором классе;

2) в этом классе решение единственно;

3) решение задачи зависит непрерывно от входных данных .

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

В случае если задача A u = f условно корректна на некотором множестве M и это множество – компакт, возможно построение оценок устойчивости нахождения обратного оператора .

Теорема. Пусть оператор А отображает компакт M непрерывно и взаимно однозначно на множество N. Тогда существует непрерывная функция (), (0) = 0, такая, что

–  –  –

Таким образом, A 0. При выводе учтено, что u(x) и w(x) обращаются в нуль на границе .

Простейшую априорную оценку для данной граничной задачи с начальными условиями получим на основе Леммы Гронуолла [2] .

Лемма. Для функции g(t), удовлетворяющей неравенству

–  –  –

Это значит, что для задачи с обратным временем имеет место непрерывная зависимость решения от начальных данных при 0 t T в классе ограниченных решений .

Задачи:

1. Покажите самосопряженность и неотрицательность дифференциального оператора одномерного параболического уравнения .

2. Докажите лемму Гронуолла .

3. Докажите теорему об априорной оценке решения задачи с обратным временем .

4. Докажите единственность решения задачи с обратным временем .

5. Докажите устойчивость решения задачи с обратным временем .

Вопросы для самопроверки

1. В каких практических проблемах возникают обратные задачи?

Поясните суть обратных задач .

2. Перечислите основные типы обратных задач. В каких задачах практики они возникают? Приведите примеры .

3. Сформулируйте определение корректной (по Адамару) задачи .

4. Что подразумевается под понятием некорректной задачи?

5. Какая связь между обратными и некорректными задачами?

6. Приведите примеры некорректных задач .

7. Объясните суть некорректности в примере Тихонова .

8. Объясните суть некорректности в примере Адамара .

9. Объясните суть некорректности в задаче с обратным временем .

10. Сформулируйте определение условно корректной (по Тихонову) задачи .

5.2. СПОСОБЫ ПРЕОДОЛЕНИЯ НЕКОРРЕКТНОСТИ .

МЕТОДЫ РЕГУЛЯРИЗАЦИИ

Определение. Оператор R ( f, ) : F U называется регуляризующим для уравнения Au = f, u U, f F или регуляризатором, если он обладает следующими свойствами [3]:

1) 1 0 : R( f, ) определен для [0, 1 ] и f F, удовлетворяющих неравенству

–  –  –

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

5.2.1. Метод квазирешений Рассмотрим некорректную задачу Au = f. Предположим, что для точной правой части f существует единственное решение, принадлежащее некоторому компакту M. Однако элемент f неизвестен, заданы его приближение f и величина погрешности, такие, что f f. Требуется, зная эту информацию, построить приближенное решение и, которое стремится к точному u при 0 [3]. Введем в рассмотрение множество U = {u U : Au f } .

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

–  –  –

Введем в рассмотрение понятие квазирешения. Пусть в задаче Au = f оператор А: U F непрерывен, а M – компакт в U .

Определение. Квазирешением уравнения Au = f называется элемент иk, минимизирующий невязку

–  –  –

на множестве M .

Из определения квазирешения следует, что оно существует f F и в силу непрерывности оператора А функционал J достигает на компакте M своей точной нижней грани; в случае когда f А (M), иk совпадает с обычным решением, поскольку min J [иk] = 0 .

Имеет место следующая теорема .

Теорема. Пусть А: U F – линейный непрерывный оператор, такой, что уравнение Au = 0 имеет единственное решение на M (M – выпуклый компакт в U); тогда для любой функции f F квазирешение существует, единственно и непрерывно зависит от f .

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

–  –  –

Необходимо найти приближенное решение уравнения (5.1) с приближенно заданной правой частью f. Это приближенное решение обозначим u, где параметр = () связан с величиной погрешности в задании правой части .

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

–  –  –

в которой оператор A обладает лучшими свойствами, чем А .

В вариационных методах (в методе квазирешений) вместо решения уравнения (5.3) минимизируется норма невязки r = Av – f, т.е. минимизируется функционал невязки

–  –  –

В (5.4) 0 – параметр регуляризации, величина которого согласуется с погрешностью задания правой части. Для выделения ограниченного решения в функционал невязки добавлен стабилизирующий функционал v .

Сходимость метода регуляризации Тихонова. Для доказательства сходимости данного алгоритма представим приближенное решение в операторном виде

u = R () f, (5.6)

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

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

–  –  –

Близость R () f к u установим для класса ограниченных функций (5.7). Покажем, что для любого 0 можно указать () такое, что z2 для всех 0 (). Для функций из класса (5.7) ряд

–  –  –

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

Если оператор А самосопряженный, положительный, то можно ограничиться возмущением самого оператора:

Au + u = f. (5.15)

Задача (5.15) соответствует использованию алгоритма упрощенной регуляризации .

Фактически, помимо вариационных методов решения некорректных задач типа (5.4), (5.5), можно выделить второй класс приближенных методов (5.14), (5.15), которые характеризуются возмущением оператора исходной или преобразованной задачи .

Вместо (5.7) потребуем выполнения более сильных ограничений на точное решение (5.1):

–  –  –

В случае самосопряженного и положительного оператора А промежуточное между (5.7) и (5.16) положение занимает класс априорных ограничений на решение типа

–  –  –

В условиях (5.16), (5.17) конкретизируем зависимость s () в оценке погрешности типа (5.13) для метода регуляризации .

Для A* = A 0 соотношение (5.14) примет вид

–  –  –

получим 21 1/2 2 z f f + A u .

С учетом (5.2), (5.17) имеем оценку (5.20) .

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

Выбор параметра регуляризации. Рассмотрим некоторые способы выбора параметра регуляризации, который согласуется с погрешностью входных данных: = (). Наибольшее распространение получили: выбор параметра регуляризации по невязке, обобщенной невязке (при учете погрешности в задании не только правой части, но и оператора А), квазиоптимальный выбор и т.д. Выбор оптимального параметра во многом определяет работоспособность вычислительного алгоритма .

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

Кроме того, оптимальный параметр регуляризации зависит от класса априорных ограничений на точное решение. Например, для ограниченных решений (класс (5.7)) полученная оценка (5.13) для погрешности приближенного решения не позволяет указать оптимальное значение параметра регуляризации .

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

opt =, (5.22) M

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

–  –  –

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

Метод невязки. При выборе параметра регуляризации по невязке в качестве определяющего уравнения выступает равенство

–  –  –

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

Обозначим

–  –  –

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

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

–  –  –

и вычисления проводятся, начиная с k = 0 до некоторого k = K, при котором равенство (5.26) выполняется с приемлемой точностью. При таком определении параметра регуляризации требуется К + 1 вычислений невязки (решений вариационных задач типа (5.4), (5.5) или уравнений Эйлера (5.13)) .

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

Поэтому для решения уравнения () = можно использовать итерационный метод Ньютона [4]:

–  –  –

который сходится при любом начальном приближении 0 0. Можно применять модификацию метода Ньютона (метод секущих), не требующую вычисления производной функции () [4]:

–  –  –

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

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

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

–  –  –

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

5.2.3. Метод регуляризации на компактных множествах Вновь рассмотрим операторное уравнение Au = f в случае некорректной задачи. Из метода квазирешений известно, что для процедуры построения квазирешения необходимо выбирать компактное множество M .

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

Критерий Арцела [5]. Множество M C[a, b] компактно, если функции из M равномерно ограничены и равностепенно непрерывны .

Рассмотрим несколько типичных способов выбора компакта M .

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

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

1) Пусть Mc – множество ограниченных монотонно невозрастающих функций:

–  –  –

Тогда справедлива теорема: Mc компактно в L2 [a, b] .

Отметим, что определенные ниже множества также являются компактами в L2 [a, b] .

)

2) M c – множество ограниченных выпуклых вверх функций:

–  –  –

и представляет собой компактное множество .

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

–  –  –

Здесь В, k+1, k = 0, 1, … – итерационные параметры (различные для различных итерационных методов), влияющие на скорость сходимости .

Итерационный метод может интерпретироваться как итерационный метод решения вариационной задачи минимизации невязки

–  –  –

В случае задачи (5.3) для принятых допущений о свойствах оператора А min = 0, max = 1, т.е. = 0. Поэтому невозможно конкретизировать скорость сходимости итерационного метода .

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

Итерационное решение некорректной задачи. Следующая теорема формулирует условия, при которых метод простых итераций дает приближенное решение задачи (5.1), (5.2) .

Теорема 4. Пусть в итерационном методе (5 .

31) – (5.33) число итераций n () и n () 0 при 0. Тогда un () u 0, если 0 .

Доказательство .

Обозначим через zn = un – u погрешность на n-й итерации. Из (5.31) непосредственно получаем n 1 un = ( E A)n u0 + ( E A) k f, (5.34) k =0 где u0 – некоторое начальное приближение .

Для точного решения можно воспользоваться аналогичным представлением

–  –  –

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

С учетом (5.34) для погрешности получаем

–  –  –

причем z0 = u0 – u – начальная погрешность. Первое слагаемое в (5.35) является стандартным для итерационных методов, второе – связано с учетом погрешности в задании правой части уравнения (5.1) .

При сформулированных ограничениях (5.33) на итерационный параметр имеем

–  –  –

где s (n) 0 при n. Из полученной оценки (5.38) непосредственно следует доказываемое утверждение .

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

Оценка скорости сходимости.

Для оценки скорости сходимости метода простых итераций потребуем более жесткие, чем (5.33), ограничения на итерационный параметр :

0. (5.39) Теорема 5 .

Пусть точное решение задачи (5.3) принадлежит классу

–  –  –

Согласно оценке (5.43) скорость сходимости приближенного решения находится в прямой зависимости от погрешности в задании правой части и гладкости точного решения (параметр р) Обобщения. Выше отмечалась возможность оптимального выбора числа итераций в случае известного класса априорных ограничений на точное решение задачи (постоянная М в (5.40)). При практическом использовании итерационных методов решения некорректных задач класс априорных ограничений часто не удается сформулировать явно. Поэтому вместо (5.42) оптимальное число итераций выбирается по невязке, т.е. итерационный процесс продолжается до тех пор, пока не будет выполнено неравенство

–  –  –

Для решения СЛАУ рассматривались итерационные методы вариационного типа (применимые и для решения других классов задач). Они характеризуются использованием явных расчетных формул для итерационных параметров. Например, метод скорейшего спуска (наименьших невязок) задается итерационной формулой

–  –  –

Не останавливаясь на регуляризирующих свойствах таких итерационных методов, отметим только, что для вариационных итерационных методов (5.46), (5.47) и (5.48), (5.49) для решения некорректной задачи (5.1), (5.2) получены аналогичные результаты, что и для метода простых итераций (5.31) .

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

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

5.2.5. Метод усеченных сингулярных разложений Метод усеченных сингулярных разложений связан с ортогональными разложениями, в частности с рядами Фурье [3] .

Определение. Сингулярным разложением оператора А: U F называется его представление в виде Au = k (u, uk ) f k, (5.50) k =1

–  –  –

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

Анализируя структуру сингулярного разложения, можно понять природу некорректности. Пусть f аппроксимирует f, т.е. f f F. Если известно только f, то для коэффициентов разложения R f из теоремы можно получить оценку 1 ( f, f k ) 1 ( f, f k ) 1 .

k k k Из этой оценки ясно, что вклад слагаемых uk в разложении A+f для малых k (практически k ) нельзя вычислить с приемлемой точностью .

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

Отметим, что число N слагаемых, участвующих в представлении (5.54), определяется из критерия невязки:

–  –  –

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

5.2.6. Проекционный метод Одним из эффективных способов построения приближенных решений операторных уравнений является проекционный метод, согласно которому приближенное решение операторного уравнения ищется в виде линейной комбинации некоторой линейной независимой системы функций – базисных функций [3]. Таковыми являются хорошо известные методы взвешенных невязок Бубнова–Галеркина, метод моментов (Петрова–Галеркина) и т.д., широко применяемые для решения дифференциальных уравнений .

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

Вновь рассмотрим операторное уравнение Au = f, где А – линейный непрерывный оператор А: U F, причем U, F – сепарабельные пространства и ker A = {0} .

Будем искать его решение в виде N u N = ak k, (5.55) k =1 причем {k} – заданная линейно независимая система в V. Представление вида (5.55) дает существенную априорную информацию о точном решении и; размерность представления N и коэффициенты этого разложения требуется определить. Отметим, что если при решении корректных задач на основе такого подхода с увеличением числа N базисных функций uN, как правило, приближается к точному решению, то для некорректных задач с неточно заданной правой частью f это не так .

Рассмотрим способы определения uN, если f f. Определим в U систему множеств

–  –  –

и набор сингулярных чисел 1 2... k... N .

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

Теорема. Если базис из собственных элементов упорядочен по невозрастанию сингулярных чисел, то при 0

–  –  –

причем матрица системы есть матрица Грама для системы элементов {k}, которая невырождена, поскольку det[ B ] 0 в силу линейной независимости {k} .

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

Замечание. Из доказательства теоремы следует конструктивный способ выбора N* = N (), при котором выполняются неравенства ( f, m ) 2 q 2 2, ( f, m ) 2 q 2 2 .

k = N () k = N () +1

Задачи:

1. Сформулируйте и докажите теорему о сходимости метода простых итераций для решения некорректных задач .

2. Объясните особенности применения итерационных методов вариационного типа для решения некорректных задач .

3. Поясните суть метода усеченных сингулярных разложений .

4. Поясните суть проекционного метода решения некорректных задач .

Вопросы для самопроверки

1. Для чего применяются методы регуляризации?

2. Что называется квазирешением некорректной задачи .

3. Сформулируйте идею метода квазирешений .

4. В чем суть метода регуляризации Тихонова?

5. Чем метод Тихонова отличается от метода квазирешений?

6. Сформулируйте и докажите теорему о сходимости метода регуляризации Тихонова .

7. Сформулируйте и докажите теорему об априорных оценках приближенного решения .

8. В чем суть выбора параметра регуляризации в методе Тихонова?

9. В чем состоит идея применения итерационных методов для решения некорректных задач?

5.3. НЕКОТОРЫЕ ПРИМЕРЫ РЕШЕНИЯ ОБРАТНЫХ ЗАДАЧ

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

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

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

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

–  –  –

Сформулируем обратную ретроспективную задачу: найти начальное распределение температуры ( x) по некоторой информации о поле температур.

При этом различают две основных постановки обратной задачи:

1. Известно поле температур в некоторый момент времени Т:

–  –  –

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

–  –  –

Отметим, что ядро К1 (х, ) является симметричным, а оператор К1 – самосопряженным в L2[0, l]. Собственные функции ядра –1, cos nx – представляют собой полную ортогональную систему в L2[0, l], а сингулярными числами оператора К1 являются

–  –  –

Анализируя представление (5.62), (5.63), заметим, что для сходимости ряда в представлении (х), в силу наличия у элементов быстро растущих множителей exp( 2 aT ), необходимо требовать довольно жестких n

–  –  –

а число N (в данном случае параметр регуляризации) выбирается согласовано с погрешностью измерения функции f ( f f ). Найдем, из каких условий оно определяется. Оценим погрешность в норме L2[0, l]:

–  –  –

представляет собой норму остатка сходящегося ряда и, следовательно, стремится к 0 при N .

Таким образом, погрешность регуляризованного решения оценивается через величину m( N,) = M ( N ) + C ( N ), и по заданному всегда найдется N* (), доставляющее min m (N, ).

Если необходимо построить решение обратной задачи с погрешностью, не превышающей, при погрешности задания входных данных, то порядок величины N таков:

–  –  –

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

Замечание. В принципе для решения ретроспективной задачи достаточно сделать замену переменных = Т – t, где Т – известное время наблюдения объекта исследования, и решить задачу с обобщенным временем. Соответствующая краевая задача относительно функции v() = (T ) имеет вид

–  –  –

Для определения v (x, Т) = (x) (температуры стержня в момент Т) необходимо решить краевую задачу (5.66) .

Таким образом, искомая функция находится либо из уравнения Фредгольма I рода с гладким ядром (5.60), либо при непосредственном решении (5.66) методом разделения переменных в виде ряда. Члены этого ряда, как указывалось выше, содержат экспоненциально растущие сомножители, и для его сходимости нужны очень жесткие ограничения на функцию f (х) в виде очень быстрого убывания коэффициентов ее ряда Фурье .

Далее покажем, как такому решению придать смысл .

–  –  –

Для любого 0 оператор R определен на всем пространстве L2[0, l] и непрерывен, если его рассматривать действующим из L2[0, l] в L2[0, l] в силу сходимости ряда справа при 0. Отметим, что при таком подходе функцию (x) = Rf можно рассматривать в качестве приближенного решения уравнения (5.69) .

Этот подход весьма эффективен при решении многих других некорректных задач .

Вторая постановка обратной ретроспективной задачи. Рассмотрим краевую задачу для уравнения теплопроводности:

–  –  –

где х0 – некоторая точка внутри отрезка [0, l], что соответствует установке датчика температуры внутри тела .

Воспользуемся решением прямой задачи в виде ряда (5.59), (5.60), удовлетворяя условию (5.72), получим для определения функции () стандартную некорректную задачу – интегральное уравнение Фредгольма

I рода с гладким ядром:

K 2 = K 2 (t,)()d = g ( x ), t [t1, t2 ], (5.74)

–  –  –

В представлении ядра (5.75) х0 – точка наблюдения (измерения температуры), и ее выбор имеет большое значение при исследовании вопроса о единственности решения обратной задачи. Отметим, что (5.74) уравне

–  –  –

Итак, ядро интегрального оператора в (5.74) – непрерывная функция .

Таким образом, рассматриваемая обратная задача сведена к проблеме решения уравнения Фредгольма I рода с непрерывным ядром, которая рассмотрена выше, и требует регуляризации при решении. Отметим, что оператор K2, в отличие от первой постановки обратной задачи, является несамосопряженным. Кроме того, уравнение (5.74) в рассматриваемом случае может иметь единственное решение, а может и не единственное. Оказывается, это в значительной степени зависит от выбора точки х0 [0, l], где происходит измерение температуры.

При этом возможны следующие ситуации:

1. Точка х0 совпадает с одним из концов отрезка [0, l] (например, х0 = 0);

в этом случае решение единственно. Действительно, пусть t [t1, t2]. Введем функцию комплексного переменного l ( z ) = K ( z,)()d,

–  –  –

Возьмем параметр (0, t1), тогда в полуплоскости Re z (z) – аналитическая функция, так как она представима равномерно сходящимся рядом из аналитических функций. В области Re z [ t1, t2] (z) = 0, тогда по теореме о единственности продолжения аналитической функции (z) = 0 всюду в полуплоскости Re z. Следовательно, ()cos( )d = 0, n = 0,1, 2.. .

n Так как система функций {1, cosn} полна на [0, l], то () 0 на этом отрезке и имеет место единственность решения обратной задачи .

2. Пусть точка наблюдения отлична от конца стержня, например l x0 =. Выясним, будет ли единственным образом восстановлено начальное распределение температуры в этом случае. Пусть

–  –  –

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

Регуляризованные разностные схемы. Важный класс методов решения некорректных задач связан с возмущением исходного уравнения .

В методе квазиобращения такое возмущение проводится для оператора исходной дифференциальной задачи. Более естественный подход основан на возмущении сеточной задачи. По этому принципу строятся регуляризованные разностные схемы [2] .

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

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

При этом безусловно устойчивые разностные схемы реализуются следующим образом [2]:

1. Для исходной задачи строится простейшая разностная схема (производящая разностная схема), не обладающая необходимыми свойствами, т.е .

схема является условно устойчивой либо даже абсолютно неустойчивой .

2. Разностная схема записывается в канонической форме, для которой условия устойчивости известны .

3. Устойчивость разностной схемы улучшается за счет возмущения оператора разностной схемы .

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

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

Рассмотрим двумерную задачу

–  –  –

Будем считать, что коэффициент k достаточно гладкий и k ( x), 0, x .

Для дифференциальной задачи (5.77)–(5.79) запишем соответствующую дифференциально-разностную задачу, полученную после дискретизации по пространству на равномерной сетке с шагами h, = 1, 2. Пусть – множество внутренних узлов .

На множестве сеточных функций y (x), таких, что y (x) = 0, x, определим сеточный оператор:

y = (a y x ) x =1

–  –  –

Построим безусловно устойчивые двухслойные разностные схемы для (5.81), (5.82) на основе принципа регуляризации .

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

–  –  –

В соответствии с (5.86) повышение устойчивости может достигаться либо за счет увеличения энергии (By, y) оператора B, либо за счет уменьшения энергии (Аy, y) оператора А .

Рассмотрим вначале аддитивную регуляризацию (регуляризация за счет добавления операторных слагаемых в операторы А и В). Начнем с аддитивного возмущения оператора В, т.е. с перехода B B + R, где R – регуляризирующий оператор, а – параметр регуляризации .

Поскольку, в данном случае В = Е, положим

–  –  –

В случае R = A условие устойчивости имеет вид /8. Другой пример такой регуляризации соответствует попеременно-треугольному методу, когда A = R* + R и /2 .

Аналогично проводится мультипликативная регуляризация за счет возмущения оператора А. С учетом неравенства (5.86) можно осуществить преобразование A A (E + R)–1 или A (E + R)–1 A. Для простейших двухслойных схем такая регуляризация может рассматриваться как новая редакция регуляризации оператора В. Для того чтобы остаться в классе схем с самосопряженными операторами, достаточно выбрать R = (А) .

Большие возможности предоставляет регуляризация А (E + R*)–1 А (E + R)–1 .

В этом случае регуляризирующий оператор R может напрямую не связываться с оператором А .

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

–  –  –

которое отличается от (5.77) только знаком при производных по пространству (соответствует замене t на –t). Граничные и начальные условия остаются прежними:

u (x, t ) = 0, x, 0 t T,

–  –  –

Определение. Разностная схема называется -устойчивой (равномерно устойчивой) по начальным данным в Н, если существуют постоянная 0 и постоянная m1, не зависящие от, n, такие, что при любых n и при всех ynH для решения yn+1 разностного уравнения справедлива оценка

–  –  –

= 1 + М. (5.95) Доказательство. Этот результат следует из общих условий -устойчивости двухслойных операторно-разностных схем. Разностная схема (5.84), (5.85) с самосопряженными и постоянными операторами B 0, A, будет -устойчива в НВ при

–  –  –

и с учетом (5.94) имеет место при выборе согласно (5.95) .

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

Исходя из явной схемы (5.84), (5.92) для задачи (5.82), (5.91) запишем регуляризованную схему в каноническом виде (5.85), где B = E + R, A =. (5.97) Теорема. Регуляризованная схема (5.85), (5.97) -устойчива в НВ с

–  –  –

Это неравенство будет выполнено при, удовлетворяющем (5.99) .

Аналогично строятся и другие регуляризованные разностные схемы .

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

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

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

Коэффициентная обратная задача для эллиптического уравнения. Постановка прямой задачи имеет следующий вид: определить функцию и (х), х = (х1, х2), удовлетворяющей уравнению

–  –  –

Младший коэффициент с зависит только от переменной х2 .

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

Решение задачи ищется в прямоугольной области

–  –  –

Покажем, что равенства (5.107)–(5.109) верны только при v (x) = 0, (x2) = 0, x .

Задача (5.107)–(5.109) рассматривается при заданных с1 (х2) и u2 (x) .

Это обратная (линейная) задача по определению пары v (x), (x2), x – задача идентификации правой части. Сформулируем некоторые достаточные условия, при которых можно гарантировать v (x) = 0, (x2) = 0, x .

Будем считать решение и коэффициент достаточно гладкими, кроме того, в (5.101)–(5.103) решение u (x) является знакопостоянной функцией

–  –  –

(5.113)–(5.115) стандартная краевая задача с единственным решением .

Осталось только указать вид функции (х). По найденной из (5.113)–(5.115) функции w (x) для этого решается краевая задача

–  –  –

с граничными условиями (5.111) .

Сеточная обратная задача. Запишем разностный аналог коэффициентной обратной задачи (5.101)–(5.103). В области введем равномерную по каждому направлению сетку; определим множество внутренних узлов сетки:

–  –  –

В обратной задаче сеточная функция с(x2), x22 неизвестна, но заданы дополнительные условия, которые соответствуют (5.103). При переходе к дискретной задаче дополнительные условия можно сопоставить с заданием сеточного решения в узлах *. Это имеет место при использовании простейшей двухточечной аппроксимации направленными разностями с погрешностью аппроксимации первого порядка. Тогда дополнительные условия можно привести к виду y (x) = (x), x *. (5.118) При решении прикладных задач входные данные задаются с погрешностью. Будем считать, что по сравнению с ними погрешностями аппроксимации можно пренебречь (используя достаточно подробные сетки) .

Здесь будем считать, что основные погрешности вносятся измерениями нормальной производной на границе. Поэтому при приближенном решении обратной задачи с погрешностями во входных данных граничное условие (5.117) оставим без изменения, а вместо (5.118) положим y ( x) ( x ), x *, (5.119) где параметр задает уровень погрешностей .

Итерационное решение обратной задачи. Для приближенного решения обратной (5.116), (5.117), (5.119) задачи будем использовать итерационные методы. В этом случае искомая сеточная функция y ( x ) ( x ), x * уточняется на каждом итерационном шаге исходя из критерия минимизации функционала невязки (точности выполнения условия (5.119)) .

Определим невязку в виде [2]

–  –  –

При получении выражения для градиента функционала J (c) примем, что приращению с соответствуют приращения J и у для функционала (5.120) и решения задачи (5.116), (5.117) .

С точностью до членов второго порядка малости из (5.116), (5.117) получим

–  –  –

Для того чтобы преобразовать правую часть (5.123), домножим уравнение (5.121) на некоторую сеточную функцию (х) h1h2, х и просуммируем его по всем узлам :

–  –  –

Для получения представления для J’(c) первый член (5.126) связывается с правой частью равенства (5.123) .

Пусть функция (х) определяется как решение уравнения

–  –  –

Его вычисление связано с решением краевой задачи (5.116), (5.117) для основного состояния (сеточная функция y (x)) и краевой задачи (5.125), (5.127) для сопряженного состояния (сеточная функция (x)) .

При использовании двухслойного градиентного итерационного метода уточнения ck (x2), где k – номер итерации, проводится по схеме

–  –  –

Необходимо только учитывать, что решаемое уравнение J’(c) = 0 является нелинейным. Для выхода из итерационного процесса (5.129) можно использовать критерий невязки .

Задачи:

1. Запишите регуляризованную схему для простейшей явной разностной схемы. Сформулируйте условия ее устойчивости .

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

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

Вопросы для самопроверки

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

2. Сформулируйте две постановки обратной ретроспективной задачи .

В чем различия с физической точки зрения?

3. Получите решение обратной задачи теплопроводности (первая постановка) в виде интегрального уравнения Фредгольма I рода. Объясните, в чем проявляется некорректность полученного решения .

4. Объясните идею применения метода усеченных сингулярных разложений для регуляризации решения задачи с обратным временем (в первой постановке). Что является параметром регуляризации и как он выбирается?

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

6. Получите решение обратной задачи теплопроводности (вторая постановка) в виде интегрального уравнения Фредгольма I рода. Объясните, в чем проявляется некорректность полученного решения .

7. Как влияет выбор точки наблюдения температуры (дополнительное условие для решения обратной задачи) на единственность решения обратной задачи?

8. В чем суть принципа регуляризации разностных схем?

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

СПИСОК ЛИТЕРАТУРЫ К ГЛАВЕ 5

1. Кабанихин С.И. Обратные и некорректные задачи. – Новосибирск:

Сибирское научное изд-во, 2009. – 457 с .

2. Самарский А.А., Вабищевич П.Н. Численные методы решения обратных задач математической физики: учеб. пособие. – изд. 3-е. – М.:

ЛКИ, 2009. – 480 с .

3. Ватульян А.О. Обратные задачи в МДТТ. – М.: ФИЗМАТЛИТ, 2007. – 224 с .

4. Самарский А.А., Гулин А.В. Численные методы: учеб. пособие для вузов. – М.: Наука, 1989. – 432 с .

5. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. – изд. 2-е. – М.: Наука, 1979. – 285 с .

6. Латтес Р., Лионс Ж.-Л. Метод квазиобращения и его приложения. – М.: Мир, 1970 .

7. Бакушинский А.Б., Гончарский А.В. Некорректные задачи. Численные методы и приложения. – М.: Изд-во Моск. ун-та, 1989. – 199 с .

–  –  –

__________________________________________________________

Подписано в печать 5.12.2012. Формат 6090/16 .

Усл. печ. л. 29,5. Тираж 10 экз. Заказ № 261/2012.

Pages:     | 1 ||



Похожие работы:

«Заседание Координационного совета по промышленности и научнотехнической политике Межрегиональной ассоциации "Сибирское соглашение" 9 апреля 2013 г. г. Новосибирск СИБИРСКОЕ ОТДЕЛЕНИЕ РАН: РЕЗУЛЬТАТЫ ФУНДАМЕНТАЛЬНЫХ ИССЛЕДОВАНИЙ...»

«РУКОВОДСТВО ПОЛЬЗОВАТЕЛЯ IP ВИДЕОКАМЕРЫ серий PN2X-XX-XXXX PD2X-XX-XXXX PS2X-XX-XXXX PQ2X-XX-XXXX PN7X-XX-XXXX PD7X-XX-XXXX PS7X-XX-XXXX PQ7X-XX-XXXX LC-XXXX             Спасибо, что выбрали товар нашей торговой марки. Внимание! Диза...»

«Столбов М.И. Финансовый рынок как передаточный механизм нестабильности в открытой экономике / М.И. Столбов // Экономические науки. 2008. № 8. С. 337-340. М.И. Столбов, к.э.н. (специальность 08.00.01-Экономическая теория), п...»

«МИНИСТЕРСТВО СТРОИТЕЛЬСТВА И ЖИЛИЩНО-КОММУНАЛЬНОГО ХОЗЯЙСТВА РОССИЙСКОЙ ФЕДЕРАЦИИ СВОД ПРАВИЛ СП 360.1325800.2017 КОНСТРУКЦИИ СТАЛЕФИБРОБЕТОННЫЕ Правила проектирования Издание официальное Москва Стандартинформ сертификат клапан СП 360.1325800.2017 Предисловие...»

«RG-P58D Руководство пользователя POS-принтер печати чеков RG-P58D Руководство пользователя RG-P58D Руководство пользователя ОГЛАВЛЕНИЕ Меры безопасности Глава I. Внешний вид и модели 1.1 Внешний вид 1.2 Маркировка модификаций моделей Глава II. Особенности 2.1 Технические харак...»

«Организация Объединенных Наций A/HRC/WG.6/6/KHM/1 Генеральная Ассамблея Distr.: General 16 September 2009 Russian Original: English Совет по правам человека Рабочая группа по универсальному п...»

«Институт Строительства и Архитектуры Е. Ю. БАГИНА ЛАНДШАФТ: КОМПОЗИЦИОННЫЕ АСПЕКТЫ Учебное пособие МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ УРАЛЬСКИЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ ИМЕНИ ПЕРВОГО ПРЕЗИДЕНТА РОССИИ Б. Н. ЕЛЬЦИНА Е. Ю. Багина ЛАНДШАФТ: КОМПОЗИЦИОННЫЕ АСПЕКТЫ Учебное п...»

«J Версия 1.3 RU БЛАГОДАРИМ ЗА ВЫБОР ИЗДЕЛИЯ COWON. Желаем вам приятного использования концептуального устройства серии “Digital Pride”. Данное руководство поможет вам ознакомиться с функци...»

«СОДРУЖЕСТВО НЕЗАВИСИМЫХ ГОСУДАРСТВ МЕЖГОСУДАРСТВЕННЫЙ СОВЕТ ПО СТАНДАРТИЗАЦИИ, МЕТРОЛОГИИ И СЕРТИФИКАЦИИ ПРОТОКОЛ № 53-2018 пятьдесят третьего заседания Межгосударственного совета по стандартизации, метрологии и сертификации г. Ташкент 26 28 июня 2018 г. СОСТАВ МЕЖГОСУДАРСТВЕННОГО СОВЕТА ПО СТАНДАРТИЗАЦИИ, МЕТРОЛОГИИ И СЕРТИФИКАЦИИ Полномочным...»

«Министерство образования и науки Российской Федерации федеральное государственное автономное образовательное учреждение высшего образования "НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ТОМСКИЙ ПОЛИТЕХНИЧ...»

«МЕТОДЫ НЕЙРОСЕТЕВОГО И МУЛЬТИАГЕНТНОГО УПРАВЛЕНИЯ В РОБОТОТЕХНИКЕ И МЕХАТРОНИКЕ А.В.Тимофеев Санкт-Петербургский государственный университет аэрокосмического приборостроения Санкт-Петербургский институт информатики и автоматизации РАН, ВВЕДЕНИЕ Развитие машиностроения, авиации, ракетной техники и наукоёмких технологий приве...»

«TKF-12 УКАЗАТЕЛЬ ПРАВИЛЬНОСТИ ЧЕРЕДОВАНИЯ ФАЗ Руководство по эксплуатации Версия 1.11 Серийный номер №_ 1 ВВЕДЕНИЕ 2 ТЕХНИЧЕСКИЕ ДАННЫЕ 3 КОМПЛЕКТАЦИЯ 4 ОПРЕДЕЛЕНИЕ ЧЕРЕДОВАНИЯ ФАЗ 5 ОБСЛУЖИВАНИЕ УКАЗАТЕЛЯ 6 ГАРАНТИИ 6.1 Общие положения гарантийного...»

«Министерство образования и науки Российской Федерации федеральное государственное автономное образовательное учреждение высшего образования "НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ТОМСКИЙ ПОЛИТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ" Инженерная школа новых производственных т...»

«RX 000 рублей AWD от 3 677 ЦВЕТА КУЗОВА1 Белый металлик (083) Искрящийся белый металлик (085) Серый металлик (1H9) Серебристый металлик (1J4) Светло-серый металлик (1J7) Черный неметаллик (212) Черный металлик (223) Красный металлик (3R1) Бежевый металлик (4U7) Темно-синий металлик (8X5) Коричневый металлик (4X2) Лак...»

«Министерство образования и науки Российской Федерации ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ОБРАЗОВАНИЯ "САРАТОВСКИЙ НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ИМЕНИ Н.Г. ЧЕРНЫШЕВСКОГО" Кафедра математического...»

«Секция 2: Инновационные технологии получения и обработки материалов в машиностроении Условия протекания реакции между указанными компонентами представлены минимальной температурой (2500 К) и наименьшей внутренней энергией (69693 кДж/кг). Список литературы 1. Тен Э.Б. Осно...»

«0512938 ТЕХНИЧЕСКИ ЦЕНТР БЕЗОПАСНОСТИ 'S\ КАТАЛОГ ОБОРУДОВАНИЯ Технический Центр Пожарной Безопасности основан 19 октября 1995 года, Продукция Техцентра патентованные порошковые зарядные станции СЗП-02ГУ и огнезащитный порошок ПО-ТМ15, огнезащитный пропиточный состав ОГЗС-1, углекислотные зарядные ста...»

«Международный инновационный центр PERSPEKTIVA PLUS Культура, просвещение, литература MSBooks Publishing, г. Виннипег, Канада ISBN: 978-0-9877600-9-8 Составители: Козлова Л.М., член Союза писателей СССР, г. Бийск, Россия Тимохин Н. Н., член Союза писателей России, г. Семей, Казахстан Реда...»

«ISSN 1811-1807 РЫЛЫМИ ЖУРНАЛ С ТОРАЙГЫРОВ АТЫНДАРЫ s.:X s m aekettik em JfiP L Павлодар УНИВЕРСИТЕТ! ФКАААКР ИАТ ТЛСИ З Е ИЫИ НМ К Е М ПМУ ХАБАРШЫСЫ ВЕСТНИК ПГУ УДК 539.3:534.2 Н. А. Испулов, А....»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ ФГБОУ ВПО "АЛТАЙСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ" ТЕОРИЯ И ПРАКТИКА АРХЕОЛОГИЧЕСКИХ ИССЛЕДОВАНИЙ №1 (7) ISSN 2307-2539 Журнал основан в 2005 г. Выходит 2 раза в год Главный редактор: А.А. Тишкин Редакционная коллегия:...»

«Ru Специальное дополнение к руководству по монтажу и эксплуатации на ротаметры Н250 /Н54 Электронный конвертор M10 взрывозащищенного исполнения EEx-d PTB 01 ATEX 1154 Ротаметры Вихревые расходомеры Контроллеры расхода Электромагнитные расхо...»

«Ермашов Алексей Олегович ГЕОМЕХАНИЧЕСКОЕ ОБОСНОВАНИЕ РАСЧЕТОВ ОСЕДАНИЙ ЗЕМНОЙ ПОВЕРХНОСТИ ПРИ ДОБЫЧЕ КАЛИЙНО-МАГНИЕВЫХ РУД (НА ПРИМЕРЕ ВЕРХНЕКАМСКОГО МЕСТОРОЖДЕНИЯ КАЛИЙНО-МАГНИЕВЫХ СОЛЕЙ)....»

«Овтина Евгения Алексеевна РЕДУПЛИКАЦИЯ НА УРОКАХ РУССКОГО ЯЗЫКА КАК ИНОСТРАННОГО Статья посвящена представлению редупликации на уроках русского языка как иностранного. Новизна исследования заключается в том, что редупликация рассматривается как активный механиз...»







 
2019 www.librus.dobrota.biz - «Бесплатная электронная библиотека - собрание публикаций»

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