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

Pages:   || 2 |

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

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

Министерство образования и науки Российской Федерации

Федеральное государственное бюджетное

образовательное учреждение

высшего профессионального образования

«Пермский национальный исследовательский

политехнический университет»

И.Ю. Зубко, Н.Д. Няшина

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ:

ДИСКРЕТНЫЕ ПОДХОДЫ

И ЧИСЛЕННЫЕ МЕТОДЫ

Утверждено

Редакционно-издательским советом университета в качестве учебного пособия Издательство Пермского национального исследовательского политехнического университета УДК 539.3 З-91

Рецензенты:

д-р физ.-мат. наук, профессор А.А. Роговой (Институт механики сплошных сред УрО РАН, г. Пермь);

д-р техн. наук, профессор В.Н. Аптуков (Пермский государственный национальный исследовательский университет) Зубко, И.Ю .

З-91 Математическое моделирование: дискретные подходы и численные методы : учеб. пособие / И.Ю. Зубко, Н.Д. Няшина. – Пермь :

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

ISBN 978-5-398-00947-7 В учебном пособии на примерах исследования физико-механических свойств кристаллических материалов (образцы с ОЦК-, ГЦК-, ГПУ-решетками, графит, графен) продемонстрированы возможности и обоснованы границы применимости как физического дискретного подхода, так и математического дискретного подхода. Достоинства и недостатки первого иллюстрируются на примере атомистических методов исследования термомеханических свойств кристаллических твердых тел и на примере метода клеточных автоматов. Особенности второго подхода демонстрируются на примерах численных методов решения уравнений механики сплошных сред – методе конечных элементов, методах решения обратных и некорректных задач .

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

Изучение материала, включенного в учебное пособие, предусмотрено в 10–11 семестрах учебного плана профиля магистратуры «Математическое моделирование физико-механических процессов» по направлению 010400.68 Прикладная математика и информатика .

УДК 539.3 © ПНИПУ, 2012 ISBN 978-5-398-00947-7 ОГЛАВЛЕНИЕ Введение

Глава 1. ОСНОВНЫЕ СВЕДЕНИЯ ИЗ ТЕНЗОРНОГО АНАЛИЗА И МЕХАНИКИ СПЛОШНОЙ СРЕДЫ

1.1. Основы векторной и тензорной алгебры

1.2. Основные сведения из механики

1.3. Балансовые уравнения и законы механики сплошных сред

1.3.1. Основные соотношения кинематики сплошной среды

1.3.2. Основные соотношения динамики сплошной среды

1.4. Практическое введение в пакет «Mathematica»

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

Глава 2. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПОВЕДЕНИЯ

ДЕФОРМИРУЕМЫХ ТВЕРДЫХ ТЕЛ НА АТОМАРНОМ УРОВНЕ .

ДИСКРЕТНЫЙ ПОДХОД

2.1. История дискретного подхода. Метод молекулярной динамики .

Сравнение различных потенциалов. Численные методы

2.2. Метод атомарной статики. Определение упругих модулей материалов с различной кристаллической решеткой

2.3. Учет температуры. Зависимость механических свойств конденсированных сред от температуры

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

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

3.1. Моделирование с использованием имитационного подхода

3.2. Метод клеточных автоматов. Основные определения





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

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

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

4.1. Проекционный и вариационный подходы к построению разрешающих соотношений МКЭ

4.1.1. Проекционные методы. Основные понятия.

4.1.2. Вариационные принципы

4.2. Виды конечных элементов и их классификация

4.2.1. Одномерные пробные функции

4.2.2. Двумерные базисные функции высших степеней

4.2.3. Трехмерные базисные функции

4.2.4. Разбиение области на элементы. Требования к конечным элементам

4.3. Отображение и численное интегрирование

4.3.1. Параметрическое отображение

4.3.2. Субпараметрические, изопараметрические и суперпараметрические элементы

4.3.3. Численное интегрирование

4.4. Разрешающие соотношения МКЭ для квазистатических задач теории упругости

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

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

4.5.2. Разрешающие соотношения МКЭ для пространственной и временной конечно-элементной аппроксимации

4.6. Физически нелинейные задачи. Пластичность, ползучесть

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

4.6.2. Методы начальных (дополнительных) напряжений

4.6.3. Методы начальных (дополнительных) деформаций

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

4.7. Геометрически нелинейные задачи

4.7.1. Общие положения

4.7.2. Общий случай больших деформаций и напряжений

Выводы по главе 4

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

Глава 5. ОБРАТНЫЕ И НЕКОРРЕКТНЫЕ ЗАДАЧИ

5.1. Основные понятия и примеры

5.1.1. Обратные и некорректные задачи

5.1.2. Понятие условно корректной задачи

5.2. Способы преодоления некорректности. Методы регуляризации

5.2.1. Метод квазирешений

5.2.2. Метод регуляризации Тихонова

5.2.3. Метод регуляризации на компактных множествах

5.2.4. Итерационные методы решения некорректных задач

5.2.5. Метод усеченных сингулярных разложений

5.2.6. Проекционный метод

5.3. Некоторые примеры решения обратных задач

5.3.1. Обратные ретроспективные задачи

5.3.2. Коэффициентные обратные задачи в механике деформируемого твердого тела

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

ВВЕДЕНИЕ

В математическом моделировании физико-механических систем, начиная с первых исследований простейших упругих реакций деформируемых тел, сосуществовали два подхода к теоретическому исследованию природы – дискретный (молекулярная теория Навье, Коши, Пуассона, 1827 год) и непрерывный (теория потенциала Грина, 1838 год) .

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

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

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

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

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

Цель предлагаемого пособия – продемонстрировать ряд возможностей и указать границы применимости как физического дискретного подхода, так и математического. Достоинства и недостатки первого подхода обсуждаются на примере атомистического подхода к исследованию Рippan R., Gumbsch P. Multiscale modelling of plasticity and fracture by means of dislocation mechanics // CISM. Courses and lectures. Springer Wien New York, 2010. No. 522. 394 p .

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

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

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

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

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

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

М.: Логос, 2005. 440 с .

ГЛАВА 1. ОСНОВНЫЕ СВЕДЕНИЯ ИЗ ТЕНЗОРНОГО

АНАЛИЗА И МЕХАНИКИ СПЛОШНОЙ СРЕДЫ

Для достижения замкнутости изложения материала в п. 1.1. кратко обсудим основы математического аппарата (тензорной алгебры), необходимого для изучения последующих разделов. При этом будем выделять векторы и тензоры жирным шрифтом. Векторы будем обозначать малыми латинскими буквами a, b, x, y и так далее, а тензоры – как правило, заглавными, например A, B, T, R, скалярные величины будут обозначаться обычным шрифтом, например,,, i, k. В п. 1.2 и 1.3 приводятся основные сведения из механики. Основной задачей изложения материала п. 1.3 является представление уравнений механики сплошных сред (МСС), термодинамики и получение уравнений диффузии в растворе двух и более компонент. В последующих главах при моделировании механического и теплового поведения дискретных систем осредненные результаты расчетов будут сравниваться с решениями соответствующих уравнений МСС для непрерывных систем, приведенных в п. 1.3, в гл. 4 и 5 будут рассмотрены специальные разделы теории численных методов для решения уравнений из п. 1.3 при различных граничных условиях. П. 1.4 посвящен практическому введению в вычислительную систему «Mathematica» разработчика Wolfram Research, которую предлагается использовать в качестве инструмента при освоении материала последующих глав .

1.1. ОСНОВЫ ВЕКТОРНОЙ И ТЕНЗОРНОЙ АЛГЕБРЫ

–  –  –

то есть площади параллелограмма, построенного на векторах x и y, как на сторонах, а направление z определяется по «правилу буравчика» перпендикулярно обоим векторам (рис. 1.2) – поворот от вектора x к вектору y с конца полученного вектора z видится как происходящий против часовой стрелки. При перестановке сомножителей в векторном произведении направление вектора z меняется на противоположное, и поворот от вектора x к вектору y будет представляться происходящим по часовой стрелке. В координатной форме результат векторного умножения имеет более сложный вид:

–  –  –

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

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

Так и векторы, родившись в механике, с самого начала стали наделяться вполне ясным физическим смыслом: вектор скорости (точки тела) v, вектор силы f, вектор момента M, вектор нормали к поверхности n и многие другие. Некоторые примеры оказались довольно необычными .

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

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

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

Причем, не вдаваясь в детали, будем считать, что именно единичный вектор n и задает площадку, то есть любому вектору n, проведенному из данной точки, будет соответствовать единственная элементарная площадка с заданной (очень малой) площадью S. Поэтому естественным оказывается появление векторной площадки S n S .

Другим важным объектом векторной природы является так называемый вектор напряжений или усилий (не путать с силой). Для пояснения его физического смысла можно привести пример опытов Галилео Галилея. Исследуя прочность канатов различной толщины на разрыв, он обнаружил, что величина приложенной силы f, при которой рвется канат, является недостаточно удобной характеристикой прочности. Это вполне понятно, поскольку в ней никак не учитывалась толщина каната. Отношение f к площади поперечного сечения каната S имеет размерность Па (как и давление) и называется вектором напряжений (называемый в некоторых работах также вектором усилий) p = f / S. Это понятие впервые было введено Галилеем и оказалось впоследствии очень плодотворным .

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

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

Во-первых, отметим, что все тела, окружающие нас, состоят из огромного (дискретного) набора мельчайших взаимодействующих частиц. Записывая для каждой из них законы Ньютона с учетом законов взаимодействия этих частиц (основа дискретного подхода), мы, видимо, сможем описать поведение тела, но при этом получим систему огромного количества уравнений, решить которую для макроскопических объектов (строительных конструкций, взаимодействующих деталей двигателей и других механизмов) практически невозможно. Другой подход предоставляет статистическая физика, изучающая статистические законы распределения параметров, описывающих, например, движение и взаимодействие частиц тела. Третий подход основан на использовании модели сплошной среды (или модели деформируемого континуума). Свойства и характеристики отдельных частиц (масса, температура, скорость, ускорение и все остальные) «размазываются» по объему, который занимает тело. В результате мы имеем дело не с набором частиц, а с непрерывными распределениями (функциями) их характеристик по области пространства, занимаемой телом. Как в эту модель включить взаимодействие частиц, существовавшее до «размазывания»? Ответ на этот вопрос помогает дать принцип разрезания Эйлера–Коши: если мы рассмотрим любую элементарную площадку в теле, то предполагается, что воздействия от всех частиц, оказавшихся по одну сторону, на частицы по другую сторону площадки можно заменить действием поверхностной силы f (а в некоторых теориях – и моментом сил). То есть получается, что каждой площадке в каждой точке тела соответствует своя поверхностная сила. Как мы выяснили, при исследовании деформируемых тел рассматриваются векторные площадки S = n S, то есть каждой нормали n будет соответствовать своя площадка и своя сила f. Чтобы ввести напряжение, можно разделить силу на скалярную величину площади и ввести напряжение в точке, зависящее от ориентации площадки: pn = f / S. Здесь неслучайно используется обозначение pn проекции. Но в этом случае ее результатом является не число, а вектор. Так появляется необходимость в новом математическом понятии, связывающем пару векторов, – тензоре второго ранга, а в рассмотренном случае тензоре напряжений. «Проекция» этого объекта на нормаль n = pn и есть вектор напряжений. Сам тензор напряжений в точке M определяется как совокупность всех возможных упорядоченных пар n и pn, заданных в точке M. И теперь вопрос об учете взаимодействия мельчайших частиц тела сводится к заданию связи тензора напряжений и характеристик движения деформируемого континуума. Для упорядоченных пар n и pn можно было бы ввести обозначение как для компонент вектора (упорядоченной пары чисел, только «компоненты» нового объекта сами являются векторами), например {n, pn }, но оно не вполне удобно. Для обозначения неразделимой упорядоченной пары векторов вводят символ, причем упорядоченность означает, что n pn pn n (как, впрочем, и для векторного умножения). Этот знак называют диадным или тензорным умножением, а упорядоченную пару векторов n pn называют диадой. Заметим, что в современной научной литературе знак для упрощения записи не ставят, предполагая, что по контексту можно понять, о каком умножении идет речь .

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

a b c, будет представлять собой тензор третьего ранга .

Впервые тензоры были введены Леонардом Эйлером в 1758 году (тензор инерции твердого тела, тензор поворота). Хотя термин «тензор»

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

Свойства диадного произведения. Тензоры второго ранга. Тензорное произведение может связывать в новый объект только векторы или тензоры, а если появляется произведение скалярной величины k (числа) на вектор a: k a, то результатом будет вектор ka.

При скалярном умножении диады a b на вектор с слева на первом месте в диаде появляется скаляр – результат скалярного произведения:

c a b = (c a) b = (c a)b = k1b. На втором месте в диаде при этом остается неизменным вектор b.

При скалярном умножении диады a b на вектор с справа на первом месте в диаде остается вектор а, а на втором месте появляется скаляр – результат скалярного произведения:

a b c = a (b c) = a(b c) = ak2 = k2a. Это помогает понять, почему в общем случае n p n p n n .

Все множество произвольных диад a b образует, подобно векторам, так называемое линейное пространство, то есть множество объектов, для которых выполняются правила (аксиомы):

(a b ) = ( a ) b = a ( b ) = a b, (1.4)

–  –  –

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

ab +cd ?

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

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

–  –  –

Как видно, при векторном умножении тензора второго ранга на вектор получается тензор второго ранга .

Часто встречающейся характеристикой тензора второго ранга является его след trA (от английского trace), который является скаляром вида:

–  –  –

где m, n, p – любые (и это очень важно!) единичные попарно ортогональные векторы (в частности, m m = 1, но m n = 0 и m p = 0 ), в качестве которых могут быть выбраны, например, орты i, j, k координатных осей .

Выше было показано, что в результате скалярного умножения тензора второго ранга на вектор получается новый вектор. Поэтому каждый тензор второго ранга можно воспринимать как некоторое действие, переводящее один вектор в другой. В математике такие действия (удовлетворяющие перечисленным для тензоров аксиомам) называют линейными операторами. Любой вектор x можно представить в виде суммы трех взаимно ортогональных векторов, полученных, например, с помощью упомянутых выше векторов m, n, p – x = x1 + x 2 + x3, где x1 = (x m)m, x 2 = (x n)n, x3 = (x p)p. Такие преобразования, как поворот, растяжение или их комбинация, переводящие вектор p = p1 + p 2 + p 3 ( pi попарно перпендикулярны) в вектор d = d1 + d 2 + d 3 (di попарно перпендикулярны), представляются как тензор второго ранга Q вида Q= d1 p 1 + d2 p2 + d3 p3. (1.15) p 1 p1 p2 p2 p3 p3 Действительно, Q p = d. Преобразование растяжения в k раз есть тензор kE. Также можно построить тензор, задающий поворот на угол ( ) относительно оси, направление которой определяется единичным вектором m. Структура этого ортогонального тензора задается теоремой Эйлера [7], которая говорит, что тензор поворота Q произвольного вектора x относительно оси m ( Q x ) имеет вид

–  –  –

Рис. 1.3. Поворот вектора вокруг оси m: а – поворот в пространстве;

б – на плоскости, ортогональной вектору m (вид из точки O на конце оси m) Используя теорему Эйлера, получим правило наложения поворотов .

Идея очень проста: повернем вектор x вокруг оси m на угол 1 и получим вектор Q m1 x, а затем повернем полученный вектор еще и на угол 2 относительно той же оси m, тогда

–  –  –

В преобразованиях использовались соотношения, в которых принималось, что для попарно ортогональных единичных векторов m, n, p выполняется ( m n = p, m p = n, n p = m ):

<

–  –  –

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

(Q a) (Q b) = Q (a b) .

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

–  –  –

Тензор второго ранга A называется проектором, если выполняются условия AT = A и A A = A. (1.20) Приведем примеры тензоров-проекторов: m m, m m + n n, E = m m + n n + p p (m, n, p – взаимно ортогональные единичные векторы). Название «проектор» можно «почувствовать» на этих трех тензорах .

Оно обусловлено тем, что при скалярном умножении одного из них на вектор x получается либо вектор-проекция x на вектор m, то есть на прямую, вдоль которой направлен m (в случае тензора m m ): x m m = (x m ) m = = (x m)m = xmm (вспомните определение проекции через скалярное умножение векторов), либо вектор-проекция x на плоскость векторов m и n (для тензора m m + n n ): x m m + x n n = (x m) m + (x n) n = = xmm + xnn, либо сам вектор x, как его вектор-проекция в трехмерное пространство (в случае тензора E ) .

Действие над вектором, приводящее к отражению вектора x = {x1; x2 ; x3} = x1e1 + x2e 2 + x3e3 x1 + x 2 + x3 (заметим, что все векторы x1, x2, x3 попарно ортогональны) относительно плоскости x1O x3 (перпендикулярной вектору e2, рис. 1.4), называется преобразованием отражения, и в компонентной записи оно приводит к результату {x1; x2 ; x3} = x1e1 x2e 2 + x3e3.

Такое отражение задается тензоромпроектором B более сложной структуры, чем рассмотренные:

B = E 2e 2 e 2 .

–  –  –

Проектор вида e2 (E e2 e2 ) сначала укладывает любой вектор x на плоскость X 1O X 3 (за счет действия выражения в скобках), а затем с помощью векторного умножения на e2 производит вектор, перпендикулярный e2, то есть принадлежащий плоскости X 1O X 3, и одновременно перпендикулярный лежащему в этой плоскости вектору (E e2 e2 ) x. Так реализуется поворот на 90° против часовой стрелки векторной проекции x в плоскости X 1O X 3. Раскрыв скобки, получим e 2 (E e 2 e 2 ) = e 2 E e 2 e 2 e 2 = = e2 (e1 e1 + e2 e2 + e3 e3 ) 0 = = (e 2 e1 ) e1 + (e 2 e 2 ) e 2 + (e 2 e3 ) e3 = = e3 e1 + e1 e3 .

Как видно, первый вариант представления тензора поворота (1.16) содержит три проектора m m, E m m и m ( E m m ), принцип действия которых мы рассмотрели. С их помощью структура тензора поворота становится более понятной .

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

1.2. ОСНОВНЫЕ СВЕДЕНИЯ ИЗ МЕХАНИКИ

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

То есть система отсчета является моделью, позволяющей фиксировать события, происходящие в окружающем нас мире. Далее тело отсчета мы не будем упоминать, подразумевая его существование и связь с системой координат. В механике Ньютона обычно принимается, что существует и может быть определена некая неподвижная система координат. Часто в качестве нее берется система координат с началом в центре Солнца и осями, направленными на «бесконечно удаленные» звезды. Но считать ее неподвижной, строго говоря, у нас нет никаких оснований. Даже само понятие неподвижности требует ответа на вопрос, по отношению к какой системе отсчета определяется неподвижность? Здесь мы будем исходить из несколько иных посылок (следуя работам таких ученых, как Л. Эйлер, Л.И. Седов, К. Трусделл [11–13]). Для каждого наблюдателя (в качестве которого может выступать любой из читателей или его коллега) в пространстве событий имеется своя собственная система отсчета и собственные «часы». Будем считать для простоты, что системы координат систем отсчета не деформируются (хотя достаточно, чтобы они деформировались единым образом) и «одинаковы». То есть они могут быть совмещены трансляцией (переносом в пространстве, при котором любая из координатных линий остается параллельной себе) и поворотом. Течение времени будем считать одинаковым, а масштабы (шкалы) «часов»

совпадающими, хотя абсолютное время в разных системах отсчета может и различаться (t = t + a, a = const ). При этом ни один из наблюдателей не может утверждать (абсолютную) неподвижность своей системы отсчета или какой-либо другой – речь может идти только об относительном движении (одной системы отсчета относительно другой, выбранной) .

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

Например, радиусы-векторы точки M (векторы, соединяющие начала соответствующих систем координат O и O* с точкой M) в некоторый момент времени в системах и * будем обозначать как r и r* .

Будем считать, что рассматриваемые системы отсчета движутся с такими скоростями, что расстояние между точками некоторого тела при взгляде из этих систем отсчета выглядит одинаково (исключаются релятивистские эффекты). Так, например, расстояние от одного угла стола до другого угла не изменится от того, смотрит ли на стол сидящий за ним наблюдатель (со своей системой отсчета) или наблюдатель, идущий мимо стола. Если мы обозначим радиус-вектор одного из углов стола в системе как r0, а радиус-вектор другого выбранного угла в той же системе отсчета как r, то модули векторов r r0 и r r0 совпадут.

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

r r0 = Q (r r0 ), (1.21) где Q – собственно ортогональный тензор поворота вида (1.16). Понять этот факт поможет рис. 1.5, на котором в качестве примера исследуемого объекта пунктирной линией изображен прямоугольник (стол). Если полностью совместить системы отсчета (их разворот и частичное смещение показаны на рис. 1.5, б), то векторы r r0 и r r0 будут отличаться на тот же поворот Q, что и оси систем координат X1OX2 и X1*O*X2* до их наложения.

Выражение (1.21) можно «обратить», умножив обе части равенства слева на тензор Q 1 :

–  –  –

Обычно принимается, что в некоторый момент времени, например, при t = 0, системы отсчета и * совпадают, то есть Q(0) = E. Точка с радиус-вектором r0 является так называемым полюсом, и обычно принимается, что она неподвижна относительно системы. Относительно времени t = t + a, как правило, принимается, что a = 0, то есть t = t .

Функция r0 r0 (t ) по сути задает трансляционное движение системы отсчета * относительно, а функция Q(t ) – вращательное движение * относительно .

В связи с произвольностью выбора системы отсчета возникают некоторые вопросы: правомочно ли использование физических законов и соответствующих математических соотношений, полученных в одной системе отсчета, при переходе к другой системе отсчета? Или они должны трансформироваться? Каким требованиям должны удовлетворять эти преобразования, какие параметры необходимы для подобной трансформации? Иными словами, возникает вопрос об объективности описания наблюдаемых процессов. Под объективностью мы будет понимать возможность предсказания результатов наблюдения в системе отсчета * по известным наблюдениям в системе и закону относительного движения систем. Этот закон, как мы видели, задается функциями r0 r0 (t ) и Q(t ) .

Никакой дополнительной информации для предсказания не должно требоваться. Например, расстояние между двумя точками твердого тела не должно зависеть от того, какой наблюдатель смотрит на тело, то есть вектор r r0 на рис. 1.5, в отличие от самого вектора r, не должен зависеть от наблюдателя (выбора системы отсчета). Для такого вектора мы получили закон преобразования (1.21), который можно обобщить на все не зависящие от выбора системы отсчета векторы. Итак, вектор a будем называть не зависящим от выбора системы отсчета, если при замене наблюдателя для него выполняется

a = Q a = a Q T. (1.23)

Примерами не зависящих от выбора системы отсчета векторов являются векторы базиса систем координат X1OX2 и X1*O*X2*, векторы сил взаимодействия между двумя телами (или частицами), а также векторы моментов сил .

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

–  –  –

Законы преобразования тензоров произвольного ранга определяются аналогичным образом. Но в общем случае ситуация сложнее, чем в соотношениях (1.23)–(1.25): во-первых, системы отсчета движутся одна относительно другой, во-вторых, тензорные объекты представляют собой физические характеристики изучаемых исследователем процессов, определяемые по известным правилам (причем одинаковым) каждым из наблюдателей. В этом случае законы, связывающие исследуемые тензорные характеристики, могут оказаться существенно сложнее приведенных соотношений, куда будут входить не только мгновенные значения Q(t), но и производные по времени от функций r0(t) и Q(t). Такая ситуация будет иметь место, например, при определении скоростей, ускорений, кинетической энергии частиц, для которых мы сейчас определим законы преобразования .

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

Скорость движения точки M определяется как производная по времени t векторной функции r(t):

–  –  –

Но мы знаем, что r = r0 + Q(t ) (r r0 ) = r0 + (r r0 ) Q т (t ), то есть скорость будет определяться следующим соотношением (при его получении используются правила дифференцирования суммы и произведения функций):

–  –  –

Появившееся произведение тензоров Q Q т называется тензором спина. Этот тензор описывает скорость вращения тела .

Для вычисления тензора спина найдем производную тензора поворота (1.16) вокруг оси m = e1, используя правило дифференцирования суммы и произведения:

–  –  –

Как видим, в общем случае скорость не является объективным вектором, поскольку закон ее преобразования (1.28) отличается от (1.23). Но для частных случаев движения наблюдателя * относительно скорость не будет зависеть от выбора системы отсчета. Это будет справедливо при v = 0 и Q(t ) = 0, то есть при отсутствии относительного движения систем отсчета (они могут быть сдвинуты одна относительно другой на постоянный вектор трансляции и развернуты на неизменный угол).

Далее нам понадобится выражение для вектора v v 0, которое получается из (1.28) аналогично получению выражения для r r0 :

–  –  –

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

Для получения закона преобразования ускорения точки w = dv(t ) / dt = = d 2r (t ) / dt 2 r (t ) надо найти производную по времени от выражения (1.28):

–  –  –

Итак, ускорение, так же как и скорость, в общем случае не является объективным вектором и зависит от выбора системы отсчета. Для частных случаев движения наблюдателей * относительно вектор ускорения может оказаться не зависящим от выбора системы отсчета в узком классе движения систем отсчета. Это будет выполняться при w = 0 и (t ) = 0, то есть при переходе в любую систему отсчета *, движущуюся относительно поступательно, равномерно и прямолинейно. Такие системы отсчета * называются галилеевыми, или системами отсчета, порождаемыми Ф .

Рассмотрим теперь преобразование уравнения движения материальной точки (второй закон Ньютона). Запишем его в системе отсчета, для которой примем, что уравнение движения приводится к виду mw = F, (1.30) где m – масса точки, w – ее ускорение, а F – действующая на нее сила .

Масса является характеристикой только рассматриваемой точки (или твердого тела) и не зависит от наблюдателя: m = m. Сила характеризует взаимодействие рассматриваемой точки или тела с другими точками или телами и зависит только от их взаимного расположения, которое, как мы обнаружили, не зависит от наблюдателя: F = Q(t ) F. А вот ускорение зависит от наблюдателя и преобразуется согласно (1.29). Выразим m, w и F через m*, w* и F* и подставим их в (1.30). Для силы имеем F = Q 1 (t ) F, тогда на первом шаге из (1.30) получим

–  –  –

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

Они появились вследствие перехода в движущуюся систему отсчета * .

Их называют силами инерции. Рассмотрим эти силы на примере задачи об учете собственного вращения Земли при описании движения тел по ее поверхности (скорости и ускорения тел мы обычно определяем именно относительно земной поверхности). Пусть система отсчета, не вращаясь, движется вместе с Землей (соответствующие оси в любых ее положениях параллельны самим себе), а система * вращается вместе с Землей относительно, центры систем отсчета совпадают с центром Земли. В качестве полюса также выберем этот центр. Тогда r0 = 0, v = 0, w = 0. Закон движения (1.36) примет вид

–  –  –

Рассмотрим пример. Пусть единичный вектор e, направленный из центра в сторону Северного полюса, задает ось вращения Земли, а постоянная угловая скорость вращения равна. Тогда тензор спина, как мы знаем, записывается в виде (t ) = e E = const .

Очевидно, что (t ) = 0, а 2 (t ) = 2 ( e E ) ( e E ) = 2 ( E e e ) .

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

–  –  –

Сила 2m e v называется силой Кориолиса. Она направлена перпендикулярно вектору скорости тела v* и оси вращения Земли e. Если тело движется по экватору, то скорость его, направленная по касательной к земной поверхности, будет перпендикулярна вектору e. Модуль силы Кориолиса при таком движении будет наибольшим. Эта же сила, действуя на частицы воды в текущих вдали от экватора реках, стремится отклонить направление их движения, что приводит к тому, что река постоянно подмывает один из своих берегов. В результате он становится более крутым, а другой берег, от которого река отходит, становится более пологим .

Вторая из оставшихся сил инерции m2 ( r ree ) называется центробежной силой. Вектор r можно представить в виде r = Rl, где R – радиус Земли, а l – единичный вектор, задающий направление от центра Земли к рассматриваемой материальной точке. Если тело находится на экваторе, то e l и центробежная сила максимальна по модулю и равна m2 R l (направлена перпендикулярно земной поверхности). На северном полюсе Земли r = Re, на южном – r = Re, то есть на любом из полюсов центробежная сила равна нулю. К этому же выводу можно прийти, заметив, что модуль r ree равен расстоянию от оси вращения Земли до рассматриваемой точки на ее поверхности. Для экватора это расстояние максимально, для полюсов оно равно нулю. Заметим, что сила инерции m r, не вошедшая в разобранный пример, называется эйлеровой и характеризует изменение угловой скорости вращения * относительно .

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

1. Как записывается тензор-проектор на плоскость с единичной нормалью n?

2. Чему равно произведение суммы диад m m + n n + p p на тензор поворота, если m, n, p – взаимно ортогональные единичные векторы?

3. Если при замене системы отсчета тензор T преобразуется по закону T = Q T Q т, то как выражается тензор T через T*?

4. Как преобразуется диада v v (v – скорость) при замене системы отсчета?

5. Запишите выражение для тензора спина, описывающего скорость поворота вектора скорости v плывущего вдоль экватора корабля относительно оси вращения Земли m .

1.3. БАЛАНСОВЫЕ УРАВНЕНИЯ И ЗАКОНЫ МЕХАНИКИ СПЛОШНЫХ СРЕД

Рассмотрим кратко на уровне основных идей смысл уравнений механики деформируемого континуума, справедливых для всех сред. Подробное обсуждение аксиом и физико-математических оснований, на которых базируется система получаемых уравнений, можно найти в работах [4–6, 8–13]. Договоримся, что начиная с этого раздела символ «» для обозначения тензорного умножения больше не будет использоваться. Для двух стоящих рядом векторов или тензоров будет приниматься, что между ними действует тензорное умножение, если не стоят знаки скалярного «» или векторного «» умножения .

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

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

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

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

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

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

Движение (деформация) континуума представляет собой изменение со временем взаимного расположения его частиц, то есть изменение конфигурации тела. Поэтому для выяснения, произошла ли деформация тела или насколько тело деформировалось, необходимо рассматривать две (или более) конфигурации тела. Для математического описания процесса деформирования необходимы меры, учитывающие различие конфигураций тела. Будем называть конфигурацию тела в момент времени t = t0 отсчетной или начальной конфигурацией и обозначать K0, а конфигурацию в произвольный момент времени t t0 – текущей или актуальной и обозначать Kt. Когда в курсе теоретической механики изучалось движение системы материальных точек, для возможности их различать точкам приписывались номера (или буквы), не меняющиеся в процессе движения .

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

Примем, что в момент t = t0 лагранжева система и пространственная система координат совпадают. Выберем в качестве пространственной системы координат декартову ортогональную систему, для пространственных координат примем обозначение xi, а для лагранжевых координат частиц – X i, i = 1,3. Координаты xi называют эйлеровыми координатами. Если все поля движущегося материального континуума рассматриваются как функции лагранжевых координат X i и времени t, то такой способ описания движения называется лагранжевым.

В этом случае закон движения записывается как xi = xi ( X 1, X 2, X 3, t ) или в векторной форме:

x = x( X, t ), (1.32) причем в силу принятого совпадения лагранжевых и эйлеровых координат при t = t0 <

–  –  –

При изучении движения твердых тел обычно принимаются следующие гипотезы: 1) представительный объем состоит из одних и тех же материальных частиц, 2) соседи частиц не меняются, 3) близкие материальные частицы при деформировании остаются близкими (в частности, не происходит нарушения сплошности в виде трещин, внедрений материала и других дефектов). Иначе говоря, две материальные точки не сольются в одну, а из одной не появятся две, три касательных вектора к материальным координатным линиям вектора в некоторой точке не окажутся результате деформации в одной плоскости и так далее. Математически первая гипотеза соответствует существованию взаимнооднозначного отображения (1.32) материального объема из отсчетной конфигурации на объем в текущей конфигурации. Вторая гипотеза означает непрерывность этого отображения (вспомните определение непрерывности на языке окрестностей). Третья гипотеза при выполнении двух первых соответствует требованию непрерывной дифференцируемости отображения при условии, что его якобиан нигде в рассматриваемой области не обращается в ноль и конечен. Все вместе эти гипотезы формулируются математически как требование диффеоморфизма отображения (1.32). Время при этом рассматривается как параметр. Напомним, что якобианом называют определитель матрицы Якоби частных производных, который для рассматриваемого случая имеет вид

–  –  –

В теории функций многих переменных доказывается теорема (принцип сохранения области) о том, что при выполнении условий непрерывной дифференцируемости функций xi = xi ( X 1, X 2, X 3, t ), i = 1,3 и xi / X j 0, xi / X j, i, j = 1,3 отображение сохраняет область, то есть в данном случае она остается трехмерной. В рассматриваемой области существует, и притом единственно, обратное отображение, якобиан которого также нигде не обращается в ноль и конечен. Из физического смысла отображения (1.32) следует также, что якобиан должен быть положительным, иначе последовательность переменных изменяется (меняются местами столбцы якобиана), что означает «выворачивание» материального объема. Итак, на якобиан отображения накладывается ограничение

–  –  –

также являющееся диффеоморфизмом .

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

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

Рассмотрим отсчетную и текущую конфигурации тела. Мы уже знаем, что положение точки P до деформации определяется вектором X, а после деформации – вектором x.

Бесконечно малый материальный отрезок PQ в недеформированном состоянии характеризуется вектором dX :

dX =| dX |, (1.36) где – единичный вектор, определяющий направление отрезка PQ.

Этот материальный отрезок в силу малости останется после деформирования также малым материальным отрезком P'Q', состоящим из тех же самых материальных частиц, но переместившихся в пространстве, и будет характеризоваться вектором dx :

dx =| dx |, (1.37) где ' – единичный вектор, определяющий направление отрезка P'Q' .

Относительным удлинением линейного элемента P'Q' в результате деформирования материала будем называть величину

–  –  –

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

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

–  –  –

+ (v ) = 0. (1.60) t Хотя это уравнение и представляет собой соотношение эйлерова подхода к описанию движения, но было выведено для фиксированного материального объема, то есть справедливо для пространственного объема, движение которого совпадает с движением среды (неподвижного относительно материала). Поэтому оно не является общим при пространственном описании движения, когда может возникнуть необходимость учесть перемещение объема относительно среды. Переход от лагранжевых к эйлеровым координатам, использованный выше, имеет место только при деформировании сред без изменения локальной топологии (например, для упругодеформируемых сред). В противном случае (жидкости, газы, пластические, сыпучие и другие неупруго деформируемые материалы) в качестве исходного соотношения необходимо брать уравнения обобщенного эйлерова подхода .

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

Тогда v (x, t ) = 0 и (x, t ) / t + v(x, t ) (x, t ) = 0 или d(x, t ) / dt = 0 .

Решением этого уравнения будет (x, t ) = (x,0) = 0 (i (x,0)) = 0 (i ), малопонятное в общем случае, но с учетом того, какой именно пространственный объем рассматривается, принимающее очевидный смысл .

Изменение в материальном объеме V(t ) некоторой физической вели

–  –  –

Проинтегрировав последнее уравнение по пространственному объему V (t ), движущемуся вместе с материалом, и используя теорему Остроградского-Гаусса, получим альтернативную формулировку интегрального балансового уравнения (для неподвижного относительно материала объема) t Vt ) (x, t )(x, t )dv = P (x, t )dv n (I (x, t ) + (x, t ) v (x, t )(x, t ))ds .

( V (t ) S (t )

–  –  –

При постоянной скорости V (x, t ) = V = const решением этого уравнения является поле (x, t ) = 0 (x Vt ), при произвольном поле V (x, t ) решение будет иметь вид (x, t ) = 1 (x f (t )), где f (t ) = V (t ) и 1 (x(i,0) f (0)) = 0 (i ) .

Будем далее рассматривать пространственное материальное описание и покажем, какие законы механики сплошных сред могут быть представлены в форме уравнения (1.63) d() / dt = I + P .

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

1) уравнение сохранения массы (уравнение неразрывности) может быть получено при выборе скалярной величины = 1 (масса, отнесенная к массе), I = 0 (не происходит изменения числа частиц в выбранном материальном объеме), P = 0 (отсутствуют источники массы):

+ v = 0, (1.69) где v – скорость центра масс элементарного объема; ( xi, t ) – плотность;

2) уравнение баланса количества движения (уравнение движения) может быть получено при = v (количество движения, отнесенное к массе), с учетом аксиом динамики сплошной среды, утверждающей, что причиной изменения количества движения материального объема являются действующие на него внешние силы. Внешние силы в механике принято разделять на объемные f e и поверхностные, которые согласно соотношению Коши записываются как t = n, где – тензор напряжений Коши,

n – нормаль к поверхности выбранного объема, то есть I =, P = f :

v = + f, (1.70)

f – массовая плотность внешних сил, действующих в объеме среды;

3) уравнения сохранения полной энергии (это такой же фундаментальный закон механики, как и уравнение баланса массы, уравнение баланса количества движения, уравнение баланса момента количества движения, а в термодинамике это первое начало без уточнения содержания потока) при = e (отнесенная к массе полная энергия материального объема), I = J e (поток энергии через поверхность малого объема), P = 0 (отсутствуют источники энергии):

e = J t, (1.71) где е – ее массовая плотность; Je – вектор потока полной энергии;

4) уравнение баланса энтропии (пока это понятие не определено, а рассматривается как некоторая физическая величина, необходимость введения которой будет обсуждаться ниже), записываемое формально в соответствии со структурой (1.63), s = J s + Ps, (1.72) где s – массовая плотность энтропии; Js – вектор потока энтропии; Ps – ее производство в элементарном объеме;

5) уравнение баланса тепла (является формой энергии наравне с другими составляющими полной энергии, такими как кинетическая или потенциальная энергия) q = J q + Pq (1.73) где q – массовая плотность количества тепла; Jq – поток тепла; Pq – объемный источник тепла (причиной его появления, помимо процессов, связанных с диссипацией энергии при интенсивном деформировании, может быть и некоторое излучение) .

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

–  –  –

или с использованием материальной производной относительно всей смеси и вводя барицентрическую скорость w n v n v (скорость относительно центра масс смеси) получим na n + n w n a n = (J n + n va n + n w na n ) + Pan + mna n. (1.80') Проведем анализ полученных соотношений в некоторых частных случаях. Если все рассматриваемые поля однородны, но изменяются во времени, то из (1.52) и (1.53) получаем

–  –  –

Рассмотрим частные случай применения записанных соотношений .

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

–  –  –

где скорость реакции m в общем случае является некоторой функцией концентрации. При рассмотрении не слишком больших интервалов времени течения реакции можно принять, что m = kc2 = k (1 c1 ), где k 0 соответствует распаду второй компоненты и ее превращению в первую .

В этом случае имеем c2 = c20 e kt /, c1 = 1 c20 e kt / = 1 (1 c10 )e kt /, где ci0 – начальные концентрации компонент. В более сложных случаях скорость реакции может быть пропорциональна уже появившемуся веществу, тогда m = k c1 (1 c1 ) и решение принимает вид

–  –  –

где f n – массовая сила, вызванная взаимодействием компонент смеси .

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

–  –  –

где da e + dU – удельная работа внешних сил над рассматриваемым объемом (включая его квазитвердое перемещение); dq** – дополнительный приток энергии за счет причин, отличных от работы макроскопических сил и теплообмена dq (например, за счет обмена массой или химических реакций), и переходя к материальной производной с использованием уравнения баланса тепла (1.73) и (1.81), получим

u = a e J q + Pq + q** v v. (1.90)

По теореме живых сил в действительном движении дифференциал кинетической энергии v v материального объема равен сумме элементарных работ внешних массовых, внутренних массовых, внешних поверхностных и внутренних поверхностных сил на перемещениях dr = vdt, то есть (получается с использованием уравнения движения (1.70))

–  –  –

где f vdt – удельная мощность внешних массовых сил, 1 ( v ) – удельная работа внешних поверхностных сил, а 1 : (v )Т – удельная работа внутренних поверхностных сил, то есть удельная работа внешних сил есть

–  –  –

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

s – энтропии:

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

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

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

Будем рассматривать среду без внутренних моментов, тензор напряжения Коши, входящий в записанные выше соотношения, симметричен, то есть : (v )T = : D, где D – тензор деформации скорости, тогда

u = J q + Pq + : D + q**. (1.95)

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

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

N u = q + : + µ n cn. (1.96) n =1 где µn – химический потенциал n-й компоненты. Здесь он вводится как коэффициент при скорости изменения концентрации. В таком виде химические потенциалы ввел Гиббс, чтобы учесть зависимость внутренней энергии от количества веществ. Из (1.96) следует, что µ n = u / cn, то есть, задав внутреннюю энергию, получим и вид химического потенциала. Часто бывает, что для конкретных процессов известен вид µn, тогда при задании внутренней энергии необходимо учесть знание о химическом потенциале .

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

Следовательно, во второе слагаемое записи внутренней энергии должна входить только обратимая (упругая) часть. Пусть тензор деформаций представим в виде разложения на упругую el и неупругую in составляющие = el + in. Тогда представление (1.96) примет вид N u = q + : el + µ n cn. (1.97) n =1

–  –  –

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

–  –  –

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

–  –  –

Заметим, что для записи линейных определяющих соотношений, связывающих обобщенные силы и переменные одного тензорного ранга, в свертке : D выделяют составляющие скалярной природы – следы тензоров и D и величины тензорной природы – девиаторы этих тензоров .

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

–  –  –

Учитывая связь концентраций c1 + c2 = 1, перейдем к одной концентрации c = c1, одному потоку J c = J c1, химическому потенциалу µ = µ1 µ 2 и скорости относительно центра масс w = w1 w 2 = v1 v 2. При этом (1.111) примет вид

–  –  –

Эти условия могут рассматриваться как частный случай условий на поверхностях разрыва в сплошных средах, подробный вывод которых можно найти в [4–5, 11–12] .

При решении задачи диффузии или теплопроводности граничные условия могут быть заданы несколькими способами. В теории диффузии различают граничные условия 1-, 2-, 3- и 4-го родов (граничные условия 4-го рода также называют условиями сопряжения). В зависимости от задачи может оказаться, что число граничных условий превышает число реальных границ. Например, в задаче для сферических или цилиндрических тел необходимо задать условия как на внешней поверхности (границе), так и в центре, где исследуемые поля не должны принимать бесконечных значений. Часто граничные условия задаются «на бесконечности» .

1. Граничные условия 1-го рода В граничных условиях 1-го рода задается распределение концентрации диффундирующего вещества по поверхности S тела, как функция координат и времени .

–  –  –

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

2. Граничные условия 2-го рода

В условиях 2-го рода задается распределение плотности потока концентрации примеси для каждой точки поверхности:

–  –  –

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

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

Jc S = 0 .

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

3. Граничные условия 3-го рода

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

–  –  –

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

Заметим, что на практике часто встречаются однородные граничные условия, когда правые части (1.124)–(1.127) равны нулю .

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

1c1 = 2c2, (1.128)

D1 n c1 S = D2 n c2 S, (1.129)

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

Начальные условия.

Задание начального условия для уравнения диффузии заключается в том, что для некоторого момента времени t = t0 (обычно полагают t0 = 0 ) должна быть известна функция c(r, t0 ) = c0 (r ) пространственных координат:

c(r,0) = c0 (r ). (1.130)

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

c(r,0) = 0 (r ). (1.131) t Причем начальные условия (1.130)–(1.131) должны быть согласованы с краевыми условиями .

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

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

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

–  –  –

1.4. ПРАКТИЧЕСКОЕ ВВЕДЕНИЕ В ПАКЕТ «MATHEMATICA»

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

В настоящее время существует масса специализированных программных комплексов, позволяющих довольно эффективно применять самые разные математические методы для решения возникающих при выполнении научных исследований проблем. Разумеется, никакой пакет программ никогда не сможет дать готовое решение задачи, представляющей научную ценность. Он является лишь инструментом в руках исследователя, упрощающим проведение некоторых промежуточных этапов, выполнение математических выкладок и позволяющим концентрировать внимание на принципиально важных вопросах об установлении причинно-следственных связей в исследуемых процессах, выяснении и проверке механизмов анализируемых явлений и в конечном счете их формализации на языке математики. В качестве примера такого инструмента выбран пакет «Mathematica» производителя Wolfram Research, наиболее приспособленного, по мнению авторов, к выполнению научных исследований. От версии к версии этого пакета синтаксис применения некоторых функций и значения опций, используемых во встроенных функциях, несколько меняется. Изложение данного раздела ориентировано на версию 7.0 .

Задача использования Пакет «Mathematica» предоставляет среду для одновременной отладки алгоритмов, используемых при построении дискретных и непрерывных математических моделей, графического представления и обработки результатов их работы, быстрого выполнения необходимых аналитических преобразований, численного исследования упрощенных вариантов дискретных и континуальных моделей .

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

Пакет «Mathematica» может использоваться на разных уровнях сложности решаемых задач или владения пакетом. Его можно применять как:

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

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

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

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

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

«Save» или «Save as» (сохранение) и «File» «Open» (открытие) файла .

Для вызова окна палитры используется кнопка «Palettes», после нажатия на которую в появившемся окне надо выбирать раздел «Other», а в нем – «Basic Math Input». Другая кнопка, которую нужно регулярно использовать, – это кнопка «Help», позволяющая получить справку по многим вопросам о работе в пакете. В частности, рекомендуется познакомиться с «Help» «Virtual Book». В дальнейшем интерфейс пакета и содержимое меню обсуждаться практически не будут. Любой пользователь персонального компьютера с начальным опытом может разобраться в них самостоятельно. Более подробную информацию о работе в пакете можно найти, например, в [3] или воспользоваться встроенной книгой помощи .

При вводе любого символа в рабочем окне открывается так называемая ячейка (отмечена квадратной скобкой по правому краю рабочего окна), в котором и появляются вводимые символы. Например, введем «2+9». Для выполнения команды необходимо одновременно нажать комбинацию клавиш «Shift»+«Enter» или «Enter» на вспомогательной (цифровой) части клавиатуры. Результат вычисления появится в следующей ячейке. Обратим внимание на то, что первое вычисление в любой сессии пакета занимает большее время, чем последующие вычисления выражений аналогичной сложности, поскольку при этом подгружается ядро пакета, в котором хранятся результаты всех вычислений текущей сессии .

Ячейки ввода автоматически нумеруются и обозначаются «In[n]», где n – порядковый номер ячейки, выходные ячейки обозначаются «Out[n]» .

Даже если ячейки удаляются из рабочего окна, результаты их вычислений сохраняются в ядре и результат их вычислений можно увидеть на экране, выполнив команду «Out[n]» с нужным номером (например, можно набрать «Out[1]» и нажать «Shift»+«Enter»). Итак, «Mathematica» относится к языкам-интерпретаторам, все команды обрабатываются в ядре пакета в последовательности их исполнения (при нажатии комбинации клавиш «Shift»+«Enter»), и результат каждого вычисления не зависит от положения соответствующей ячейки в рабочем окне, а только от порядка его выполнения относительно других команд .

Перейдем к описанию работы с математическими выражениями .

Арифметика: сложение «+», вычитание «–», умножение «», деление «/», возведение в степень «^». Отметим, что операцию «» можно заменять пробелом (результат вычисления в пакете выражения «4 5» отличается от вычисления «45»), а в случае когда числовой множитель стоит перед символьным выражением, никакой знак ставить не нужно. Например, выражения «5x», «5 x» и «5x» эквивалентны, причем рекомендуется использовать последнюю запись для облегчения чтения кода. Выражения «5x» и «x5» не эквивалентны: первое есть произведение числа 5 и переменной x, а второе означает независимую переменную x5, для перемножения в последнем случае необходимо вставить пробел «x 5». Выполните для сравнения команды «3x + 5x» и «x3 + x5». Часто для улучшения читаемости вводимых выражений используется традиционная математическая нотация. Возведение в степень при этом задается не как «x^3», а с помощью комбинации клавиш «Ctrl»+«^» после ввода «x». Операция деления представляется в виде дроби с помощью комбинации «Ctrl»+«/». Квадратный корень вводится с помощью комбинации клавиш «Ctrl»+«@» .

Обратите внимание, что пакет «Mathematica» выделяет введенные символы различными цветами – числа отображаются черным цветом, переменные «x», «x2», «x5» – синим. Это означает, что для пакета величины «x», «x2», «x5» (выделенные синим цветом) не известны или не определены, а числа являются известными величинами .

Для большинства операций в пакете возможны различные способы их реализации. Например, сложение можно задавать, как было показано выше, с помощью «+», а можно с помощью функции «Plus[]» – выполните «Plus[2,9,–1,x]». Аналогично умножение можно представить иначе, как «Times[3,x,y]», возведение в степень как «Power[x,3]» и так далее .

Пакет «Mathematica» позволяет работать с логическими выражениями (относительно которых можно сказать, что они истинны «True» или ложны «False»). Для сравнения выражений применяются операторы «==» – равны ли выражения, «! =» – не равны, «» – меньше, «» – больше, «=» – меньше или равно, «=» – больше или равно. Сравните результаты выполнения команд «53», «7==7», «41», «12=19». Выясните с помощью логических операций пакета, что больше, 13 12 или 7 6, 27 + 6 + 1 или 48, 3 или 21/2, расположите в порядке возрастания 5 2, 4 3, 3 6. Для манипулирования логическими выражениями используются операторы логического «И» – «&&» и логического «ИЛИ» – «||». Отметим, что последовательность символов «! =» автоматически преобразуется пакетом в знак «», «=» – в знак «», «=» – в знак «» и так далее .

В пакете реализована работа с точными и приближенными типами данных. К первым относятся числа «», « », « », « 3 » и так далее .

Заметим, что «» и « » отображаются черным цветом, хотя являются символами. Это означает, что пакету известны их значения .

После выполнения команды по вводу точных выражений они могут быть упрощены сокращением дроби (выполните команду «14/20»), вычислением точного значения функции, но приведенные выше выражения не преобразуются. Для перехода к их приближенному значению в числовых выражениях необходимо ставить десятичную точку, например, « 3.0 » или просто « 3. », но что делать с «»? Можно умножить или разделить на «1.0», можно добавить «0.0», но существует более удобный способ – использовать встроенную функцию пакета. Все встроенные функции начинаются с заглавных букв, а аргумент указывается в квадратных скобках (как, например, уже встречавшаяся функция «Out[1]») .

Для приближенного вычисления значения выражений используется функция «N[]» (от английского numerical), например, «N[]». Многие функции пакета имеют дополнительные необязательные аргументы и могут вызываться с разным числом аргументов. Функция «N[]» имеет второй необязательный аргумент, задающий количество цифр для вывода приближенного значения: выполните команду «N[,150]», в результате в ячейке вывода появится 150 значащих цифр числа .

Для ввода греческих букв и специальных символов с клавиатуры (часто это оказывается быстрее, чем делать вставку соответствующего знака из открытого окна палитры) используется следующая последовательность действий: нажимается клавиша «Esc» (на экране появится символ « »), вводится соответствующая буква латинского алфавита, например, «p» и снова нажимается «Esc» (последовательность « p » автоматически преобразуется в «»). Для ввода экспоненты надо набрать «Esc», «ee», «Esc» (если набрать только один раз «e», то результатом будет греческая буква «»), для ввода мнимой единицы – «Esc», «ii», «Esc» (на экране она отображается как « ») .

В арифметических выражениях используются такие типы данных, как Integer, Rational, Real. Для выяснения типа выражения применяется функция «Head[expr]» (заголовок). Выясните, какое значение она возвращает в следующих случаях: Head[5], Head[7/3], Head[3.14159], Head[]. «Mathematica» дает каждому выражению и каждой его части свой тип, например, выражение «2+» имеет заголовок «Plus», выражение «Sin[+3]» – заголовок «Sin» .

В пакете доступны стандартные математические функции: тригонометрические – «Sin[]», «Cos[]», «Tan[]»,»Cot[]», обратные тригонометрические функции –»ArcSin[]», «ArcCos[]», «ArcTan[]», «ArcCot[]», извлечение арифметического квадратного корня – «Sqrt[]», экспонента – «Exp[]», логарифм по основанию b – Log[b, ]» (натуральный логарифм «Log[]»). Вычислите, например, «Sin[ ]», «Log[5e]» .

Пакет «Mathematica» позволяет работать и с комплексными числами .

Для ввода мнимой единицы, как упоминалось выше, используется последовательность нажатия клавиш «Esc», «ii», «Esc». Найдите квадрат мнимой единицы « ». Вычислите также приближенное значение комплексного числа вида. Функции, позволяющие находить действительную и мнимую части комплексного числа, записываются соответственно как «Re[]» и «Im[]». Примените их к результату предыдущего вычисления .

«Mathematica» предоставляет массу средств рисования графиков функций. Универсальная функция, которая для этого используется, – это функция «Plot[]». Она имеет два обязательных аргумента: саму функциональную зависимость, график которой надо отобразить, и интервал изменения аргумента, например, для рисования параболы необходимо выполнить команду «Plot[x2,{x,–2,3}]». Второй аргумент {x,–2,3} имеет следующий смысл – x есть переменная, которая меняется от –2 до 3 .

Постройте график функции y = sin( x) на отрезке [–; 2] .

Для изображения на одном рисунке графиков нескольких функций первый аргумент представляется в виде {x2,x3,Exp[2x]} (такая конструкция называется списком), то есть в фигурных скобках перечисляются исследуемые функции: «Plot[{x2,x3,Exp[2x]},{x,–2,3}]». Результат, полученный в «Mathematica 7.0», будет выглядеть, как на рис. 1.6, а. Обратите внимание, что графики перечисленных функций изображаются различными цветами, выбираемыми пакетом автоматически .

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

Для того чтобы вывести все возможные точки графика функции, этой опции надо придать значение «All»: Plot[{x2,x3,Exp[2x]},{x,–2,3}, PlotRange–All]» (пакет автоматически преобразует последовательность знаков «–» в стрелку «»). Для задания нужного размера рисунка по оси ординат опции «PlotRange» надо передать диапазон значений, например (рис. 1.6, б), “Plot[{x2,x3,Exp[2x]},{x,–2,3},PlotRange{–10,15}]” .

Рис. 1.6. Результаты выполнения команд: а – “Plot[{x2,x3,Exp[2x]},{x,–2,3}]”;

б – “Plot[{x2,x3,Exp[2x]},{x,–2,3},PlotRange{–10,15}]” Можно задать и окно вывода (например, для более детального изображения графиков функций вблизи некоторой точки – точки пересечения или особой точки функции и так далее) с помощью структуры «PlotRange {{xmin,xmax},{ymin,ymax}}». При этом можно предоставить пакету выбор размеров по одному из направлений с помощью слова «Automatic”, например:

“PlotRange{{xmin,xmax},Automatic}” или слова «All»: «PlotRange {{xmin,xmax},{0,All}}» .

Чтобы сохранить получаемые рисунки, надо навести указатель (стрелку) мыши на рисунок и нажать на правую кнопку. В появившемся меню нужно выбрать строку «Save Graphic As…» и нажать на левую кнопку мыши. Далее необходимо выбрать место, в котором будет сохранен рисунок, и ввести его название. Пакет предлагает для сохранения рисунка большое количество форматов. Для последующего размещения рисунков в текстовых документах удобно использовать типы файлов «GIF», «JPEG», «TIFF» .

Для задания свойств линии, которой рисуется график, применяется опция “PlotStyle”. Этой опции может быть передано несколько значений. Чтобы нарисовать утолщенную линию, необходимо передать этой опции значение “Thick”, например, “PlotStyleThick”. Для задания другой толщины линии используется функция “Thickness[]”, например, “PlotStyle Thickness[0.01]”. Аргумент задает долю, которую составляет толщина линии от высоты рисунка. Для изображения графика некоторой математической зависимости пунктирной линией опции “PlotStyle” передается значение “Dashed” или используется функция “Dashing[{r1,r2,r3,…}]”, где сегменты с длинами, соответствующими значениям ri (доли, которые составляют длины сегментов пунктирной линии от высоты рисунка), повторяются циклически. Например, можно задать “PlotStyleDashing[{0.08,0.02}]” или “PlotStyle Dashing[0.05]”. Для того чтобы одновременно задать и толщину, и пунктир, используется структура “PlotStyle{Thickness[0.01], Dashing[0.05]}”. Здесь порядок функций в фигурных скобках значения не имеет. Если при выводе нескольких кривых на одном рисунке требуется задавать свой стиль для каждой кривой, то необходимые свойства перечисляются в фигурных скобках: “PlotStyle{{Styles-1},{Styles-2},…}”, например (рис. 1.7), “Plot[{x2,x3,Exp[2x]}, {x,–2,3}, PlotStyle{{Thickness[0.01], Dashing[0.05]}, {Thickness[0.0125]}, {Dashing[{0.08,0.02}]}}]” .

–  –  –

Для удобства вывода, помимо толщины и вида пунктира, в стиле графика можно указать его цвет. Самый простой вариант – использовать названия цветов (“Black”, “Gray”, “Blue”, “Brown”, “Cyan”, “Green”, “Magenta”, “Orange”, “Pink”, “Purple”, “Red”, “Yellow”). Эти цвета можно сделать темнее с помощью функции “Darker[]”, например, “Darker[Green]” или светлее с помощью функции “Lighter[]”, например, “Lighter[Blue]” .

Также можно задать цвет с помощью функции “RGBColor[r1,r2,r3]”, в которой величины r1, r2, r3 задают доли красного, зеленого и синего цветов соответственно. Например, параболу синего цвета можно построить с помощью команды “Plot[x2,{x,–2,3},PlotStyle{Thickness[0.02], RGBColor[0,0,1]}]” .

Задавать промежуточные оттенки с помощью этой функции не всегда удобно. Для этого можно использовать названия оттенков, список которых можно получить в “Mathematica 7.0” с помощью функции “ColorData[]” .

Эта функция позволяет получить обширный набор цветовых схем и отдельных цветов. Рассмотрим некоторые примеры ее использования. Выполнение команды “ColorData[]” приводит к появлению списка ее возможных аргументов “{Gradients, Indexed, Named, Physical}”, соответствующих наборам цветовых схем: “Gradients” дает названия схем, цвет в которых плавно изменяется от одного до другого при непрерывном изменении от 0 до 1 параметра схемы. Например, выполнив команду “ColorData["Gradients"]”, можно получить список названий предопределенных цветовых схем. Далее эти названия используются для задания конкретного цвета. Выберем для примера название “TemperatureMap”. Доступ к цветам этой схемы дает команда “ColorData["TemperatureMap"]” (рис. 1.11). Конкретный цвет из этой схемы получается заданием значения параметра. Например, выполните “Plot[x2,{x,–2,3},PlotStyle{Thick,ColorData["TemperatureMap",0.2]}]” .

Рис. 1.8. Результат выполнения команды “ColorData["TemperatureMap"]” (приведен вид цветовой схемы, крайние оттенки цвета которой соответствуют значениям параметров 0 и 1) Поименованные цветовые схемы можно получить, выполнив команду “ColorData["Named"]”, которая приводит к появлению на экране списка имен схем: {Atoms, GeologicAges, HTML, Legacy, WebSafe}. Например, схема “Atoms” содержит цвета, соответствующие атомам из периодической таблицы Д.И. Менделеева. Список названий цветов этой схемы можно получить с помощью команды “ColorData["Atoms","Range"]”.

Используются они, например, таким образом (взят «цвет меди» – Cu):

Plot[x2,{x,–2,3},PlotStyle{Thick,ColorData["Atoms","Cu"]}]” .

Чтобы увидеть все предопределенные в пакете «цвета атомов», необходимо выполнить команду “ColorData["Atoms", "Image"]” (результат приведен на рис. 1.9). Набрав команду “ColorData["Legacy","Range"]” и нажав “Shift + Enter”, получим доступные названия оттенков цветов.

Теперь эти названия могут использоваться для рисования графиков, например:

“Plot[x2,{x,–2,3},PlotStyle{Thick,ColorData["Legacy", "SteelBlue"]}]” .

–  –  –

Другой полезной опцией является задание названий осей координат .

Для этого применяется опция “AxesLabel”, например:

“Plot[x2,{x,–2,3},PlotStyle {Thick,Darker[Green]},AxesLabel{"X", "Y"}]” .

Для задания названия рисунка используется опция “PlotLabel”, например:

“Plot[x2,{x,–2,3},PlotStyle{Thick,Darker[Green]}, AxesLabel{"X", "Y"},PlotLabel"Fig.1"]” .

Для того чтобы отменить изображение меток на осях координат, задать свои метки в нужных местах на осях применяется опция “Ticks”. Опция “Ticks False” убирает все метки на обеих осях, Ticks{Automatic, False}” отменяет метки на оси ординат и оставляет автоматическое расположение меток на оси абсцисс.

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

–  –  –

Функция “Plot[]” содержит ряд других опций, информацию о которых можно найти в пакете. Для получения основной информации о функции “Plot[]” надо выполнить команду “? Plot”, для получения более подробной информации необходимо выполнить “?? Plot” либо подвести курсор к слову “Plot” и нажать клавишу F1. Аналогично можно посмотреть справку о любой другой функции пакета или используемого символа .

На основе функции “Plot[]” сформировано несколько других функций пакета для графического отображения информации. Полезной модификацией является функция “ParametricPlot[]”, которая строит на плоскости график функции, заданной параметрически. Структуру этой функции можно представить как “ParametricPlot[{x[t],y[t]}, {t,t1,t2}]”. Здесь в качестве дополнительных аргументов можно использовать те же опции, которые приводились выше при обсуждении функции “Plot[]”. Например, построим график “ParametricPlot[{2Cos[t],Sin[t]},{t,0, 7}, PlotStyleThick]” .

Другая модификация функции “Plot[]” – это функция “DensityPlot[]”, которая отображает на координатной плоскости значения функции двух переменных, используя различные цветовые оттенки. Использование этой функции аналогично использованию функции “Plot[]”. Например, выполним “DensityPlot[Exp[–x2 – y2], {x, –2, 2}, {y, –2, 2}]” .

В результате получим плоскую область, точки которой закрашены в соответствии с выбранной по умолчанию цветовой схемой, отображаемые оттенки будут соответствовать значениям заданной функции “Exp[–x2 – y2]”. Опция “PlotStyle” для данной функции не определена .

Для выбора другой цветовой схемы применяется опция “ColorFunction”, которой передается название этой схемы, например, “ColorFunction "TemperatureMap"” (рис.

1.11, а):

DensityPlot[Exp[–x2 – y2], {x, –2, 2}, {y, –2, 2}, ColorFunction"TemperatureMap"]” .

Часто использование функции “DensityPlot[]” с заданными по умолчанию значениями дополнительных аргументов приводит к получению грубого изображения. Улучшить рисунок помогает опция “PlotPoints”, задающая число точек, используемых для построения изображения, например (рис. 1.11, б), DensityPlot[Exp[–x2 – y2], {x, –2, 2}, {y, –2, 2}, ColorFunction"TemperatureMap", PlotPoints100, FrameTicks False, Mesh2]” .

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

а б

Рис. 1.11. Пример работы функции “DensityPlot[]” с параметрами:

а – заданными по умолчанию; б – заданными пользователем Функция “DensityPlot[]” имеет опцию, позволяющую выводить на экран область, определяемую некоторыми условиями: “RegionFunction Function[{x, y}, 1 x2 + y2 5]”. В данном примере область вывода определяется неравенством 1 x2 + y2 5. Цвет границы отображаемой области задается опцией “BoundaryStyle Red”. Результат приведен на рис. 1.12 .

–  –  –

Другую модификацию функций, используемых для графического представления результатов, представляет добавление «окончания» “3D”:

“Plot3D[]”, “ParametricPlot3D[]” и другие. Например, выполните команду (рис. 1.13, а) <

–  –  –

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

Для построения трехмерных кривых, заданных параметрически, применяется функция “ParametricPlot3D”. Структуру этой функции можно представить как “ParametricPlot3D[{x[t],y[t],z[t]}, {t,t1,t2}]”.

Например, построим график винтовой линии на цилиндре с эллипсом в основании:

“ParametricPlot3D[{2Cos[t],Sin[t],0.3t}, {t,0,7}, PlotStyleThick, TicksFalse]” .

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

Далее они будут использоваться в различных примерах, здесь же рассмотрим, как построить изображение сферы. Сфера задается функцией “Sphere[{x1,x2,x3},r]”, где тройка чисел {x1,x2,x3} задает положение центра сферы в декартовой ортонормированной системе координат, а r – радиус сферы. Для построения изображения сферы надо применить функцию “Graphics3D[]”, например, “Graphics3D[Sphere[{1,0,0}, 0.5]]” или, задавая цвет поверхности, “Graphics3D[{Brown, Sphere[{1,0,0}, 0.5]}]” .

Также можно построить изображение нескольких сфер:

“Graphics3D[{{Brown, Sphere[{1,0,0},0.5], Sphere[{2,0,0},0.5]}, {Darker[Green], Sphere[{3,0,0}, 0.5]}}]” .

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

H2O. Две сферы будут изображать атомы водорода соответствующего цвета:

{ColorData["Atoms","H"], Sphere[{0, 0, 0}, 0.7], Sphere[{1.4, 0, 0}, 0.7]} и одна сфера – атом кислорода O (если радиус не указан, то он равен 1):

{ColorData["Atoms", "O"], Sphere[{0.7, 0, 0.7}]} .

Для одновременного изображения всех указанных атомов выполните команду (результат выполнения приведен на рис. 1.14, а) “Graphics3D[{{ColorData["Atoms","H"], Sphere[{0, 0, 0}, 0.7], Sphere[{1.4, 0, 0}, 0.7]}, {ColorData["Atoms", "O"], Sphere[{0.7, 0, 0.7}]}}]” .

Отметим, что внутренние фигурные скобки здесь можно не ставить, а функции, задающие цвет, действуют на все перечисленные после них объекты. Для управления изображением можно добавить функцию “Specularity[White, 50]”, которая задает степень зеркальности (блеска) изображаемой поверхности, а также опцию “Lighting”, которая задает способ освещенности объекта (положение источника света, его цвет, размер, форму). Значение опции “Lighting"Neutral"” соответствует источнику белого цвета. Более подробную информацию об этих дополнительных параметрах можно найти в книге помощи, встроенной в пакет. Рассмотрим пример их использования для молекулы H2O (пример из встроенной книги документации пакета «Mathematica 7.0», результат выполнения команды приведен на рис.

1.14, б):

“Graphics3D[{Specularity[White, 50], ColorData["Atoms", "H"], Sphere[{0, 0, 0}, 0.7], Sphere[{1.4, 0, 0}, 0.7], ColorData["Atoms", "O"], Sphere[{.7, 0,.7}]}, Lighting"Neutral"]” .

а б

Рис. 1.14. Пример работы функции “ Graphics3D[]” с различными параметрами:

а – стандартные значения; б – добавлен блеск и освещенность Пакет «Mathematica 7.0» предоставляет доступ к данным о химических веществах, которые подгружаются через Интернет и могут быть отображены в различном виде – как химические формулы, как схемы соединений, как трехмерные изображения молекул. Доступны химические и физические свойства этих веществ. Эта информация может быть получена с помощью функции “ChemicalData[]” .

–  –  –

в Рис. 1.15. Изображение молекул: а – структурная формула Fe5O12Y3;

б – структурная формула C8H10N4O2; в – трехмерное изображение C8H10N4O2 Чтобы увидеть разделение доступных веществ по классам, применяется команда “ChemicalData["Classes"]”. Далее имя класса можно использовать для получения списка веществ этого класса, например, выполнив команду “ChemicalData["Nanomaterials"]”, получим список химических соединений, молекулы которых используются в нанотехнологиях как наночастицы. Используя их названия, можно узнать их химическую формулу, ее структурное представление, физические свойства. Например, для получения химической формулы надо выполнить команду “ChemicalData["YttriumIronOxide", "HillFormulaDisplay"]”, в результате получим Fe5O12Y3. Структурная формула этой молекулы имеет вид, получаемый с помощью функции “ChemicalData["YttriumIronOxide"]” (рис. 1.15, а). Аналогично можно получить формулу молекулы кофеина:

“ChemicalData["Caffeine", "HillFormulaDisplay"]”, как в виде C8H10N4O2 или структурной формулы (рис. 1.15, б) с помощью функции “ChemicalData["Caffeine"]”, так и в виде трехмерного изображения молекулы (рис. 1.15, в) с помощью функции с дополнительным параметром “ChemicalData["Caffeine", "MoleculePlot"]” .

Переменные и функции. Для создания символьной переменной выполняется команда, например, “a=5”, которая автоматически создает в ядре пакета переменную a с текущим значением 5, и во всех выражениях, содержащих a, которые будут выполнены позже, вместо a будет подставлено число “Integer[5]”. То есть появляется глобальная переменная с заданным значением. Отметим, что после выполнения команды присваивания переменная a стала изображаться черным цветом, а не синим, как до выполнения команды. Для того чтобы очистить переменную a, например, для ее использования в символьных преобразованиях, выполняется команда “Clear[a]” или “a=.”. После этого цвет переменной a на экране снова становится синим .

–  –  –

В «Mathematica» может использоваться несколько форм выполнения некоторой операции. В частности, для присваивания может использоваться функциональная форма “Set[]” – “Set[a,19]” .

Если придать какой-либо переменной новое значение, а затем выполнить записанные ранее выражения с этой переменной, то будут получены другие ответы. Часто бывает необходимо изменить значение переменной с использованием ее предыдущего значения, например “a = a + 7”. Для упрощения аналогичных команд в пакете реализованы операции: “+=”, “–=”, “*=”, “/=”. Если требуется изменять значение переменной на единицу (например, при реализации циклов), то можно использовать удобные короткие формы: “++” и “––”, которые используются как “i++”, “i––”, “++i”, “––i” .

Пример: a = 7 --- 7 ++a --- 8 a++ --- 8 a --- 9 Выполните приведенную последовательность команд и объясните, как получаются такие результаты .

«Mathematica» позволяет вычислять выражении без присваивания переменной некоторого значения. Для этого используются так называемые подстановки. При этом в конец выражения, которое требуется вычислить, дописывается последовательность символов вида “/.x5”, в которой вместо числа 5 может задаваться любое значение переменной x, в том числе и другая переменная. Набор символов “x5”, идущий после “/.”, называется правилом замены переменной или правилом подстановки. При преобразованиях выражений последовательно можно применять любое число правил замены. Для ввода стрелки “” используется последовательность символов “–” и “” .

Если в вычислениях возникает необходимость многократного использования некоторой последовательности операций – это хороший повод создать собственную функцию. Синтаксис задания функции требует указания ее аргумента, который определяется знаком подчеркивания после имени переменной, и задания тела функции, например: “myFunc[x_]: =x2+x+1”. Имя собственной функции пользователя может начинаться с любой буквы, не обязательно заглавной. Подчеркивание указывает, что переменная с именем x является локальной, и при вызове функции на ее месте может оказаться любое выражение или число. При вводе приведенной команды видно, что пакет «Mathematica» автоматически изображает переменные курсивом зеленого цвета. Операция “:=” называется отложенным присваиванием и означает, что правая часть будет вычисляться только при вызове функции, то есть тогда, когда функции f[x] будет передано конкретное выражение для x. А обычное присваивание “=” по аналогии можно назвать немедленным присваиванием. Итак, в пакете мы встретили уже три различных «знака равенства»: “==”, “=”, “: =” .

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

Рассмотрим обычные ошибки, допускаемые при задании функций .

–  –  –

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

Пример: x = 5; bad2[x_] = x2 + x – 1;

bad[y] --- 29 Здесь при инициализации функции произошло вычисление ее правой части x2 + x – 1 при заданном значении x, то есть была введена функция с постоянным значением 29. В таком случае использование немедленного присваивания “=” приводит к ошибке .

Рассмотрим другой пример. Для нахождения производной функции по некоторой переменной используется встроенная функция “D[]”, например “D[x2,x]”. Создадим функцию, которая будет возвращать значение производной выражения x2 + x – 1 при различных значениях переменной x .

Пример: bad3[x_]: = D[x2 + x – 1, x] bad3[5] --- 5 is not a valid variable .

Произошло не то, что ожидалось.

Из-за отложенного присваивания в правую часть тела функции D[x2 + x – 1, x] вместо x было подставлено значение 5 – D[29, 5], но производной не существует! В этом примере при задании функции необходимо использовать немедленное присваивание, чтобы сначала была найдена производная, а потом функция возвращала ее значения при различных значениях переменной x:

x =.; g[x_] = D[x2 + x – 1, x];

g[5] --- 11 Помимо функций, задаваемых выше с помощью явного указания имен функций и переменных и поэтому называемых поименованными функциями, пакет «Mathematica» позволяет работать с так называемыми анонимными функциями. Они удобны при однократном применении некоторого выражения. Вместо имени функции при использовании анонимной функции указывается конструкция, оканчивающаяся знаком “&” .

Поскольку эта конструкция заменяет имя функции, то после нее можно ставить квадратные скобки с указанием аргумента (или нескольких аргументов). Конструкция, о которой идет речь представляет собой тело функции, в котором вместо имени переменной используется символ “#” для функции одного аргумента или набор вида “#1”, “#2” и так далее для функции нескольких аргументов. Например, для вычисления значения выражения x2 + x – 1 может быть создана анонимная функция “ (#2 + # – 1) &”. Ее применение выглядит как “ (#2 + # – 1) &[x]”. Если в теле функции используются логические операции, то результатом будет оценка истинности получаемого высказывания, например, “ (#12 #22) &” будет сравнивать квадраты двух передаваемых функций аргументов. Для вычисления производной с помощью анонимной функции необходимо заставить пакет предварительно выполнить символьное преобразование производной. Сделать это можно с помощью функции “Evaluate[]”. Например, в ситуации, аналогичной рассмотренной в приведенном выше примере, запись “Evaluate[D[#2 + # – 1, #]]&” дает корректный результат .

Списки. Объекты этого типа уже встречались выше, например, когда задавались свойства линии при рисовании графика функции. Список представляет собой перечень элементов в фигурных скобках. Причем элементом списка может быть другой список, элементы которого называются элементами второго уровня, а список – многоуровневым. Заголовок этого объекта называется “List”, и задавать список можно функцией “List[]”, перечисляя в качестве ее аргументов все элементы списка “List[a,b,c,3,–1,x]” .

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

Полученный список сохраним в переменной atomChain:

“atomChain={{1,0,0},{2,0,0},{3,0,0},{4,0,0}}” .

Функция “Length[]” возвращает длину списка, что в нашем случае соответствует количеству атомов:

“ nAtoms=Length[atomChain]” .

Для обращения к элементу списка используется имя списка и номер элемента, который указывается в двойных квадратных скобках, например “atomChain[[3]]”. Если в двойных квадратных скобках указано отрицательное число, например “atomChain[[–1]]”, то номер элемента будет отсчитываться с конца списка (в приведенном примере будет взят последний элемент). Для выбора случайного атома необходимо получить случайное число от 1 до nAtoms. Случайные числа в пакете генерирует функция “Random[]”. Ее первый аргумент указывает тип возвращаемого числа, второй – интервал, из которого выбирается число.

В нашем случае нужны целые случайные числа:

“Random[Integer,{1, nAtoms}]” .

“atomChain[[Random[Integer,{1, nAtoms}] ]]” .

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

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

Ранее список задавался перечислением его элементов. В «Mathematica» существуют функции для генерации списков. Одним из наиболее универсальных способов создания списков является использование функции “Table[expr,iterator]”. Первым аргументом стоит выражение expr, которое в результате повторных исполнений будет задавать элементы списка, второй аргумент – итератор, который определяет характер повторения expr. Итератор всегда задается в фигурных скобках и может иметь следующий вид: “{n}” – выражение expr будет повторено n раз (если оно содержит функцию, генерирующую случайные числа, то результаты исполнения будут различными), “{i,n}” – итератор содержит переменную итерирования i, которая принимает значения от 1 до n, “{i,n1,n2}” – переменная i принимает значения от n1 до n2, “{i,n1,n2,st}” – переменная i принимает значения от n1 до n2 с шагом st. Для задания цепочки атомов заданного их количества n выполним “Table[{i,0,0},{i,n}]”, предварительно задав значение n, например, “n=20”. Итераторов может быть несколько – для одного итератора получаемый список соответствует некоторому условному вектору, для двух итераторов – плоской матрице, для трех итераторов будет получаться трехмерный массив и так далее .

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

Функция “Join[]” соединяет списки, не выполняя упрощений. Функция “Union[]” объединяет любое количество списков в теоретико-множественном смысле, удаляя повторяющиеся элементы и сортируя оставшиеся элементы списка в алфавитном порядке (применение к одному списку приводит к удалению повторов и сортировке его элементов). Функция, которая задает теоретико-множественное пересечение списков, имеет название и синтаксис применения вида “Complement[list,list1,lit2,…]”, где из списка list удаляются все элементы, встречающиеся в списках list1, lit2 и так далее .

Функции “Append[list,el]” и “Prepend[list,el]” добавляют элемент el к списку list. Первая добавляет новый элемент в конец списка, вторая – в начало. Функция “Delete[list,k]” удаляет элемент с номером k из списка list. Если число k является отрицательным, то отсчет номера удаляемого элемента ведется с конца списка. Функция “Flatten[list]” удаляет все внутренние фигурные скобки, то есть делает многоуровневый список одноуровневым. Дополнительный аргумент задает номер уровня, до которого необходимо проводить упрощение списка: “Flatten[list,L]”.

Например, генерируя с помощью функции “Table[a{i,j,k},{i,n},{j,n},{k,n}]” набор атомов с межатомным расстоянием a, соответствующий простой кубической решетке, получим список сложной структуры:

“n=2; atomList = Table[a{i,j,k},{i,n},{j,n},{k,n}]” --a,a,a},{a,a,2a}},{{a,2a,a},{a,2a,2a}}}, {{{2a,a,a},{2a,a,2a}},{{2a,2a,a},{2a,2a,2a}}}}” .

Применение функции “Flatten[atomList]” приводит к такому результату:

“{a,a,a,a,a,2a,a,2a,a,a,2a,2a,2a,a,a,2a,a,2a,2a,2a,a,2a,2a,2a}”, когда структура списка полностью утрачивается. Применение функции с дополнительным аргументом “Flatten[atomList,1]” приводит к списку “{{{a,a,a},{a,a,2a}},{{a,2a,a},{a,2a,2a}}, {{2a,a,a},{2a,a,2a}},{{2a,2a,a},{2a,2a,2a}}}”, в котором наборы координат атомов сгруппированы по парам. Применение “Flatten[atomList,2]” приводит к двухуровневому списку “{{a,a,a},{a,a,2a},{a,2a,a},{a,2a,2a}, {2a,a,a},{2a,a,2a},{2a,2a,a},{2a,2a,2a}}”, который представляет собой набор координат атомов. Это наиболее подходящая для последующего исследования структура списка положений атомов. Итак, в качестве списка положений атомов в простой кубической решетке будем рассматривать результат выполнения команды “n=4; atomList = Flatten[Table[a{i,j,k},{i,n},{j,n},{k,n}],2]” .

Для работы с векторами и матрицами в пакете «Mathematica» реализованы операции скалярного умножения – десятичная точка “.”: “a.b”, где списки a и b могут быть векторами одной длины или матрицами с соответствующими размерностями. Функциональная форма этой операции – функция “Dot[a,b]”. Векторное произведение двух векторов из двумерного или трехмерного пространства задается с помощью команды “ab”, где символ “” вводится с помощью последовательности “ Cross ” (“ ” появляется после нажатия клавиши Esc). Соответствующая функция – “Cross[a,b]” .

Для квадратных матриц может быть найден определитель: “Det[A]”, обратная матрица: “Inverse[A]”, транспонированная матрица любой размерности получается при выполнении команды: “Transpose[A]”. В ряде преобразований надо задавать единичную матрицу, для этой цели в пакете предусмотрена функция “IdentityMatrix[n]”, где n – ее размерность .

Матрица компонент тензора поворота задается с помощью функции “RotationMatrix[]” в двумерном случае, где – угол поворота вокруг оси, ортогональной плоскости, либо с помощью “RotationMatrix[,e]”, где переменная e задает ось поворота (здесь это не обязательно единичный вектор, при вызове указанной функции этот вектор автоматически нормируется). Для получения системы собственных чисел тензора второго ранга используется функция “Eigenvalues[A]”, для получения системы собственных векторов, соответствующих этим числам предусмотрена функция “Eigenvectors[A]”. полная система собственных чисел и векторов находится с помощью функции “Eigensystem[]” .

Для задания тензорного умножения применяется частный вариант функции “Outer[f,list1,list2,…]”, которая возвращает список значений функции f нескольких аргументов на всех возможных комбинациях аргументов, когда первый берется из первого списка list1, второй аргумент – из второго списка list2 и так далее. Например, выполним команду “Outer[f,{a,b,c},{x1,x2}]” --f[a, x1], f[a, x2]}, {f[b, x1], f[b, x2]}, {f[c, x1], f[c, x2]}}” .

Здесь для получения конкретного результата вычислений функция f должна быть задана заранее. В частности, для получения тензорного произведения двух и более векторов в качестве f надо указать имя функции “Times”, реализующей умножение своих аргументов. Функция “Outer[f,list1,list2,L]” имеет дополнительный аргумент L, задающий для многоуровневых списков уровень, элементы которого воспринимаются как единое целое и подставляются целиком в аргументы функции f .

Какой операции тензорного анализа соответствует команда “Outer[D,v,{x,y,z}]”, где v представляет собой векторное поле?

Для нахождения суммы, произведения элементов списка или результата применения какой-либо другой функции ко всему списку сразу используется функция “Apply[f,list]” или в компактной форме “f@@list” .

Эта функция заменяет заголовок “List” списка на имя f. То есть для нахождения суммы элементов списка {a,b,2a,–1,c,3,b} выполняется команда “Plus@@{a,b,2a,–1,c,3,b}” .

При этом происходит замена заголовка списка и вычисляется значение функции вида “Plus[a,b,2a,–1,c,3,b]”, что приводит к результату “2+3a+2b+c”. Для получения произведения элементов надо применить имя функции “Times” .

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

Реализуется эта процедура с помощью команды:

“cMass=Plus@@ atomList/Length[atomList]” .

Для смещения массива атомов в пространстве так, чтобы его центр масс совпал с началом координат, необходимо от радиус-вектора каждого атома отнять вектор, задающий положение центра масс. Команда “atomList – cMass”, очевидно не даст желаемого результата. В данном случае к каждому элементу списка atomList необходимо применить анонимную функцию смещения “ (# – cMass) &”. Применение некоторой функции f к каждому элементу списка list реализуется с помощью команды “Map[f,list]”.

В нашем примере:

“Map[ (# – cMass) &, atomList]” .

Короткий вариант применения функции “Map[]” задается символами “/@”:

“ (# – cMass) & /@ atomList” .

Функция “Map[f,list]” также может использоваться для визуализации сгенерированного кристаллического образца. Действительно, команда “Graphics3D[{Specularity[White,80],ColorData["Atoms","Lu"], Sphere[#, 0.35]} & /@ (atomList /. a 1), Lighting "Neutral"]” позволяет получить геометрическую модель образца с простой кубической решеткой (рис. 1.16, а) .

Часто возникают задачи выбора элементов из списка по некоторому критерию. Например, выясним, как выбрать атомы, принадлежащие одной из граней кубического образца с простой кубической решеткой (их координаты содержатся в переменной atomList). Для выбора элементов списка list по заданному критерию crit применяется функция “Select[list, crit]”. Критерий представляет собой имя функции одной переменной. Выбираются те элементы, для которых вычисленная функция критерия возвращает значение “True”. Рассмотрим грань, ортогональную оси Ox1 для набора атомов atomList, центр масс которого еще не смещен в начало координат. На ребре образца находится n атомов, следовательно, абсцисса атомов нужной грани равна a n, тогда команда “rightFace = Select[atomList, #[[1]] == n a&]” даст требуемый список атомов одной грани (рис. 1.16, б). Получить рисунок с выделенной боковой гранью с использованием функции “Map” .

–  –  –

Рис. 1.16. Образец с простой кубической решеткой:

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

“Complement[atomList, rightFace]” .

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

“force = Plus@@Flatten #1 #2 f[ (#1 #2).(#1 #2) ] &, [Outer[ (#1 #2).(#1 #2) rightFace, Complement[atomList, rightFace],1],1]]” .

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

1. Предопределение: любому действию, которое требуется выполнить, соответствует специальная функция. Пример – классический научный калькулятор. В нашем случае примером предопределения будет команда Factorial[9] --- 362880 .

Можно ввести и “9!” .

2. Программирование на основе списков: для вычисления необходимо предварительно составить список целых чисел (для приведенного примера от 1 до 9) и применить к нему операцию перемножения элементов. Создать список натуральных чисел можно с помощью уже известной функции “Table[]”, но существует и более простая в использовании функция создания числовых списков “Range[]”:

Range[9] --- {1,2,3,4,5,6,7,8,9} Как применять к некоторому списку те или иные функции, мы рассмотрим позже, а пока приведем только один пример. Отметим, что полученный список хранится в ядре пакета в виде конструкции “List[1,2,3,4,5,6,7,8,9]”, а применение к нему функции перемножения элементов (или другого действия) представляется как замена названия функции “List” на новое, в нашем случае – это имя (функция) “Times”. Для замены названия функции на другое применяется функция “Apply[f,expr]”, где f – новое название функции, применяемое к выражению expr. Итак, создадим функцию listFact[n_]: =Apply[Times,Range[n]] listFact[9] --- 362880

3. Рекурсия: идея заключается в том, что сама создаваемая функция выполняет только один оператор, который требуется повторять для решения задачи. То есть код программы сдержит небольшую выполняемую часть и использует результаты своих вычислений, сделанных при других значениях аргументов. Факториал 9 при этом подходе определяется как факториал 8, умноженный на 9:

recurFact[n_]: =n recurFact[n–1]; recurFact[0]=1;

recurFact[9] --- 362880

4. Процедурное программирование: используются операторы цикла, условные операторы и прочее. В качестве примера рассмотрим один из операторов цикла – цикл “Do[expr,iterator]”. Синтаксис этой функции аналогичен синтаксису функции “Table[]”, но операторы цикла, выполняя повторы предписанным образом, не возвращают никакого результата, результаты вычислений можно сохранять в переменных и обращаться к ним после завершения работы цикла. Создадим новую функцию:

loopFact[n_]: =Module[{fact}, fact=1; Do[fact*=i,{i,2,n}]] loopFact[9] --- 362880 Проверьте, чему будет равно вычисленное выражение loopFact[0], и объясните результат. Отметим, что помимо цикла “Do[]” «Mathematica»

дает возможность работать с циклами “For[]” и “While[]”. Здесь на особенностях их использования останавливаться не будем .

Преобразование символьных выражений. Для преобразования символьных выражений используются следующие основные встроенные функции: 1) “Expand[expr]” – раскрытие скобок в выражении, для примера выполните команду “e1=Expand[(x–7)5(2x+3)4(8–x)]”; 2) “Factor[expr]” – разложение на множители, выполните “e2=Factor[e1]”; 3) для тригонометрических выражений используются специальные функции “TrigExpand[expr]” и “TrigFactor[expr]”; 4) для степенных выражений – “PowerExpand[expr]”; 5) общая функция, служащая для упрощения любых выражений, – это функция “Simplify[expr]”. При использовании функции “FullSimplify[expr]” пакет применяет алгоритмы более глубокого упрощения, в том числе специальных функций, что в ряде случаев позволяет получить более простые конечные выражения .

Стандартные математические функции, реализованные в пакете, – это “Power[x,n]” (возведение x в степень n, которая может быть дробной), тригонометрические функции “Sin[x]”, “Cos[x]”, “Tan[x]”, “Cot[x]”, обратные тригонометрические функции “ArcSin[x]”, “ArcCos[x]”, “ArcTan[x]”, “ArcCot[x]”, логарифмические функции “Log[x]” (натуральный логарифм x), “Log[a,x]” (логарифм x по основанию a), экспонента “Exp[x]” .

Найти производную позволяет функция “D[f[var],var]”, где первый аргумент указывает функцию, а второй переменную, по которой берется производная. Пример: “D[f[x,y],x]” – первая частная производная f по x, “D[f[x,y],x,x]” – вторая частная производная f по x, “D[f[x,y],{x,2}]” – также вторая частная производная по x (после запятой в фигурных скобках указывается любой порядок производной), “D[f[x,y],x,{y,3}]” – первая частная производная по x и третья по y .

Интегрирование реализуется с помощью “Integrate[f[var],var]”, во втором аргументе можно также указать и пределы интегрирования: “{var,a,b}” для определенных интегралов. Численное интегрирование реализуется с помощью функции “NIntegrate[f[var],{var,a,b}]”. При этом необходимо задавать числовые значения всех параметров в подынтегральной функции .

Решение уравнений. Основной функцией для точного решения уравнений является функция “Solve[eqn,x]”, где eqn – уравнение, которое может содержать любые параметры, а x – неизвестная величина. Уравнение в пакете задается с помощью логической операции “==”.

Например:

“sol=Solve[x + 6 – a/x == 0, x]” --x 3 9 + a},{x 3 + 9 + a}} ” .

Как видно, решение представляется в виде подстановок. Полученное в примере решение сохраняется во введенной переменной sol. Для графического исследования поведения корней уравнения в зависимости от параметра a построим график (рис.

1.17):

“Plot[x /. sol, {a, -10, 0}, AxesLabel {"a", "x"}, PlotStyle Thick, LabelStyle {Black, FontFamily "Times", 16}]” .

Для решения системы уравнений используется та же функция со списком уравнений и неизвестных: “Solve[{eqn1,eqn2,…},{x1,x2,…}]”. Модификация этой функции “NSolve[]” дает приближенные значения корней уравнения (примените ее к решению записанного выше уравнения). Читателям предлагается решить точно и приближенно уравнение “x4 – 2 x3 – x2 – 2 x + 1 == 0” .

Рис. 1.17. Зависимость значений корней x1, x2 уравнения “x + 6 – a/x == 0” от параметра a Отметим, что функции “Solve[]” и “NSolve[]” предназначены для решения полиномиальных уравнений .

Другая функция, позволяющая получить более полное решение ряда уравнений (в том числе некоторых трансцендентных), – это функция с таким же набором аргументов “Reduce[{eqn1,eqn2,…},{x1,x2,…}]” .

Кроме уравнений эта функция позволяет решать и неравенства, в также системы уравнений и неравенств.

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

–  –  –

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

Пример решения неравенства (выполнить самостоятельно):

–  –  –

Для численного поиска корней трансцендентных уравнений в пакете имеется функция “FindRoot[eqn,{x,x0}]”, которая находит ближайший к x0 корень уравнения eqn. Для нахождения нуля функции вместо уравнения eqn подставляется имя функции f. Эта функция использует метод Ньютона поиска корня. При задании двух начальных значений “{x,x0,x1}” функция будет применять метод секущих. Эта функция может находить и комплексные корни, если задается начальное комплексное значение .

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

Для нахождения минимума функции одного или нескольких аргументов используются функции: 1) “Minimize[f[x1,x2,…],{x1,x2,…}]”, причем дополнительно можно задавать ограничения cons в виде равенств и неравенств “Minimize[{f[x1,x2,…], cons},{x1,x2,…}]”, поле dom (Reals, Integers), над которым ведется поиск решения “Minimize[f[x1,x2,…],{x1,x2,…}, dom]”, 2) “NMinimize[f[x1,x2,…],{x1,x2,…}]”, причем можно также задавать ограничения, выбирать метод численного решения и управлять его параметрами, 3) “FindMinimum[f[x1,x2,…],{x1,x2,…}]”, здесь можно задавать метод численного решения, управлять его параметрами, выбирать начальную точку для поиска минимума, задавать ограничения. Отметим, что первая и вторая функции “Minimize[]” и “NMinimize[]”ведут поиск глобального минимума, а функция “FindMinimum[]” ищет локальный минимум. Метод поиска решения задачи минимизации задается опцией “Methodname”, где name – название метода, принимающее одно из следующих значений: "Automatic", "Gradient", "ConjugateGradient", "InteriorPoint", "QuasiNewton", "Newton", "LinearProgramming", "QuadraticProgramming", "LevenbergMarquardt", названия которых говорят сами за себя. Значение "Automatic" принято по умолчанию и означает автоматический выбор одного из перечисленных методов в соответствии с некоторым предварительным анализом постановки .

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

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

1. Босс В. Интуиция и математика. – М.: КомКнига. 2007. – 192 с .

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

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

3. Воробьев Е.М. Введение в систему «Математика»: учеб. пособие. – М.: Финансы и статистика, 1998. – 262 с .

4. Гольдштейн Р.В., Городцов В.А. Механика сплошных сред. Ч. 1. Основы и классические модели жидкостей. – М.: Наука. Физматлит, 2000. – 256 с .

5. Димитриенко Ю.И. Нелинейная механика сплошной среды. – М.:

Физматлит, 2009. – 624 с .

6. Дьярмати И. Неравновесная термодинамика. Теория поля и вариационные принципы. – М.: Мир, 1974. – 304 с .

7. Жилин П.А. Векторы и тензоры второго ранга в трехмерном пространстве. – СПб.: Нестор. 2001. 276 с .

8. Жилин П.А. Актуальные проблемы механики // Актуальные проблемы механики: материалы междунар. летней школы-конференции. Т. 1. Институт проблем машиноведения РАН. – СПб., 2006. – 306 с .

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

1986. – 232 с .

10. Пригожин И., Кондепуди Д. Современная термодинамика. От тепловых двигателей до диссипативных структур. – М.: Мир, 2002. – 461 с .

11. Седов Л.И. Механика сплошной среды: в 2 т. Т. 1. – СПб.: Лань, 2004. – 528 с .

12. Седов Л.И. Механика сплошной среды: в 2 т. Т. 2. – СПб.: Лань, 2004. – 560 с .

13. Трусделл К. Первоначальный курс рациональной механики сплошных сред. – М.: Мир, 1975. – 592 с .

ГЛАВА 2. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

ПОВЕДЕНИЯ ДЕФОРМИРУЕМЫХ ТВЕРДЫХ ТЕЛ

НА АТОМАРНОМ УРОВНЕ. ДИСКРЕТНЫЙ ПОДХОД

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

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

2.1. ИСТОРИЯ ДИСКРЕТНОГО ПОДХОДА .

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

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

Одним из способов теоретического исследования термомеханических свойств конденсированных сред является дискретный подход, основанный на прямом моделировании движения частиц, вызванного их взаимодействием и приложенными внешними воздействиями. История дискретного подхода в статической постановке начинается в 1827 году с работ Клода Луи Мари Анри Навье (1785–1836) и Огюстена Луи Коши (1789–1857) об определении структуры линейного упругого закона и числа независимых упругих модулей. Навье вывел уравнение равновесия упругих тел, исходя из предположения, что идеальное упругое тело состоит из молекул, между которыми возникают силы взаимодействия. При выводе он принимал, что силы взаимодействия направлены по прямым, соединяющим частицы тела, и пропорциональны расстояниям между ними (первое приближение для потенциала межатомного взаимодействия). Для изотропных тел его уравнения содержали только одну упругую постоянную. Примерно в то же время и Коши исследовал строение упругого закона .

Он изначально ввел два упругих модуля для изотропного материала. В самом общем случае принимал, что каждая из 6 компонент тензора напряжений является однородной линейной функцией от 6 компонент тензора деформаций, то есть для описания упругого поведения твердого тела необходимо не более 36 упругих постоянных. Принимая молекулярную теорию, Коши снизил число упругих постоянных в обобщенном законе Гука с 36 до 15. В изотропном случае он получил, как и Навье, только одну упругую постоянную. Молекулярной теории при анализе упругих свойств материалов придерживался и Симеон Дени Пуассон (1781–1840). В 1829 году он опубликовал работу, в которой на основе молекулярной теории показал, что при простом растяжении стержня осевое удлинение сопровождается поперечным сужением на величину = 1 / 4.

И если E модуль Юнга, то модуль G в опыте на простой сдвиг изотропного материала должен быть равен:

E G= = 0,4 E .

2(1 + 1 / 4) Проведенные впоследствии эксперименты показали результаты, расходящиеся с представлениями одноконстантной теории. Хотя всегда оставалось сомнение, что испытуемый материал в действительности является изотропным. Поэтому представление о возможности полностью оценить упругие свойства изотропного тела одной упругой постоянной на ранних стадиях развития теории упругости пользовалось всеобщим признанием (Навье, Коши, Пуассон, Ламе, Клапейрон…) .

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

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

Работа Грина явилась отправным пунктом долгой дискуссии о количестве упругих модулей и привела к формированию двух школ в теории упругости: «упругость по Коши», «упругость по Грину». Таким образом, первое противопоставление дискретного и непрерывного подхода возникло со спора о количестве упругих постоянных для изотропной линейно упругой среды .

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

Построим на этой площадке в направлении, противоположном n, циРис. 2.1. Цилиндр для определения напряжений Коши в упругом теле линдр высотой l, которая превышает расстояние, на котором еще существенно взаимодействие между парой частиц (рис. 2.1). Вектор напряжений на площадке n dS вызван силами, действующими на частицы с радиусом-вектором b цилиндра со стороны всех частиц с радиусом-вектором a из «положительного»

(см. рис. 2.1) полупространства .

Принимаются следующие предположения:

1. В результате деформирования из начальной ненапряженной конфигурации частицы a и b перемещаются в положения a' и b'. Поле вектора перемещений u(x) будем считать непрерывным однозначным на всем упругом пространстве. Тогда a ' = a + u(a), b ' = b + u(b), a ' b ' = a b + u(a) u(b) = a b + uab. (2.1) Значения вектора перемещений u(a) и u(b) определены для точек упругого континуума, в которые попадают рассматриваемые частицы .

2. Предположим, что сила парного взаимодействия fab между частицами a и b направлена по линии a ' b ', а ее величина зависит от расстояния между частицами и пропорциональна их массам ma, mb :

fab = ma mb f ( a ' b ' ) (a ' b ') / a ' b '. (2.2)

3. Функцию f (r ) будем считать непрерывной вместе со всеми производными любого порядка .

4. Функция f (r ) с ростом расстояния r настолько быстро стремится к нулю, что в рассмотрении остаются значения f (r ) только для бесконечно малых r .

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

–  –  –

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

В силу предположения 6 число частиц N в цилиндре высотой l = (a ' b ') n при заданной плотности среды

–  –  –

В общем случае анизотропии это представление дает 15 независимых упругих модулей, для тетрагональной сингонии материала – 6 упругих модулей, для кубической сингонии – 3 модуля. В изотропном случае вследствие этого соотношения остается только 1 упругий модуль. Это стало причиной возникновения спора о количестве упругих постоянных .

Напомним, что сейчас мы работаем с двумя упругими модулями для изотропного случая .

Рассмотрим далее определяющие соотношения линейно-упругой среды в представлении Коши.

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

= C: .

В силу симметрии тензоров и тензор жесткости C симметричен относительно перестановок внутри первой и последней пар индексов:

Cijkl = C jikl = Cijlk .

Плотность внутренней энергии для линейно-упругой среды

u = : = :C:

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

Cijkl = Cklji .

Тензоры четвертого ранга, удовлетворяющие этим условиям, называются полусимметричными. Можно показать, что в трехмерном случае тензор C имеет 21 независимую ненулевую компоненту, а в двумерном случае – 6 .

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

Рассмотрим закон преобразования компонент тензора C при наложении ортогонального преобразования (например, при повороте материала):

–  –  –

где Q – ортогональный тензор. В случае, когда ортогональное преобразование Q согласовано с симметрией материала, получаем C* = C .

При этом говорят, что преобразование Q принадлежит группе симметрии тензора линейно-упругих свойств: Q LC.

Записанное условие соответствует системе дополнительных ограничений, накладываемых на компоненты C :

–  –  –

в которой все индексы пробегают значения от 1 до n, где n – размерность пространства. В работе [12] установлен изоморфизм 21 мерного пространства упругих модулей для 3-мерной среды и эрмитова пространства той же размерности. Приведенные соотношения позволяют проводить анализ структуры тензора C для различных классов симметрии материала .

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

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

Уравнения движения имеют вид (например, [6]) N N mrk = F(rki )rki + (rki, v ki )rki + (rk ) + (rk, v k ), (2.12) i =1 i =1 где m – масса частицы, rki rk ri – разность радиусов-векторов частиц, rki = rki, v ki v k v i – разность скоростей частиц, члены F(rki ) и (rki, v ki ) описывают консервативную и неконсервативную составляющие силы взаимодействия между частицами, векторзначные функции (rk ) и (rk, v k ) описывают внешнее консервативное и неконсервативное силовое поле соответственно. Приняты обозначения: F (r ) = f (r ) / r, f (r ) = '(r ), f (r ) – скалярная сила взаимодействия частиц, (r ) – потенциал взаимодействия .

Неконсервативное взаимодействие (rki, v ki ) предназначено для описания внутренней диссипации в теле. Внешние силовые поля (rk ) и (rk, v k ) обычно используются для двух целей – для задания внешних массовых силовых воздействий (гравитационного, электромагнитного) и для задания силовых граничных условий. В первом случае эти силы распределены во всем объеме пространства, где проводится расчет, во втором случае они локализованы вблизи некоторых поверхностей, часто являющихся границами расчетной области. Кроме того, неконсервативное воздействие (rk, v k ) часто используют для описания отвода энергии из системы посредством внешней диссипации, простейшим вариантом которой является сила вязкого трения (rk, v k ) = B v k, B 0. Эта сила также используется для внешнего контроля над тепловым движением атомов в системе, в этом случае коэффициент В может быть знакопеременным и зависеть от уровня тепловой энергии всей системы .

В простейшем случае в отсутствие неконсервативного взаимодействия и внешних сил получим уравнение движения вида

–  –  –

Система (2.12) представляет собой задачу Коши. Начальными условиями являются распределение частиц в пространстве rk и их скоростей v k .

Эти распределения могут быть самыми разными в зависимости от целей исследования, но обычно для исключения движения рассматриваемой системы частиц как жесткого целого должны выполняться требования N N mk v k = 0, m (r r0 ) v k = 0. (2.14) k k k =1 k =1 Прежде чем перейти к описанию потенциалов взаимодействия частиц, перечислим кратко основные типы связей в кристаллических телах (рис. 2.2). Отметим, что в металлах формируется обобщенный электронный газ, приводящий к сферической симметрии взаимодействия атомов .

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

Схематичное представление распределения электронов в твердых кристаллических телах:

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

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

Для выражения законов движения частиц используется уравнение Шредингера, которое связывает энергию электронной системы с волновой функцией. Орбиталь – пространство вокруг ядра, в котором наиболее вероятно нахождение электрона. Атомные орбитали отличаются энергией, размером, формой и положением в пространстве относительно ядра. Согласно квантово-механическим расчетам s-орбитали имеют форму шара (рис. 2.3, а), p-орбитали – форму гантели (рис. 2.3, б), d-орбиталь в зависимости от характеризующих ее квантовых чисел может принимать две различные формы (рис. 2.3, в), f-орбиталь – четыре различные формы (рис. 2.3, г). При образовании ковалентных кристаллов главную роль играет гибридизация орбиталей .

а б в

–  –  –

Гибридизация орбиталей – это изменение формы некоторых орбиталей при образовании ковалентной связи для достижения более эффективного перекрывания орбиталей. Выделяют sp-гибридизацию (рис. 2.4, а), в которой участвуют атомные орбитали одного s- и одного p-электронов .

В ряде соединений sp-гибридизации подвергаются атомы углерода, соединяющиеся между собой тройными связями. При этом гибридные орбитали атомов углерода образуют две -связи с соседними атомами, негибридные орбитали атомов углерода образуют две -связи. В sp2-гибридизации (рис. 2.4, б) участвуют атомные орбитали одного s- и двух p-электронов .

Атомы углерода, находящиеся во втором валентном состоянии (sp2-гибридизация), связаны друг с другом двойными химическими связями (решетки типа графена). В sp3-гибридизации участвуют атомные орбитали одного s- и трех p-электронов (рис. 2.4, в), подобные связи для углерода обеспечивают существование решеток типа графита (или терморасширенного графита) и алмаза .

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

–  –  –

Потенциал Леннард-Джонса. Двухчастичный потенциал ЛеннардДжонса и соответствующая сила взаимодействия имеют вид 12 12

–  –  –

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

Потенциал Ми. Потенциал и сила взаимодействия Ми имеют вид

–  –  –

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

Модифицированный потенциал. Часто требуется изменить степень дальнодействия потенциала, сохранив его основные свойства.

Для этого используется следующий модифицированный потенциал:

(r ) = (k (r ) + ), (2.18) где – некоторый произвольный потенциал взаимодействия. При k = 1 исходный и модифицированный потенциалы совпадают, при k 1 модифицированный потенциал «сжат» по сравнению с исходным, при k 1 – он «растянут». Отметим, что () = (), следовательно, сжатие и растяжение происходят относительно точки равновесия r =. Тогда сила взаимодействия

f (r ) = k f (k (r ) + ). (2.19)

Отметим, что модифицированный потенциал, полученный из потенциала Морзе, является потенциалом Морзе с измененным значением = k. Модифицированный потенциал, полученный из потенциала Леннард-Джонса при k = / 6, эквивалентен потенциалу Морзе, определяемому параметром. Под эквивалентностью здесь понимается, что у потенциалов могут совпадать три важнейших размерных параметра: равновесное расстояние, жесткость связи и энергии связи. Таким образом, использование модифицированного потенциала позволяет расширить взаимодействие Леннард-Джонса, сделав его трехпараметрическим по аналогии с потенциалом Морзе .

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

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

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

r (t + ) = 2r (t ) r (t ) + w (t )2, (2.20)

где – шаг интегрирования, w (t ) – ускорение частицы, получаемое подстановкой рассчитанных значений r (t ) в правую часть уравнений (2.12) .

Данная схема не требует вычисления скоростей и удобна, если в уравнениях (2.12) отсутствуют неконсервативные силы .

Метод центральных разностей, видимо, впервые примененный в МД Виньярдом, определяется уравнениями v (t + / 2) = v(t / 2) + w (t ), (2.21) r (t + ) = r (t ) + v (t + / 2) .

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

Более точного расчета позволяют добиться методы типа предикторкорректор, в частности метод Нордзика. Методы предиктор-корректор позволяют получить значительно более высокую точность при малых шагах интегрирования, чем метод Верле и родственные ему методы. Однако для достаточно больших шагов интегрирования, которые часто применяются в методе частиц, метод Верле может оказаться точнее, чем метод Нордзика порядка 3, 4 и 5 .

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

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

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

–  –  –

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

–  –  –

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

Учитывая, что любой потенциал лишь приближенно описывает характер взаимодействия атомов, введем аналог потенциала Леннард-Джонса – периодический потенциал «6–12»:

–  –  –

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

Задания:

1) Построить график зависимости значений периодического потенциала от расстояния между атомами и сравнить с графиком для потенциала Леннард-Джонса при учете нескольких образов области моделирования .

2) Провести расчет динамического поведения нескольких частиц на ячейке периодичности с потенциалом (2.26), исследовать изменение кинетической и потенциальной энергий .

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

1. Какие типы связей могут образовываться в кристаллических телах?

2. Гибридизации какого типа подвержены атомы углерода в структурах графена и графита?

3. Что качественно описывают потенциалы взаимодействия атомов?

4. Получить линеаризацию потенциалов Леннард-Джонса, Ми и Морзе для вычисления упругих модулей методом Коши при условии малых смещений .

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

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

2.2. МЕТОД АТОМАРНОЙ СТАТИКИ. ОПРЕДЕЛЕНИЕ УПРУГИХ

МОДУЛЕЙ МАТЕРИАЛОВ С РАЗЛИЧНОЙ КРИСТАЛЛИЧЕСКОЙ РЕШЕТКОЙ

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

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

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

Если исследуемый материал является линейно-упругим по Коши с несимметричным тензором напряжений, удовлетворяющим обобщенному закону Гука с несимметричными мерами: = C : u, T, Cijkl C jikl, Cijkl Cijlk, Cijkl = Cklij, где – оператор градиента из отсчетной конфигурации, то полученные коэффициенты при первых степенях параметров деформации – компонент тензора дисторсии u T, – будут компонентами тензора линейно-упругих свойств С. Если материал является (нелинейно) упругим по Грину, то из разложения упругого потенциала в степенной ряд по дисторсии получается упругий закон также в виде ряда, первый член которого совпадает с правой частью записанного несимметричного обобщенного закона Гука .

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

Материал с кубической решеткой имеет три взаимно ортогональные оси симметрии четвертого порядка [12], совпадающие с кристаллографическими осями (рис. 2.6, а). С учетом всех симметрийных и анизотропных свойств тензор линейно-упругих свойств C имеет четыре независимые компоненты, в качестве которых выбираются компоненты с наименьшими значениями индексов в указанных кристаллографических осях: C1111, C1122, C1212, C1221. Остальные компоненты либо равны нулю, либо выражаются через них: C1133 = C3311 = C2233 = C3322 = C2211 = C1122, C1313 = C2121 = C2323 = = C3131 = C3232 = C1212, C1331 = C3113 = C2332 = C3223 = C2112 = C1221, C2222 = C3333 = C1111 .

Ограничения на значения упругих модулей, следующие из условия положительной определенности тензора линейно-упругих свойств C ( e : C : e 0, e 0 ; e : C : e = 0 e = 0 ), получаются из неравенства

–  –  –

Из анализа этого выражения при различных eij или с помощью критерия Сильвестра, примененного к матрице 9 9 с компонентами C{i, j}{k,l }, где каждая пара индексов {i, j}, {k, l} пробегает значения из набора пар вида {{1,1},{1,2},{1,3},{2,1},{2,2},{2,3},{3,1},{3,2},{3,3}}, выводятся следующие ограничения на компоненты:

C1111 0, 1 2C1111 C1122 C1111, C1212 0, C1212 C1221 C1212. (2.28) Четвертое (двойное) неравенство (*) является нестрогим, хотя и обеспечивает строгую положительную определенность тензора C.

Это связано с тем, что последние два слагаемых в неравенстве (2.28) в случае равенств модулей вида C1221 = ± C1212 после приведения подобных дают строго положительную при ненулевых компонентах несимметричного тензора e сумму полных квадратов сумм или разностей:

–  –  –

Простой сдвиг монокристалла. Пусть n – единичная нормаль к плоскости сдвига, m – единичный вектор, задающий направление сдвига. Закон движения в лагранжевой форме имеет вид: r = R + ( n R ) m, где R – радиус-вектор центров атомов в отсчетной конфигурации, а r – радиус-вектор атомов в текущей конфигурации. Тогда деформационный градиент F = ( r ), описывающий сдвиг на величину, имеет вид F = E + mn. РасT смотрим верхнюю и правую боковую грани. Площадь верхней грани при описанном сдвиге не меняется и равна N 2 a 2, так как на каждый атом на грани приходится ее часть площадью a 2, где a – параметр решетки, N – количество атомов вдоль ребра куба (на рис. 2.6 N = 4 ). Площадь деформируемой грани меняется согласно соотношению ds = J (N U 2 N)1/2 dS, где N – единичный вектор нормали к недеформированной грани, U 2 = F 1 F T, F 1 = E mn, J = det F, dS – площадь грани до деформирования, ds – площадь грани после деформирования [9]. Для правой боковой грани получим площадь после сдвига, равную N 2 a 2 1 + 2 .

Выражение для суммарной силы, действующей на любую из граней, для произвольной силы f ( r ) r / r взаимодействия пары атомов, соединенных вектором r, могут быть найдены аналитически с помощью любого пакета символьных вычислений, например пакета Wolfram Research «Mathematica». Для произвольных значений межатомного расстояния a может оказаться, что в состоянии до приложения сдвига суммарные силы на гранях исходного куба отличны от ноля, то есть в качестве начальной конфигурации берется не естественное состояние (ненапряженное и недеформированное), а некоторое равноосно деформированное состояние. Для получения естественного состояния необходимо решать задачу о выборе естественного периода решетки. Эта проблема учитывается далее при нахождении упругих модулей .

Нормаль n к деформированной грани определяется выражением 2 1/2 T

n = (N U N) F N [9]. Например, для правой боковой грани с нормалью N R = {1,0,0} в недеформированной конфигурации после рассматриваемого простого сдвига, получим нормаль в текущей конфигурации:

–  –  –

где индекс “R” (Right) относится к правой боковой грани, индекс “Up” (Upper) к верхней, индекс “Fr” (Front) к передней грани. В частности, полученные выражения дают возможность сравнивать недиагональные компоненты тензора напряжений 12 ( ) и 21 ( ). Для определения упругих модулей функции 12 ( ) и 21 ( ) раскладываются в степенные ряды в окрестности ноля ( = 0 ), коэффициенты при первых степенях в которых будем называть модулями сдвига G12 и G21 : 12 ( ) = 12 (0) + G12 + 2! 12 ( ) =0 2 +.. .

Значения компонент тензора напряжений Коши в начальной конфигурации 12 (0) и 21 (0) не оказывают влияния на значения исследуемых модулей .

Для линейно-упругого материала найденные таким образом касательные (упругие) модули будут совпадать с секущими G12 ( 12 ( ) 12 (0) ) / и G21 ( 21 ( ) 21 (0) ) / модулями. Далее рассматриваются только касательные модули .

Рис. 2.6. Монокристаллический образец с ГЦК-решеткой:

а – до сдвига; б – после простого сдвига в плоскости X1OX2 в направлении X1 (атомы в вершинах ячеек изображены серым цветом, на гранях – коричневым)

Закон Гука (2.27) для рассматриваемого простого сдвига принимает вид:

–  –  –

Выберем в качестве потенциала взаимодействия атомов (двухчастичный) потенциал Леннард-Джонса. Используя конкретный вид силы, можно в явном виде получить набор коэффициентов разложения в ряды Тейлора компонент 12 ( ) и 21 ( ) тензора напряжений Коши в окрестности ноля ( = 0 ). Заметим, что при аргументе r = сила взаимодействия атомов обращается в ноль, то есть параметр есть равновесное расстояние для изолированной пары атомов, но a, так как параметр решетки a определяется из условия равновесия большого количества атомов. Действительно, даже для одномерной цепочки из 3 атомов крайние атомы приблизятся к центральному атому, и равновесное расстояние a в такой системе будет отличаться от величины .

Для ГЦК-решетки, начиная с количества атомов N = 2, были получены коэффициенты разложения в ряд Тейлора при первых девяти степенях (табл. 2.1–2.4). Для коэффициентов степенного ряда компоненты 12

–  –  –

Оказалось, что коэффициенты при первых степенях при разложении компонент 12 ( ) и 21 ( ) совпадают. Коэффициенты при четных степенях на 15–18 порядков меньше коэффициентов при первых степенях, и ими можно пренебречь. Коэффициенты при третьих, пятых и более высоких нечетных степенях для разложений в ряды Тейлора 12 ( ) и 21 ( ) оказались различными. Но при малых значениях сдвига (до величины = 0,1 ) члены разложения с третьей и выше степенями дают пренебрежимо малый вклад в выражения для компонент 12 ( ) и 21 ( ). Дадим оценку вклада этих членов, принимая для простоты, что a (см. табл. 2.4) .

–  –  –

Для сдвига на величину = 0,1 и числа атомов N = 12 на ребре исходного куба от слагаемого третьего порядка получим добавку 932,389 / a 3 к значению компоненты 12 и добавку 932,153 / a 3 к значению компоненты 21. Относительное различие в полученных значениях составит примерно 1,72 103 %. Для величины сдвига = 0,05 различие будет значительно меньшим. Таким образом, отклонение тензора напряжений от симметричного вида для монокристалла с ГЦК-решеткой в условиях малых упругих деформаций мало, определяется членами ряда с более высокими степенями параметра сдвига и зависит от величины .

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

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

Например, графическое представление зависимостей от числа атомов N первого и второго числовых коэффициентов C1G12 ( N ) и C2 12 ( N ), стоящих G

–  –  –

не сводятся к таким, не имеющим горизонтальной асимптоты функциям, как логарифмическая и степенная функция с меньшим единицы положительным показателем. Окончательно зависимости искались в классе функций f ( x) = a ( x x0 ) k + b. Методом наименьших квадратов получены аппроксимирующие выражения C1G12 ( N ) = 1343,13( N + 0,47) 1,01 811,86 N коэффициент C1G12 стремится к значению –811,86) и (при C2 12 ( N ) = 23608,90( N + 0,41) 1,008 + 15367,82 (при N этот коэффиG циент стремится к значению 15367,82). Показатели степени близки к –1, то есть кривые на рис. 2.7 являются ветвями гипербол с центрами в точке с абсциссой, примерно равной –0,45. Итак, упругий модуль G12 представляется в виде

–  –  –

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

В упрощенном случае, когда принимается приближенное равенство периода решетки и параметра потенциала a, графики полученных зависимостей недиагональных компонент тензора напряжений 12 ( ) и 21 ( ) в виде рядов с предельными ( N ) значениями числовых множителей приведены на рис. 2.8, а. Относительное несовпадение этих компонент в зависимости от величины сдвига показано на рис. 2.8, б. Как видно из полученных результатов, относительное различие между недиагональными компонентами тензора напряжений 12 ( ) и 21 ( ) в опыте на простой сдвиг ГЦК-монокристалла при упругом деформировании составляет менее одной сотой процента и им можно пренебречь. Поэтому можно говорить только об одном модуле сдвига ГЦК-материалов G = G12 = G21. Другие анизотропные модули C1111 и C1122 определяются в опыте на чистое растяжение-сжатие .

–  –  –

Сами зависимости компонент 12 ( ) и 21 ( ) от величины сдвига для ГЦК-решетки являются нелинейными, а полученные разложения в степенные ряды справедливы только в окрестности ноля. Точное выражение для этих компонент, применимое для произвольных значений параметра сдвига, хотя и выведено, но не может быть представлено в полном виде. Использование пакета символьных вычислений позволяет представить зависимость для него в графическом виде (например, на рис. 2.9 приведен график для случая N = 12 ) .

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

–  –  –

в пределе при N с использованием потенциала межатомного взаимодействия Леннард-Джонса все упругие модули будут выражены через два параметра и .

Приведенная на рис. 2.10 зависимость описывается соотношением

–  –  –

То что период ГЦК-решетки a оказался больше, чем параметр потенциала равновесного расстояния для изолированной пары атомов, объясняется строением ГЦК-кристалла, на гранях элементарных ячеек которого расположены дополнительные атомы, «раздвигающие» атомы, лежащие в вершинах ячеек. Для простой кубической решетки период a будет меньше, чем параметр – в аналогичных расчетах получено значение a = 0,9561, для ОЦК-решетки получена связь a = 1,1251 .

Исследуем вопрос об однородности поля напряжений в изучаемом объеме атомов. Рассмотрим начальную конфигурацию монокристаллического куба для случаев с различным числом атомов N на ребре. Период решетки a в зависимости от N берется различным согласно соотношению (2.36), чтобы обеспечить нулевые напряжения на поверхности куба. Для полученной (отсчетной) конфигурации строится набор сечений по атомным плоскостям, параллельным граням куба. Часть атомов, отделенных от тела проведенным сечением, отбрасывается, а ее влияние заменяется эквивалентными ей компенсирующими силами. В каждом таком сечении вычисляется вектор напряжений. Для этого находится сумма всех сил, действовавших на атомы сечения со стороны всех отброшенных атомов тела, и делится на площадь сечения. Показано, что получаемый вектор напряжений в рассматриваемой (отсчетной) конфигурации всегда направлен строго по нормали к сечению. Получена зависимость единственной ненулевой компоненты вектора напряжений от положения сечения в объеме куба. Для примера приводятся результаты для правой боковой грани куба и вычисляется зависимость 11 от k, k = 1, N. Оказалось, что при постоянном периоде решетки однородность (отличие от ноля) напряжений нарушается только в тонком поверхностном слое, толщина которого составляет порядка 2–3 межатомных расстояний (рис. 2.11, а), и в нем возникают сжимающие напряжения. Это говорит о том, что в этом поверхностном слое атомные плоскости стремятся сблизиться и образовать пленку с большей плотностью, чем у материала в объеме .

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

–  –  –

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

max

–  –  –

При реализации простого сдвига распределение значений нормальных компонент ii тензора напряжений Коши не меняется. На рис. 2.12 приведена зависимость значений касательных напряжений 12, определяемых на внешней грани тела в опыте на простой сдвиг на величину = 105, от числа атомных слоев k, взаимодействие с которыми учитывается при расчете 12 (k откладывается от поверхности в объем тела). Таким образом, получаемое напряженное состояние можно считать однородным .

–  –  –

Приведенные соотношения позволяют найти зависимость упругих модулей кристалла от числа атомов N на ребре исследуемого объема или при известном периоде решетки a от размеров образца L = Na. Графическое представление зависимости модуля сдвига и периода решетки от числа атомов N для меди приведено на рис. 2.13 .

–  –  –

G10 nm = 4,28 1010 (Па). Оказалось, что с уменьшением размеров моноCu кристалла его модуль сдвига также уменьшается, и существует такой “переходный” размер образца (рис. 2.13, а), что при дальнейшем сокращении размеров тела падение модуля происходит очень быстро .

С уменьшением размеров тела равновесное межатомное расстояние растет, следовательно, уменьшается плотность частиц с ГЦК-решеткой ( ( N ) = 4ma / ( a ( N ) ), где ma – масса атома) .

Чистое растяжение-сжатие монокристалла. Далее исследуется поведение материала в опыте на чистое растяжение-сжатие исходного объема с ГЦК-решеткой (рис. 2.14) при сохранении размеров в двух других направлениях (ортогональных оси растяжения-сжатия). Пусть m – единичный вектор, задающий направление оси растяжения-сжатия (для деформации, показанной на рис. 2.14, m = {1,0,0} ). Деформационный градиент, описывающий соответствующую деформацию с кратностью удлинения, имеет вид F = E + ( 1)mm .

Рассмотрим верхнюю, переднюю и правую боковую грани куба, задаваемые соответственно нормалями N Up = {0,1,0}, N Fr = {0,0,1}, N R = {1,0,0} в отсчетной конфигурации. Площадь правой боковой грани при описанном деформировании не меняется и равна N 2 a 2. Площади деформируемых (верхней и передней) граней равны между собой и меняются согласно соотношению ds = J (N U 2 N)1/2 dS [9], где N – единичный вектор нормали к грани в отсчетной конфигурации; U 2 = F 1 F T, F 1 = E + ( 1 1)mm ;

J = det F = ; dS – площадь грани до деформирования, ds – площадь грани после деформирования. Для верхней и передней граней получим одинаковые значения их площади dsup = N 2 a 2.

Нормали к граням при растяn Up = N Up = {0,1,0}, n Fr = N Fr = {0,0,1}, жении-сжатии не меняются:

n R = N R = {1,0,0} .

Для нахождения матрицы компонент тензора напряжений Коши ij в приведенном на рис. 2.14 базисе кристаллографической системы координат используется соотношение Коши n = t n. Записывая это соотношение для верхней, передней и правой боковой граней, в рассматриваемом случае получим выражения для компонент этого тензора: 11 () = ( t R )1, 12 ( ) = ( t R )2, 13 = ( t R )3, 21 ( ) = ( t Up ), 22 ( ) = ( t Up ), <

–  –  –

ругих модулей функции 11 ( ), 22 () и 33 ( ) раскладываются в степенные ряды по параметру в окрестности точки = 1. Коэффициент при первой степени в разложении 11 ( ) будем обозначать как модуль растяжения-сжатия E11, а коэффициенты при первых степенях в 22 () и 33 ( ) разложении – как модули растяжения-сжатия E22 и E33. При этом значения компонент тензора напряжений в начальной конфигурации не будут влиять на значения введенных модулей. Для линейно-упругого материала найденные таким образом касательные (упругие) модули будут совпадать с секущими модулями. Из компонентной формы записи обобщенного несимметричного закона Гука (2.29) для рассматриваемого опыта на чистое растяжение-сжатие получим

–  –  –

Для ГЦК-решетки при использовании потенциала Леннард-Джонса в опыте на чистое растяжение-сжатие были получены коэффициенты разложения в ряд Тейлора при нескольких степенях (часть результатов приведена в табл. 2.5–2.9). Для коэффициентов степенного ряда (упругих модулей) рассматриваемых компонент тензора напряжений вводятся обозначения вида

–  –  –

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

–  –  –

Все найденные зависимости имеют практически равные показатели степени k 1, то есть являются ветвями гипербол, центры которых имеют почти одинаковое значение абсциссы x0 0,4 .

Для выведенного ранее соотношения a = 1,3918 (2.36) получена зависимость “макроскопических” компонент 11, 22 от кратности удлинения, представленная графически на рис. 2.15. Для произвольной величины растяжения (рис. 2.15, а, правая ветвь) наблюдается падение напряжения, связанное с разрушением образца вследствие такого расхождения атомных слоев, что сила взаимодействия между ними стремится к нулю. При сжатии, напротив, сила отталкивания слоев неограниченно возрастает по модулю из-за невозможности “бесконечно близко” сдвинуть атомы .

Макроскопические модули E11 и E22 = E33 при найденной для ГЦКматериалов связи a = 1,3918 выражаются как:

–  –  –

вольно точно попадает в интервал приведенных экспериментальных значений макроскопического модуля Юнга меди (отличие для литой меди составляет порядка 3.5 %). Найденные значения упругих модулей (2.44), в частности, для меди (2.45), удовлетворяют условию (2.28) положительной определенности тензора C .

–  –  –

Коэффициент Пуассона получился одинаковым для всех ГЦК-материалов и согласно (2.30) и (2.44) = 0,343. (2.46) Найденное значение коэффициента Пуассона для металлических монокристаллов с ГЦК-решеткой (2.46) также соответствует условию положительной определенности тензора линейно-упругих свойств C в виде (2.31). Для меди экспериментально определенное значение коэффициента Пуассона Cu = 0,35, для алюминия – Al = 0,34. Результаты применения предложенной в работе методики теоретического расчета величины этого коэффициента дают совпадение с его экспериментальными значениями с точностью 2 % .

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

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

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

Другие типы решеток. В качестве примеров кристаллических решеток материалов, широко используемых в современной технике и технологиях, приведем образец с объемно-центрированной (ОЦК) решеткой, характерной для многих металлов (рис. 2.16, а), образец с решеткой графита а б в

Рис. 2.16. Примеры кристаллических решеток:

а – ОЦК-решетка; б – графит; в – графен (рис. 2.17, б) и лист графена (рис. 2.17, в), который является отдельным слоем графита. Рисунки построены в пакете «Mathematica» по аналогии с рассмотренным ранее алгоритмом. Лист графена можно представить состоящим из двух плоских «треугольных» решеток, вложенных одна в другую. Для наглядности они изображены шариками разных цветов, хотя химически представляют собой одинаковые атомы углерода с ковалентными связями .

–  –  –

ГПУ-решетка. Другим распространенным типом кристаллической решетки металлов наряду с ГЦК- и ОЦК-решетками является гексагональная плотноупакованная (ГПУ) решетка. Такое строение имеют титан, цинк, цирконий, бериллий, магний, кобальт и ряд других металлов. Плотноупакованными называются решетки, в которых при заданном минимальном расстоянии между центрами атомов достигается максимальное число атомов в единице объема. Такие решетки изображают с помощью набора шаров, уложенных слоями. В металлических кристаллах связи сферически симметричны [3]. Из этого следует, что чем выше для металлов плотность упаковки атомов (каждый атом окружен максимальным числом соседей), тем меньше их потенциальная энергия. Большинство металлов имеют наиболее плотно упакованные ГЦК- или ГПУ-решетки .

Используя представление решеток в виде набора шаров, можно представить, что ГПУ-решетка обладает двухслойной периодичностью. Последовательность ее слоев записывают в виде …ABABAB…, причем слои, условно обозначаемые как A и B, соответствуют наиболее плотному периодическому расположению шаров на плоскости (их центры лежат в узлах плоской периодической сетки из правильных треугольников). Структуры этих слоев не отличаются, но один слой сдвинут относительно другого так, что атомы слоя A расположены напротив центров пустот слоя B (см. рис. 2.17, б, в). ГЦК-решетка является трехслойной и представляется последовательностью слоев …ABCABC…, в которой структуры слоев A, B, C также не отличаются. Отличие расположения атомов ГЦК-решетки от описанного выше расположения атомов состоит в том, что если взять слой A в качестве основного слоя, то слои с такой же плотнейшей плоской упаковкой B и C будут различаться типом соответствующих им лунок на слоях типа A, в которых располагаются атомы слоев B и C. Представленный на рис. 2.17 выбор формы образца, обоснован в работе [10] .

Плотность упаковки определяется как отношение объема, занимаемого шарами в элементарной ячейке, к ее полному объему, для ГПУ- и ГЦКрешеток она одинакова. Но ГЦК-решетка является простой – атомы всех слоев находятся в одинаковом положении по отношению к соседям, а ГПУ-решетка является сложной, так как по отношению к атомам из разных слоев (A или B) окружающие их атомы в смежных слоях расположены по-разному относительно кристаллографического базиса. Атомы из слоев A и B часто считают атомами разного типа, хотя их физико-химические свойства не различимы и отличие состоит только в геометрии их окружения. Любая сложная решетка может быть представлена как несколько вложенных одна в другую простых подрешеток. Для ГПУ-решетки выделяют две простые подрешетки, они являются объединением слоев атомов типа A или B. Для материалов с ГЦК-решеткой расстояния между ближайшими атомами всегда одинаковы вне зависимости от их принадлежности одному или разным слоям, что является следствием кубической симметрии. Каждый атом ГПУ-решетки, не лежащий на поверхности образца, окружен 6 (ближайшими) атомами из слоя, которому он сам принадлежит .

Также он окружен 6 (ближайшими) атомами из двух соседних слоев (3 атома снизу и 3 атома сверху). Расстояния от выбранного атома до ближайших соседей из его собственного слоя и до ближайших соседей из других слоев различаются. Поэтому ГПУ-решетка в отличие от ГЦК- или ОЦК-решетки характеризуется не одним, а двумя межатомными расстояниями. Переменную a будем использовать для обозначения равновесного расстояния между атомами одного слоя, а переменную b для обозначения равновесного расстояния между одноименными слоями (удвоенное расстояние между слоями A и B). Для проведения численных экспериментов над ГПУ-кристаллом по исследованию его физико-механических свойств необходимо определить равновесные межатомные расстояния .

Дискретные подходы к исследованию свойств монокристаллов с ГПУрешеткой применяются многими авторами как в динамической, так и в статической постановках (например, [1, 2, 4, 5, 7, 8]) с использованием различных, в ряде работ довольно сложных, потенциалов. Рассмотрим вопрос о выборе формы образца для проведения вычислительных экспериментов с ГПУ-решеткой. За основу возьмем подход атомарной статики и простейшие степенные потенциалы: потенциал Леннард-Джонса и обобщающий его потенциал Ми .

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

Рис. 2.18. Пример двух типов гексагонального образца с ГПУ-решеткой; первый тип: а – вид сверху; б – вид спереди;

второй тип: в – вид сверху; г – вид спереди В качестве боковых граней выбираются наборы атомов из слоя обоих типов – A и B, то есть совокупность двух ближайших к геометрической боковой грани слоев, параллельных оси Ox3. Тогда второй вариант призмы с шестиугольным основанием имеет боковые грани с различным расположением атомов. Например, это видно по «передней» и «задней»

граням на рис. 2.18, г), ортогональным оси Ox2. На этом рисунке передняя и задняя грани состоят из 16 атомов типа A (серый цвет) каждая, но передняя содержит 9 атомов типа B, а задняя грань – 12 атомов типа B (коричневый цвет). Расположение атомов типа B относительно атомов типа A на этих гранях также различается. Это приводит к тому, что при нахождении равновесной конфигурации образца невозможно удовлетворить требованию равенства нулю поверхностных сил на всех гранях при сохранении формы образца. Этому требованию можно удовлетворить при искажении шестиугольника, лежащего в основании призмы, то есть, по сути, введением дополнительного параметра ГПУ-решетки. Такой путь не представляется конструктивным .

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

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

Для исследования механических свойств ГПУ-монокристалла также будем применять дискретный подход атомарной статики [4–5, 10] при использовании потенциала Леннард-Джонса. Выбор потенциала обоснован тем, что при исследовании упругих свойств не рассматриваются процессы, идущие при сверхнизких температурах или при скоростном деформировании, когда важную роль во взаимодействии атомов могут играть квантовые эффекты. Использование статической постановки, в которой положения атомов строго определены, позволяет получить аналитическое решение задачи определения равновесного межатомного расстояния для различных кристаллических решеток и образцов различного размера. При этом точные значения межатомного расстояния получаются для объемов материала с небольшим числом атомов N на ребре основания образца (от 3 до 20). Далее для перехода на макроуровень по этим точным решениям делается предельный переход при стремлении числа атомов N на ребре образца к бесконечности .

При определении параметров равновесной ГПУ-решетки для образца некоторого заданного размера требуется решать систему двух уравнений относительно двух параметров a и b. Для вывода уравнений определяется суммарная сила, действующая на атомы нижней грани, а также на атомы одной из боковых граней призмы, например, передней грани 1 (рис. 2.19). Полученные силы скалярно умножаются на вектор нормали к соответствующей грани, и результат приравнивается нулю. Заметим, что модули сил, действующих на атомы трех различных боковых граней, или произведения этих сил на соответствующие нормали, оказываются одинаковыми для всех боковых граней (см. рис. 2.19) .

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

–  –  –

r = r, где r – вектор, соединяющий центры двух атомов) и учета взаимодействия атомов любой выбранной грани со всеми остальными атомами образца, оказывается полиномиальной очень высокой степени. Аналитически эта система не может быть решена. Для численного решения была использована итерационная процедура. На первом шаге принималось, что второй параметр решетки b равен b = 1,5a и из условия равенства нулю результирующей силы на боковых гранях находилось решение для a, имеющее вид a = k0a, где – параметр потенциала Леннарда-Джонса, равный равновесному расстоянию между двумя атомами, а k0a – коэффициент, найденный при численном решении уравнения для a. Далее, при найденном значении a и произвольном b определялась результирующая сила на верхней грани. Из условия, что она также должна быть равна нулю, находилось нулевое приближение к значению b в виде b = k0. Итерационная b процедура поиска равновесных расстояний стартует с полученных нулевых приближений к искомым значениям параметров ГПУ-решетки. Найденное значение b подставляется в уравнение, обеспечивающее нулевую силу на боковой грани, и определяется новое приближение параметра a, которое подставляется в выражение для силы на верхней грани для получения следующего приближения параметра b. Такая процедура сходится достаточно быстро – через 8–10 итераций значения коэффициентов kia,b не изменяются в первых 6 значащих цифрах. Для проверки сходимости проводилось 100 итераций, подтверждающих ее монотонность. Было показано, что все значения параметров a и b ложатся на гладкие кривые. Для них были построены аппроксимирующие зависимости, которые искались в классе функций вида y ( x) = c1 ( x + x0 ) k + c2.

Получено, что показатель степени равен k = 1 с точностью до 4-го знака после запятой:

a ( N ) / = 0,024( N 1,246) 1 + 0,979, (2.47) b( N ) / = 0,060( N + 2,552) 1 + 1,599, где N – число атомов на ребре основания образца. Заметим, что с ростом N межатомное расстояние убывает, то есть образцы меньшего размера обладают меньшей плотностью и состоят с точки зрения механики из другого материала, хотя химически они состоят из тех же атомов .

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

Это позволяет идентифицировать параметр потенциала Леннард-Джонса: при N получаются теоретически рассчитанные значения равновесных параметров ГПУ-решетки:

a = 0,979, b = 1,599, b* / a* = 1,634 .

Для титана известны параметры a Ti = 2,951 (А), bTi = 4,697 (А), то есть bTi / a Ti = 1,592, что неплохо соответствует расчетам (отклонение составляет 2,6 %). Для циркония a Zr = 3,231 (А), b Zr = 5,148 (А), то есть b Zr / a Zr = 1,593. Для цинка a Zn = 2,665 (А), b Zn = 4,947 (А), то есть b Zn / a Zn = 1,856 (значительное отклонение). Для бериллия отношение составляет 1,567, для магния – 1,624, для гафния – 1,580, для рутения – 1,582, для кобальта – 1,632 (лучшее соответствие), для кадмия – 1,886 (большое отклонение). Для проведения расчетов по теоретическому определению физико-механических свойств материалов с ГПУ-решеткой, например их упругих модулей, коэффициента теплового расширения и зависимости упругих свойств от температуры, для каждого конкретного материала необходимо достичь известного соотношения параметров решетки. Потенциал Леннард-Джонса, как видно из полученных результатов, не всегда подходит для такого анализа .

Определение равновесных параметров ГПУ-решетки с помощью потенциала Ми. Потенциал Леннард-Джонса не позволяет описывать разнообразие металлов с ГПУ-решеткой. Его важным преимуществом по сравнению с часто используемым трехпараметрическим потенциалом Морзе (например, в [6–7]), является степенная форма, позволяющая получать аналитические решения, что затруднительно сделать при выборе потенциала Морзе.

В качестве альтернативы рассматривается потенциал более общего вида, чем потенциал Леннард-Джонса, но также степенного вида – потенциал Ми [1], содержащий дополнительные безразмерные параметры m и n:

–  –  –

Потенциал Ми содержит четыре параметра и дает большую свободу по сравнению с потенциалом Леннарда-Джонса при теоретическом исследовании свойств материалов. В рассматриваемой задаче использование потенциала Ми позволяет получать различные отношения параметров решетки для разных металлов с ГПУ-решеткой .

Исследовалась зависимость значений параметров ГПУ-решетки от значений параметров потенциала Ми. Из условия, что параметр m n, последовательно брались значения параметра n из интервала 2 n 51 и перебирались значения параметра m от n + 1 до 50. В результате оказалось, что отношение равновесных параметров ГПУ-решетки для макроскопического уровня изменяется в пределах [1.621; 1.637], что не позволяет описать все металлы с ГПУ-решеткой, но дает возможность выбрать параметры потенциала Ми, подходящие для ряда таких материалов .

Результаты соответствия значений параметров потенциала Ми и отношения приведены в табл. 2.10. Приведенные параметры позволяют получать равновесные межатомные расстояния ГПУ-монокристалла различных металлов и проводить дальнейшие исследования их механических свойств .

–  –  –

Таким образом, с использованием потенциала Ми появляется возможность с большей точностью описывать геометрические характеристики кристаллов с ГПУ-решеткой – равновесные межатомные расстояния и их отношение. При этом удалось определить три параметра потенциала Ми: два показателя степени m и n, а также параметр, задающий равновесное расстояние для изолированной пары атомов выбранного материала. Параметр может быть определен на основании описанного ранее подхода при исследовании упругих свойств монокристаллов ГЦК-решеткой .

Определение равновесных параметров ГПУ-решетки с помощью энергетического подхода. Использованный выше подход к нахождению равновесных значений параметров решетки не является единственным .

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

Для рассматриваемой формы образца (рис. 2.17) с помощью потенциала Леннард-Джонса была найдена энергия U взаимодействия всех атомов решетки. Были получены значения энергии для решетки с вычисленными ранее равновесными параметрами ГПУ-решетки a и b (2.47), а также для образцов с произвольными параметрами решетки, которые подлежат определению из условия минимизации энергии. Для определения значений параметров ГПУ-решетки вторым способом находился минимум энергии взаимодействия атомов U по параметрам a и b .

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

–  –  –

Исследование асимптотического поведения показало, что предельное значение потенциальной энергии взаимодействия атомов ГПУ-решетки при N (при значениях a = 0,979, b = 1,599 ) составляет

–  –  –

Также может быть определена потенциальная энергия U взаимодействия атомов при произвольных параметрах решетки a и b для различного числа атомов N на ребре основания образца с использованием потенциала Леннард-Джонса. Затем из условия минимума потенциальной энергии были определены параметры решетки для каждого числа N. Функция, аппроксимирующая минимальное значение энергии U min ( N ) при каждом N, качественно имеет такой же вид, как и на рис. 2.20. Предельные значения параметров решетки также определяются при N. Заметим, что «макроскопическое» значение энергии решетки ГПУ-металла, полученное в пределе при N, составляет U min / = 8151,16. (2.50) Отличие предельных значений (2.49) и (2.50) составляет 0,06 %. Таким образом, силовой подход к определению равновесных межатомных расстояний в ГПУ-решетке обеспечивает получение равновесного состояния образца при сохранении однородного распределения атомов в решетке, находящегося близко к состоянию с минимальной потенциальной энергией .

Зависимость параметров a и b, обеспечивающих минимум потенциальной энергии взаимодействия, от числа атомов N принимает вид a ( N ) / = 0,067( N 1,196) 1 + 0,972, b ( N ) / = 0,101( N + 2,409)1 + 1,585, предельные значения параметров решетки а = 0,972, b = 1,585 .

Значения параметра a, полученные с помощью силового подхода и на основе минимизации потенциальной энергии взаимодействия атомов образца, отличаются на 0,72 %. Для параметра b отличие составляет 0,93 %. Эти различия пренебрежимо малы, и в расчетах по определению упругих модулей образцов с ГПУ-решеткой удобнее использовать первый метод, гарантирующий точное выполнение заданных статических граничных условий .

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

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

При растяжении образца вдоль осей Ox1 и Ox2 с кратностями удлинения 1 и 2 слои A и B ГПУ-решетки деформируются одновременно, что приводит к нарушению расположения атомов типа B в лунках между атомами типа A. Этот факт продемонстрирован на рис. 2.21, а. Рассмотрим отдельный треугольник, образованный атомами типа A, в лунке между которыми (в проекции на плоскость слоя A ее центр C совпадает с центром описанной окружности) расположен атом типа B. До деформирования выбранный треугольник является равносторонним (рис. 2.21, б) .

После растяжения вдоль осей Ox1 и Ox2 – равнобедренным. Положение C' атома типа B после такого деформирования (полупрозрачный шар на рис. 2.21, в) не соответствует центру C'' описанной окружности треугольника, образованного атомами типа A. Новый центр C'' может быть найден из решения геометрической задачи (шар коричневого цвета на рис. 2.21, в). Тогда при растяжении вдоль осей Ox1 и Ox2 атомам типа B надо давать дополнительное смещение относительно слоя типа A, зависящее от кратностей удлинения 1 и 2 и приводящее к уменьшению потенциальной энергии системы атомов (рис.

2.22):

CC'' = a. (2.51) а б в

Рис. 2.21. «Неоднородная» структура деформированной ГПУ-решетки:

а – верхняя грань после растяжения вдоль осей Ox1 и Ox2 ( 2 1 );

б – «элементарный треугольник» до деформации (шары уменьшены);

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

–  –  –

В зависимости от соотношения между 1 и 2 смещение может быть как в положительном направлении Ox2, так и в отрицательном направлении .

После наложения смещений (2.51) на слой B получается однородное распределение атомов деформированной ГПУ-решетки (рис. 2.23, б) .

–  –  –

Заметим, что при растяжении-сжатии вдоль оси Ox3 однородность ГПУ-решетки не нарушается и дополнительного смещения слоев B относительно слоев A не требуется .

Задания:

1. Построить алгоритм исследования упругих свойств для материала с ГПУ-решеткой .

2. Построить алгоритм исследования равновесных межатомных расстояний и упругих свойств графена и терморасширенного графита .

Симметрийные свойства ГПУ-решетки. Для планирования численных экспериментов по исследованию упругих свойств материала с ГПУ-решеткой определим число независимых ненулевых компонент тензора линейно-упругих свойств. Будем считать исследуемый материал с ГПУ-решеткой линейно-упругим по Коши, деформируя его однородно, ограничиваясь малыми деформациями. Отклик такого материала будем описывать несимметричным тензором напряжений, удовлетворяющим обобщенному закону Гука с несимметричными мерами: = C : u, T, Cijkl C jikl, Cijkl Cijlk, Cijkl = Cklij, где – оператор градиента из отсчетной конфигурации. Разложим получаемый отклик как функцию параметров накладываемой деформации (соответствующих компонент тензора дисторсии u ) в степенной ряд по этим параметрам. Коэффициенты при первых степенях компонент дисторсии будут соответствующими компонентами тензора С линейно-упругих свойств материала. Если материал является (нелинейно) упругим по Грину, то из разложения упругого потенциала в степенной ряд по дисторсии получается упругий закон также в виде ряда, первый член которого совпадает с правой частью записанного несимметричного обобщенного закона Гука .

Для определения числа независимых ненулевых компонент тензора С для материала с ГПУ-решеткой рассмотрим ее симметрийные свойства .

ГПУ-решетка обладает одной осью симметрии 3-го порядка и ортогональной ей плоскостью симметрии (или осью симметрии 2-го порядка, ортогональной к первой оси). Пусть первая ось совпадает с осью Ox3, а плоскость симметрии содержит оси Ox1 и Ox2. Если преобразование Q принадлежит группе симметрии кристаллической решетки, то справедливо соотношение

C = Cijkl (Q ei )(Q e j )(Q ek )(Q el ) = C, (2.52)

из которого следуют уравнения для определения ненулевых независимых компонент Cijkl. В случае ГПУ-решетки с учетом общих свойств Cijkl = Cklij и симметрии решетки получается 11 ненулевых независимых компонент:

C1111, C1112, C1122, C1133, C1212, C1233, C1332, C2323, C2332, C3232, C3333. (2.53)

–  –  –

При учете несимметрии тензора напряжений Коши появляются такие модули, как C1112 или C1332.

При наложении условия симметрии для этих компонент получим:

C1112 = C1121 = C1112, то есть C1112 = 0 .

–  –  –

11 = C1111 (u)11 + C1112 ( (u)21 (u)12 ) + C1122 (u)22 + C1133 (u)33, 12 = C1112 (u)11 + (C1111 C1122 C1212 )(u)12 + C1212 (u)21 + C1112 (u)22 + C1233 (u)33, = C (u) + C (u) + C (u),

–  –  –

Проводя опыты на растяжение вдоль различных из трех взаимно перпендикулярных осей, можно получить выражения для трех коэффициентов Пуассона. Модуль всестороннего растяжения-сжатия K связывает следы тензоров напряжений и деформаций sp( ) = 3K sp(u) и также может быть найден для ГПУ-решетки .

Рассмотрим опыты, необходимые для определения компонент (2.52) .

Пусть n – единичная нормаль к плоскости простого сдвига, m – единичный вектор, задающий направление сдвига. Закон движения в лагранжевой форме, описывающий однородную деформацию образца, имеет вид r = R + ( n R ) m, R – радиус-вектор центров атомов в отсчетной конфигурации, r – радиус-вектор атомов в текущей конфигурации. Тогда аффинор, описывающий простой сдвиг на величину,

–  –  –

Можно заметить, что недиагональная часть тензора напряжений получилась кососимметричной .

Таким образом, в рассматриваемом опыте можно определить 6 из 11 компонент тензора линейно-упругих свойств материала с ГПУ-решеткой – C1111, C1122, C1112, C1233, C1133 и C3333. Остается найти 5 компонент. Заметим, что проведенные численные расчеты показали, что получаемый при трехосном растяжении-сжатии тензор напряжений является диагональным, то есть должны быть справедливы равенства C1112 = 0, C1233 = 0. (2.64) Рассмотрим далее простой сдвиг образца в плоскости, перпендикулярной оси Ox1, в направлении Ox2: (u) 21 =, (u)ij = 0.

Рассмотрим также опыт на простой сдвиг, в котором (u)12 =, (u)ij = 0 :

–  –  –

22 = 23 = 31 = 32 = 33 = 0. 22 = 23 = 31 = 32 = 33 = 0 .

Таким образом, при таком деформировании можно найти еще 1 компоненту тензора линейно-упругих свойств материала с ГПУ-решеткой – C1212. Остается найти еще 4 компоненты .

Рассмотрим простой сдвиг образца в плоскости, перпендикулярной оси Ox1 в направлении Ox3: (u)31 =, (u)ij = 0. Рассмотрим также опыт на простой сдвиг, в котором (u)13 =, (u)ij = 0. Упругие законы соответственно принимают вид

–  –  –

Следовательно, в этих двух опытах можно найти оставшиеся 4 компоненты тензора линейно-упругих свойств материала с ГПУ-решеткой – C1332, C2323, C2332, C3232 .

При простом сдвиге образца в плоскости, перпендикулярной оси Ox3, в направлении Ox2 ( (u) 23 =, (u)ij = 0 ), а также в опыте на простой сдвиг, для которого (u)32 =, (u)ij = 0, упругий закон принимает вид (все компоненты уже известны, этот опыт может использоваться для их проверки)

–  –  –

Итак, с помощью 5 экспериментов можно найти все 11 независимых ненулевых компонент тензора линейно-упругих свойств материала с ГПУ-решеткой .

Задания:

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

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

3) Реализовать перечисленные опыты для получения значений упругих модулей образца с ГПУ-решеткой .

2.3. УЧЕТ ТЕМПЕРАТУРЫ. ЗАВИСИМОСТЬ МЕХАНИЧЕСКИХ

СВОЙСТВ КОНДЕНСИРОВАННЫХ СРЕД ОТ ТЕМПЕРАТУРЫ

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

Рассмотрим, как в предложенной модификации метода атомарной статики, позволяющей провести идентификацию параметров потенциала межатомного взаимодействия для кристаллических материалов с различного типа решетками (ГЦК-, ОЦК-, ГПУ-решетками графена и графита) по их макроскопическим параметрам – равновесному межатомному расстоянию и упругому модулю сдвига. В основе этой методики лежит рассмотрение статики взаимодействующих частиц при явном задании структуры кристаллической решетки. Возможности подхода были продемонстрированы на примере использования степенных потенциалов ЛеннардДжонса и Ми. Это обосновано тем, что при исследовании механических свойств нет необходимости рассматривать процессы, идущие при сверхнизких температурах или при скоростном деформировании, когда важную роль во взаимодействии атомов могут играть квантовые эффекты, учитываемые в потенциалах специального вида. Использование статической постановки позволяет получить аналитическое решение задачи определения равновесного межатомного расстояния для различных кристаллических решеток и образцов различного размера. При этом точные значения межатомного расстояния получаются для объемов материала с небольшим числом атомов N на ребре образца (от 3 до 20). Далее для перехода на макроуровень по этим точным решениям можно сделать предельный переход, устремив число атомов на ребре образца к бесконечности. В частности, это позволяет проводить идентификацию параметров потенциала по известным макроскопическим физико-механическим свойствам материалов .

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

При стабилизации системы энергия крупномасштабного движения атомов переходит в мелкомасштабные колебания с меньшей амплитудой (рис. 2.24). Под крупномасштабным движением здесь понимается изменение со временем положений атомов относительно центра масс системы, перемещения атомов при этом превышают среднее межатомное расстояние. Мелкомасштабные колебания характеризуются перемещениями атомов относительно некоторых стационарных центров на величины, значительно меньшие среднего расстояния между такими центрами. На рисунке показаны результаты расчета изменения кинетической энергии системы, отнесенной к числу атомов, при эволюции из начального состояния. Начальные условия соответствовали нулевыми скоростями, но увеличенным по сравнению с равновесным периодом решетки межатомным расстоянием. То есть в систему вносилась дополнительная внутренняя энергия, вначале – за счет потенциальной энергии, далее вследствие консервативности системы, потенциальная энергия переходит в кинетическую. Кинетическая энергия менялась в крупномасштабном движении за счет последовательности прохождения волн сжатия и расширения. Затем крупномасштабные колебания уменьшались при сохранении среднего значения кинетической энергии системы. На графике зависимости кинетической энергии от времени расчета (см. рис. 2.24) виден момент перехода энергии от крупномасштабного движения к колебаниям отдельных атомов вблизи некоторых фиксированных во времени положений. Эти мелкомасштабные колебания относительно не меняющихся со временем центров отождествляются с тепловым движением и с помощью выражения E = kT, где E – амплитуда пульсаций кинетической энергии, k – постоянная Больцмана, определяется температура тела T. В статическом подходе нет возможности вводить температуру таким образом. Можно было бы считать, что статический подход дает возможность определять положения центров, относительно которых атомы совершают колебательное движение, независимое от перемещений атомов, связанных с механическим движением или деформированием тела. Тогда температуру в статическом подходе можно учесть с помощью переопределения параметров потенциала введением в них коэффициента термического расширения .

Рис. 2.24. Зависимость удельной кинетической энергии системы атомов образца с ГЦК-решеткой (40 атомов на ребре) от числа шагов в процессе перехода к равновесному состоянию. Полная энергия системы сохраняется Можно пойти и по другому пути – имитировать тепловые колебания атомов с некоторой заданной амплитудой. Это можно реализовать наложением на систему атомов случайных смещений с этой амплитудой при равномерном распределении направления смещений в пространстве. Частота этих смещений (колебаний) в таком подходе рассматривается как независимый параметр, требующий идентификации для конкретного материала (подобно частоте Дебая). Для корректного задания отклонений атомов необходимо гарантировать равномерное распределение направления случайных смещений в пространстве. Направление задается единичным вектором, то есть задача сводится к нахождению равномерного распределения точек на сфере единичного радиуса .

Алгоритм получения равномерного распределения точек по сфере. В сферической системе координат равномерное разбиение интервалов для угла [0;2), откладываемого от оси Ox1 в плоскости x1Ox2 и угла [ / 2; / 2], откладываемого от плоскости x1Ox2, приводит к концентрации точек на полюсах сферы (рис. 2.25, а). Связь декартовых и сферических координат имеет вид: x1 = r cos cos, x2 = r cos sin, x3 = r sin .

Для того чтобы концентрации не происходило, разбивать необходимо не интервал угла [ / 2; / 2], а интервал [1;1] значений функции sin .

Функция cos выражается через полученное случайное значение как 1 sin 2. В результате будет сгенерировано равномерное покрытие сферы точками (рис. 2.25, б) .

–  –  –

“Show[Graphics3D[Point/@positions1]], AxesTrue, TicksFalse, AxesLabel {"X1", "X2", "X3"}]” .

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

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

В пакете «Mathematica» заложена полная информация о большом количестве геометрических тел, к которой можно получить доступ с помощью функции “PolyhedronData[]”. В частности, для изображения икосаэдра “PolyhedronData["Icosahedron"]” .

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

Для этого сформируем список vertexes координат вершин икосаэдра:

“vertexes = PolyhedronData["Icosahedron", "VertexCoordinates"]” .

Эти координаты представлены для икосаэдра с единичным ребром .

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

“nodes = N[ #/ #.# ] & /@ vertexes;” .

Команда “PolyhedronData["Icosahedron", "Faces"]” возвращает список, состоящий из двух частей. Первая часть содержит список вершин, а вторая имеет структуру “Polygon[faces]”, где список faces состоит из троек натуральных чисел, каждая из которых хранит номера вершин (из первого списка), образующих ту или иную грань. Длина списка faces равна 20. Команда “facesSets = PolyhedronData["Icosahedron", "Faces"][[2, 1]]” сохраняет набор троек вершин граней в списке facesSets .

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

“pole = (sin=Random[Real, {–1, 1}]; = Random[Real, {0, 2}]; r=1;

–  –  –

Повернем каждую вершину икосаэдра на сфере на случайный угол u вокруг оси выбранного случайного полюса pole и сохраним изменения в том же списке:

“nodes = (RotationMatrix[u, pole].#) &/@ nodes) /. u Random[Real, {0, 2}];” .

Изобразим вершины полученного случайно развернутого икосаэдра на единичной сфере красными точками:

“mainPointsPlot=Show[Graphics3D[{PointSize[Large], Red,Point[#]}& /@ nodes], Graphics3D[{Opacity[0.3], Sphere[{0, 0, 0}, r]}], Axes True, Ticks False, AxesLabel {"X1", "X2", "X3"}]”, где функция “Point[]” рисует точку с заданными координатами, функция “Sphere[]” дает изображение сферы заданного радиуса r (в примере r = 1 ) с центром в заданной точке, функция “ PointSize[]” определяет размер точки при рисовании, “ Opacity[]” задает прозрачность изображаемой поверхности .

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

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

“pointsInFace = Table[{},{i,20}]; s = 1.05146;

Do[If[Max[ (positions[[ j]] #).(positions[[ j]] #) &/@nodes[[facesSets[[i]] ]]] s, pointsInFace[[i]]=Append[pointsInFace[[i]], positions[[j]] ] ], {i,20}, {j, Length[positions]}]; ” .

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

“oneFacePlot = Show[Graphics3D[{PointSize[Large], Blue, Point[#]}& /@ nodes[[facesSets[[5]] ]] ] ];

pointsInFacePlot = Graphics3D[{PointSize[Large], Orange, Point[#]}& /@ pointsInFace[[5]] ];

Show[randomPointsPlot, mainPointsPlot, oneFacePlot, pointsInFacePlot]” .

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

Многократно повторяя в цикле построенную процедуру и сохраняя в списке distributionList число точек на каждом из 20 треугольников при каждом их случайном размещении на сфере “distributionList = Join[distributionList, Length /@ pointsInFace]”, получим статистику попадания точек на сфере в правильные криволинейные треугольники, полностью покрывающие сферу .

–  –  –

В результате для приведенных на рис. 2.25 распределений 5000 точек на сфере получена статистика их попаданий в различные треугольники (рис. 2.27). При равномерном распределении в каждый из 20 треугольников должно попадать примерно 250 точек .

–  –  –

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

Для анализа полученных распределений построим с помощью пакета «Mathematica» гистограммы функции распределения вероятностей (один столбец диаграммы строится на 5 точках – второй аргумент функции):

“Histogram[{distributionList1, distributionList2}, {5}, "Probability", PlotRange {{0, Max[distributionListI] + 10}, All}]” .

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

Рис. 2.28. Распределение вероятности для двух способов размещения точек по треугольникам: левое распределение соответствует рис. 2.25, а;

правое распределение соответствует рис. 2.25, б

–  –  –

“tensorP = Table[0, {3}, {3}];

Do[tensorP+=Outer[Times,positions1[[i]],positions1[[i]]], {i, Length[positions1]}]” .

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

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

Рис. 2.29. Неравномерное распределение точек по сфере и собственные векторы соответствующего тензора P

Построить рис. 2.29 позволяет команда:

“Show[randomPointsPlot1, Graphics3D[{Opacity[0.3], Sphere[{0, 0, 0}, r]}], Graphics3D[{Green, Arrow[Tube[{{0,0,0},1.4#},0,02]]& /@ Eigenvectors[tensorP/Length[positions1]]}], Axes True, Ticks False, AxesLabel {"X1", "X2", "X3"}}]” .

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

Для материалов с кубической решеткой (ГЦК, ОЦК) исследуется кубический образец. Для материала с ГПУ-решеткой берется прямая призма с правильным треугольником в основании (см. рис. 2.17), симметрия которого соответствует симметрии ГПУ-решетки .

Для определения зависимости равновесного межатомного расстояния a* монокристалла в форме куба с ГЦК- или ОЦК-решеткой от амплитуды возмущений A (аналога температуры) для каждого заданного значения A при произвольных переменных a (параметр решетки), и (параметры потенциала Леннард-Джонса) требуется вычислить силы на гранях куба по описанному ранее алгоритму. Далее определяется минимум суммы квадратов всех найденных 6 сил. Минимум такой суммы не всегда равен нулю, поскольку для каждого отдельного распределения атомов не всегда выполняется условие равновесия, справедливое для сил, полученных осреднением по набору реализаций возмущенной конфигурации атомов (как и для произвольного момента времени в методе МД). Сила на грани образца равна (C1(i ) (n, A) 6 C2i ) (n, A) a 6 ) / a13, где индекс i = 1,6 (

–  –  –

Для каждого значения амплитуды A строились 10000 реализаций случайных конфигураций решетки, и по ним с помощью (2.69) определялось среднее значение межатомного расстояния. Также для каждой реализации при расстоянии a = a* определялась величина U потенциальной энергии всей системы атомов, отнесенной к числу атомов .



Pages:   || 2 |



Похожие работы:

«АВТОМОБИЛЬНЫЙ ВИДЕОРЕГИСТРАТОР AV-604 / РУКОВОДСТВО ПОЛЬЗОВАТЕЛЯ / СОДЕРЖАНИЕ ОБЩАЯ ИНФОРМАЦИЯ ВАЖНАЯ ИНФОРМАЦИЯ СХЕМА УСТРОЙСТВА НАЧАЛО РАБОТЫ ЗАРЯДКА АККУМУЛЯТОРА УСТАНОВКА ВИДЕОРЕГИСТРАТОРА УСТАНОВКА ВЫНОСНОЙ КАМЕРЫ ВКЛЮЧЕНИЕ /ОТКЛЮЧЕНИЕ УПРАВЛЕНИЕ ФУНКЦИЕЙ МЕНЮ НАСТ...»

«МЕЖГОСУДАРСТВЕННЫЙ СОВЕТ ПО СТАНДАРТИЗАЦИИ, МЕТРОЛОГИИ И СЕРТИФИКАЦИИ (МГС) INTERSTATE COUNCIL FOR STANDARDIZATION, METROLOGY AND CERTIFICATION (ISC) ГОСТ МЕЖГОСУДАРСТВЕННЫЙ 13047. 5СТАНДАРТ НИКЕЛЬ. КОБАЛЬТ Методы определения никеля в кобальте Издание оф ициальное М осква С т...»

«Министерство образования и науки Российской Федерации федеральное государственное автономное образовательное учреждение высшего образования "НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ТОМСКИЙ ПОЛИТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ" Инженерная школа природных ресурсов Направ...»

«Государственное управление. Электронный вестник Выпуск № 37. Апрель 2013 г. Стремоухова А.Д. Механизмы развития инноваций на базе высших учебных заведений 1. Роль высших учебных заведений в национальных инновационных системах Проблема развития инноваций является одной из самых актуальных, так как именно сегодня набл...»

«Dell Vostro 3470 Руководство по настройке и техническим характеристикам нормативная модель: D13S нормативный тип: D13S003 Примечания, предостережения и предупреждения ПРИМЕЧАНИЕ: Пометка ПРИМЕЧАНИЕ указыв...»

«TAJIKISTAN INDUSTRIAL MODERNIZATION AND COMPETITIVENESS IMPROVEMENT OF CARPET WEAVING, EMBROIDERY AND TEXTILE SECTORS SECOND COLLECTION OF UNIDO’S PROJECT PILOT BENEFICIARY ENTERPRISES ВТОРАЯ ОБЪЕДИНЕННАЯ КОЛЛЕКЦИЯ ПРЕДПРИЯТИЙ-БЕНЕФИЦИАРОВ ПРОЕКТА ЮНИДО Republic of Tajikistan UNIDO PROJECT: INDUSTRAL MODERN...»

«ИСО "Орион" Устройство оконечное системы передачи извещений по каналам сотовой связи GSM "УО-4С исп.02" Руководство по эксплуатации ЗАО НВП "Болид", Россия, 141070, Московская область, г. Королёв, ул. Пионерская, д. 4...»

«Министерство науки и высшего образования Российской Федерации Федеральное государственное автономное образовательное учреждение высшего образования "НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ТОМСКИЙ ПОЛИТЕХНИЧЕСКИЙ УНИВ...»

«Приложение № 3 к Договору на выпуск и обслуживание расчетных корпоративных банковских карт Банка "ВБРР" (АО) № от "" _ 20_ г. УСЛОВИЯ ВЫПУСКА, ОБСЛУЖИВАНИЯ И ИСПОЛЬЗОВАНИЯ РАСЧЕТНЫХ КОРПОРАТИВНЫХ БАНКОВСКИХ КАРТ БАНКА "ВБРР" (АО) 1. ОБЩИЕ ПОЛОЖЕНИЯ 1.1. Настоящие Условия определяют порядок выпуска,...»

«По всем вопросам обращайтесь в офис компании Ти-Системс Телефоны для связи:+7 (495) 777-4-788, 5007154, 55, 65, 7489127, 28, 29 Web-сайт: www.tisys.ru E-Mail: info@tisys.ru Спринклеры ELO ТЕХНИЧ...»

«ПАСПОРТ БЕЗОПАСНОСТИ ХИМИЧЕСКОЙ ПРОДУКЦИИ (8аГе1у БаСа 8Ьее1) Внесен в Регистр & РПБ № 1 1 1 1^ ** -. ^ / 1 •—1 ! — $,:,,,3* -1 '''{гг " "^4^420 11$ 1 & 5 З ' У О ф И Л Г * " 20/$ Г. /* ? // ИАЦ _ —1 Г Ъ *1 / РосстанЦа^т 4 V \ Л |. Г А II | / г ъ? / Информационно-аналип ический центр д ВЯ *Щ I,я гу п А Йк И "Безопасность веществ 1(м атериалов" Руководитель ч /А. А...»

«Круглов Евгений Владимирович Микроконтроллерная система управления сканирующим зондовым микроскопом НаноСкан 05.13.05 — Элементы и устройства вычислительной техники и систем управления Автореферат диссертации на соискание ученой степени кандидата технических наук Москва 2009 Работа...»

«МЕЖГОСУДАРСТВЕННЫЙ СОВЕТ ПО СТАНДАРТИЗАЦИИ, МЕТРОЛОГИИ И СЕРТИФИКАЦИИ (МГС) INTERSTATE COUNCIL FOR STANDARDIZATION, METROLOGY AND CERTIFICATION (ISC) ГОСТ МЕЖГОСУДАРСТВЕННЫЙ ISO 9 5 2 6 СТАНДАРТ ФРУКТЫ, ОВОЩИ И ПРОДУКТЫ ИХ ПЕРЕРАБОТКИ Определение содержания железа методом пламенной атомно-абсор...»

«АГЕНТСТВО 12 НОЯБРЯ МАКСИМОВ–КОНСАЛТИНГ 2015 ГОДА www.maksimov.info АНАЛИТИЧЕСКИЙ МАТЕРИАЛ "НОВЫЕ ВОЗМОЖНОСТИ НОВЫХ ПАРТИЙ" ВВЕДЕНИЕ В предыдущем докладе "Шансы и Техники" мы оценили потенциал непарламентс...»

«ФЕДЕРАЛЬНОЕ АГЕНТСТВО НО ТЕХНИЧЕСКОМУ РЕГУЛИРОВАНИЮ И МЕТРОЛОГИИ ФЕДЕРАЛЬНОЕ БЮДЖЕТНОЕ УЧРЕЖДЕНИЕ "ГОСУДАРСТВЕННЫЙ РЕГИОНАЛЬНЫЙ ЦЕНТР СТАНДАРТИЗАЦИИ, МЕТРОЛОГИИ И ИСПЫТАНИЙ В Г. МОСКВЕ" (ФБУ "РОСТЕСТ МОСКВА") УТВЕРЖДАЮ Заместитель генерального директора "Ростест-Москва" у А.Д. Меньшиков "2...»

«КАТАЛОГ РАЗРАБОТОК И УСЛУГ никимт-! АТОМСТРОИ Открытое акционерное общество 9 Действующие филиалы О Перспективные филиалы 4 Сосновый бор ОУдомля ОМосква ООбнинск Курчато Десногорск ОНововоронеж ОЗаречный дБелоярск ОСеверск °0зерск ОТомск Волгодонск 28 ОЖелезногорск СОДЕРЖАНИЕ Головная роль ОАО "НИНИМТ-Атомстрой" 2 Монтаж, демонтаж,...»

«Министерство образования и науки Российской Федерации Федеральное государственное автономное образовательное учреждение высшего образования "Уральский федеральный университет имени первого Президента Рос...»

«International Paint Ltd. Справочный Лист Безопасности QZB000 INTERCRYL 525 WHITE Номер редакции документа 4 Дата Последней Редакции 26/06/14 Соответствует требованиям Директивы (EC) No.1907/2006 (REACH), Приложения II. РАЗДЕЛ 1: Идентифик...»

«МРНТИ: 67.11.29 Аскар Жагпарович Жусупбеков Карлыгаш Боранбайкызы Боргекова Евразийский национальный университет имени Л.Н. Гумилева, Астана, Казахстан ГЕОТЕХНИЧЕСКОЕ СТРОИТЕЛЬСТВО И ТЕСТИРОВАНИЕ СВАЙНЫХ ФУНДАМЕНТОВ В СЛОЖНЫХ ГРУНТОВЫХ УСЛОВИЯХ КАЗАХСТАНА Аннотация. В статье приведены результаты исследования определения...»

«ГОСТ Р 51924-2002 Г О С У Д А Р С Т В Е Н Н Ы Й СТАНДАРТ Р О С С И Й С К О Й Ф Е Д Е Р А Ц И И ТРУБЫ ДВО Й Н Ы Е к о л о н к о в ы е ДЛЯ ГЕОЛОГО-РАЗВЕДОЧНОГО БУРЕНИЯ Общие те...»

«© www.complexs.ru Комплекс – С ГОСУДАРСТВЕННЫЙ СТАНДАРТ СОЮЗА ССР ФЕРМЫ ЖЕЛЕЗОБЕТОННЫЕ ТЕХНИЧЕСКИЕ УСЛОВИЯ ГОСТ 20213-89 Издание официальное (343) 378-02-78 г.Екатеринбург © www.complexs.ru Комплекс – С УДК 624.072.2:006.354 группа ЖЗЗ Г О С У Д А Р С Т В Е Н НЫ...»







 
2019 www.librus.dobrota.biz - «Бесплатная электронная библиотека - собрание публикаций»

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