Численные процедуры решения геометрически нелинейных задач механики деформируемого твердого тела Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 14, № 4, 2009

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

В. Д. Кургузов Учреждение Российской академии наук

Институт гидродинамики им. М.А. Лаврентьева СО РАН, Новосибирск, Россия

e-mail: [email protected]

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

Ключевые слова: геометрическая нелинейность, метод конечных элементов.

Введение

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

В настоящее время на рынке имеется довольно широкий выбор пакетов конечноэлементного анализа, например, зарубежные Ansys, MSC.Nastran, MSC.Marc, LS-Dyna,

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 07-01-00163).

© ИВТ СО РАН, 2009.

Abacus, отечественные Scad, WinMachine и др. Наличие пре- и постпроцессоров, полностью управляемый графический интерфейс, гибкая и автоматизированная генерация конечно-элементной сетки, увязка нагрузочных и граничных условий с геометрией, изоповерхности и другие средства визуализации — все это помогает ускорить получение результатов и облегчить работу конструкторов и расчетчиков. Однако у молодых инженеров наблюдается некая эйфория при использовании пакетов конечно-элементного анализа: “Дайте мне компьютер помощнее, и я решу любую задачу”. Проблема заключается не в получении какого-либо решения вообще, а в том, чтобы убедиться в достоверности полученного решения, а это особенно важно при решении нелинейных задач. Необходим качественный анализ поведения конструкции, проведение ряда тестовых расчетов, апробация расчетной методики на модельных задачах, сравнение с результатами решений, полученных другими методами и, в идеале, сравнение с экспериментом, что, конечно, не всегда бывает возможным.

Многие задачи механики деформируемого твердого тела сводятся к решению систем нелинейных алгебраических, трансцендентных или дифференциальных уравнений, содержащих параметр. Трудности, возникающие в геометрически нелинейных системах, видны уже при анализе процесса деформирования простейшей двухстержневой системы — фермы Мизеса [1], где зависимость сила—смещение является нелинейной. Еще большие трудности возникают при решении задач деформирования таких систем, как пологие арки и пологие оболочки. Например, для пологой арки связь между нагрузкой и усилием распора имеет вид петлеобразной кривой [2], аналогичные особенности имеют решения для пологих оболочек [3].

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

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

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

Предполагаем, что материал стержней является линейно-упругим и все стержни работают в условиях чистого растяжения-сжатия. Следовательно, при использовании МКЭ функция формы первого порядка точно описывает поле перемещений и для моделирования каждого стержня достаточно использовать по одному одномерному двухузловому конечному элементу. Внутренняя сила Д, возникающая в стержневом элементе при его деформации, определяется так: f = се, е = (I — 10)/10. Здесь с — жесткость стержневого элемента (с = ЕА, где Е — модуль Юнга, А — площадь поперечного сечения), 1о — начальная длина стержня, I — длина стержня в деформированном состоянии. В дальнейшем нам понадобится величина Д = дf/дl = с/\0.

Обозначим через ие вектор перемещений стержневого элемента:

ие = [п\ п1 и\ п\ п1 п\]т.

Здесь П — ]-я (] = 1, 2, 3) компонента вектора перемещений г-го (г = 1,2) узла стержневого элемента. Прямыми вычислениями можно показать выполнение равенства

д1 Вт,

д И'

где В = [—в\ — в2 — в3 в\ в2 в3], вк = 1к/I (к = 1, 2, 3) являются компонентами

вектора единичной длины е = 1/1, 1 = х2 — XI, I = \/1тТ — длина стержневого элемента, XI, х2 — радиус-векторы узлов стержневого элемента 1, 2 соответственно. Вектор внутренних сил элемента Ее теперь можно представить следующим образом Ее = Д Вт. С использованием стандартной процедуры метода конечных элементов [7] вектор внутренних сил стержневой конструкции Е определяется ассемблированием векторов внутренних сил Ее (1 < в < М) всех элементов, составляющих стержневую конструкцию.

Введем касательную матрицу жесткости Ке стержневого элемента [8]:

Ке = /іВТВ + у(Р - БТБ), (1)

где

Р

I -I I I

1 0 0 0 1 0 0 0 1

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

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

Линеаризованные уравнения квазистатического деформирования имеют вид [6]

* ки = *+А*К — *Е, (2)

где И = *+А*И — * И — вектор приращений перемещений, *+А*И, — вектор внешних сил, *Е — вектор внутренних сил. Вслед за [7] здесь и далее левый верхний индекс обозначает момент времени, в который рассматривается данная величина. Под временем в квазистатических задачах понимается некоторый монотонно возрастающий параметр нагружения.

I

2. Процедуры численного решения уравнений квазистатического деформирования

Для решения системы уравнений (2) воспользуемся алгоритмом пошагового интегрирования уравнений равновесия в сочетании с итерационной процедурой Ньютона—Рафсона

[6]. Наиболее известным методом итерационного уточнения решения нелинейных задач является стандартный метод Ньютона—Рафсона. В момент времени £ + Д£ итерационная процедура сводится к решению последовательности систем уравнений

*+А*К(*-1) ди<*) = *+А*Л — *+а*е(*-1) (г =1 2 ) (3)

где верхний правый индекс обозначает номер итерации, ДИ(г) = *+А*И(г) — *+А*И(г-1) — поправка вектора приращений перемещений, а начальные значения для итерационного процесса следующие:

*+А*к(о) *к *+а*и(о) *и *+а*е(о) *е

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

При мягком нагружении вектор внешних сил *+А*Л известен и решение системы уравнений (2) можно получить, используя итерационную процедуру (3). Сложнее обстоит дело, когда в решении задачи встречаются собственные состояния типа максимальной нагрузки. Трудность состоит в том, что в качестве монотонно возрастающего параметра деформирования нельзя взять внешнюю силу, действующую на конструкцию. Например, при Р > Ртах, где Ртах — максимальная нагрузка, решения задачи не существует. Для того чтобы обойти эту трудность, введем параметр внешней силы в число неизвестных величин (наряду с перемещениями), а в качестве дополнительного (контрольного) уравнения к системе (3) зададим длину дуги в (И, Л)-пространстве (И — вектор перемещений, Л — параметр внешней силы). Ниже используется метод Крисфилда с движением вектора неизвестных по сфере в (И, Л)-пространстве [9].

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

*+А*^ *+А* л тз

Л — Л Ло,

где Ло — некоторый постоянный вектор, а параметр *+А*Л характеризует интенсивность действия внешней силы и предполагается неизвестным. На первом шаге интегрирования уравнений равновесия задачу будем решать так, чтобы некоторая компонента вектора перемещений принимала заданную величину и*, т. е. получаем контрольное уравнение для первого шага по времени

А*и = и *,

где А*ик — некоторая компонента вектора перемещений на первом шаге по времени. На последующих шагах по времени в качестве контрольного уравнения используется условие движения в итерационном процессе по сфере с постоянным радиусом Д1 (приращение длины дуги интегральной кривой Д1 в (И, Л)-пространстве предполагается заданной величиной):

ф2 Л«2 + Имт Ии = Д12, где Л« = *+А*Л« — *Л, 0 < ф < 1.

Для получения начальных значений решаем систему уравнений [9]

*кИ = Ло

1Ю-

Приращение вектора перемещений И(1) находится из равенства И(1) = л(1) II, где л(1) определяется из контрольного уравнения: в момент времени £ = Д£ (первый шаг) Л(1) = и*/С/к; в моменты времени £ = 2Д£, 3Д£, ... (последующие шаги по времени)

Л(1) = ±

А/

иТ и + ф2

(4)

Для итерационного уточнения решения на каждой итерации решаются две системы уравнений [9]:

*+а*К(*-1) ДИ (*) = *+А*Л(*-1)И,о — *+а*е(*-1) *+а*К(*-1) ДИ м = Ло.

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

ДИИ = ДИ ^ + ДЛМ ДИ

Здесь ДЛ(г) определяется из контрольного уравнения в момент времени £ = Д£ (первый шаг):

ДЛ(г) = — Дикг)/ Дикг);

в моменты времени £ = 2Д£, 3Д£, . .. (последующие шаги по времени):

АЛ»

-Ь ± s/W~-

ас

(5)

где а = ф2 + ДИ «т ДИ«, Ь = ф2Л(*-1) + И(*-1)т ДИ«, с = ф2Л(*-1)2 + 'И(*-1)т И(*-1) — Д12. Здесь введен вектор И(г-1) = И(г-1) +ДИ(г). Корни в равенствах (4), (5) выбираются так, чтобы обеспечивалось максимальное значение величины 7 = И(г-1)т И(г). Условие максимальности 7 соответствует большей гладкости пути решения в (И, Л)-пространстве.

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

На нулевой итерации решается система уравнений

*ки(0) = к(0).

(6)

Здесь Ь компонент вектора приращений смещений и(0) с номерами к1, к2, ... &£ равны заданным значениям , Ц;2,... , , а остальные неизвестны. Соответственно Ь

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

' К* К„ь ' ' и„ ' Ка

Кьа Кьь . иь .

а

где Иа — неизвестные смещения, И,а — заданные силы, Иь — заданные смещения, Ль — неизвестные силы. Умножая первую строку на столбец, получим КааИа + К„ьИь = И,а,

откуда Кааиа = И,а — К„ьи&. Решая эту систему уравнений, находим вектор иа. Умножая вторую строку на столбец, получим Кь„иа + К^Ц = И*, откуда находим вектор Иь Вычисляем матрицу жесткости *+А*К(0) и вектор внутренних сил *+А*р(0), находим вектор *+А*И(0) = + И(0).

На последующих итерациях решаем систему уравнений

*+А*К(і-1) ДЦ(І) = *+А*к,(і) — *+А*р(і-і) (і =1, 2, ...), (7)

в которой Ь компонент вектора ДЦ(і) с номерами к1, к2, ... равны нулю, а остальные неизвестны. Соответственно Ь компонент вектора *+А*И(г) с теми же номерами неизвестны, а остальные равны нулю. В результате решения находим N — Ь компонент вектора Ди(і) и Ь компонент вектора *+А*И(г). Вычисляем вектор Ц(і) = Ц(і-1) + ДЦ(і) и проверяем выполнение критериев сходимости.

3. Численный пример

Рассмотрим решение задачи о деформировании простой стержневой конструкции (рис. 1), которая несколько лет назад вызвала оживленную дискуссию в сети Интернет и на страницах специализированных изданий [10-13].

Система состоит из двух линейно-упругих стержней 1 и 2, испытывающих деформации растяжения-сжатия. Под действием силы Р = 95 Н узел 1 скользит по вертикальной направляющей, вовлекая в движение узел 2, который может перемещаться только в горизонтальном направлении, узел 3 шарнирно оперт. Жесткости стержней с1 = 1000 Н, с2 = 2000 Н, геометрические размеры (в метрах) показаны на рис. 1. Хотя материал стержней подчиняется закону Гука, вследствие наличия больших смещений и поворотов задача становится геометрически нелинейной.

Проведем численный анализ решения этой задачи, используя изложенный выше алгоритм (мягкое нагружение). На рис. 2 представлен график зависимости силы Р от смещения первого узла и1. Показаны все возможные равновесные состояния системы до уровня нагружения 220 Н. Как видно из рисунка, прямая Р = 95 пересекает график

Рис. 1. Простая конструкция из двух линейно-упругих стержней

в трех точках, т. е. решение задачи не единственно. На рис. 3 показаны три равновесные конфигурации системы, отвечающие внешней нагрузке 95 Н: две конфигурации — стержень 2 прощелкивается вправо, стержень 1 — вниз (рис. 3, а, б), также оказалось, что возможно равновесное состояние без прощелкивания стержня 2 (рис. 3, в).

Д. Назаров в своих публикациях [10, 12, 13] утверждает, что левое положение узла 2 (рис. 3, в) “неправильное”, а “правильное” только правое положение (рис. 3, а). Из приведенного анализа видно, что все три положения являются правильными. Неоднозначность решения, невозможная при линейном анализе, — обычное явление при анализе нелинейном. Кроме того, решение, приведенное Д. Назаровым в [13], неверно. В этом легко убедиться, подставив найденные Д. Назаровым положения узлов 1 и 2 (у = 22.074 для узла 1, x = 13.02 для узла 2) в уравнения равновесия. На рис. 4 показаны внутренние усилия в стержнях 1 и 2 в одном из равновесных положений. Проведя несложный качественный анализ, можно убедиться, что стержень 1 растянут, а стержень 2 сжат. Условия равновесия всех сил, действующих на узлы 1 и 2 (рис. 4), имеют вид

P = N cos а, N cos в + N2 cos 7 = 0, (8)

P

Рис. 2. Зависимость силы Р от смещения первого узла щ при мягком нагружении

14

1\

а

в

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

Рис. 4. Равновесие сил, действующих на узлы конструкции

где P — сила внешнего воздействия на систему; N и N2 — внутренние силы, действующие в стержнях 1 и 2 соответственно;

у — 10 x x — 12

cos а = —-----, cos в = т“, cos 7 = —--;

/1 /1 /2

/1 = Уx2 + (у — 10)2, /2 = У(x — 12)2 + 22 — длины деформированных стержней 1 и 2

соответственно; N1 = c1e1, N2 = c2e2; e1 = (/1 — /01)//01, £2 = (/2 — /о2)//o2 —относительные

удлинения стержней 1 и 2 соответственно; /01 = 10\/2, /02 = 2\/2 — первоначальные длины стержней; x — положение узла 2, у — положение узла 1.

Система (8) представляет собой систему двух нелинейных алгебраических уравнений относительно двух неизвестных x и у. Автору не известно аналитическое решение системы (8). При P = 95 решение системы (8) может быть найдено с использованием компьютерной программы аналитических вычислений Mathematica 6.0 с любой степенью точности, если в качестве начальных условий итерационного процесса взять данные, полученные из численного решения (см. рис. 2). Кроме того, в правильности численного решения можно удостовериться прямой подстановкой найденных значений в уравнения равновесия (8). Заметим, что решение, предложенное Д. Назаровым в [13], x = 13.02, у = 22.074 не удовлетворяет уравнениям равновесия (8).

В таблице приведены результаты расчетов по алгоритму (3)-(5) (u1, u2 — перемещения узлов 1 и 2 соответственно, P — значение приложенной силы) и аналитическое решение (x, у) по программе Mathematica 6.0. Величины u1, u2 связаны с x, у соотношениями u1 = у, u2 = x — 10.

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

П\ П2 Р X У Стержень 1 Стержень 2

19.444551 3.516602 95.063244 13.5169 19.4422 Растянут Сжат

20.449667 2.414533 94.935944 12.4148 20.4513 » »

22.457239 -0.197498 94.991328 9.8025 22.4574 » Растянут

Р

Рис. 5. Зависимость силы Р от смещения первого узла П\ при жестком нагружении

и

конструкции при жестком нагружении. Зададим перемещение узла 1 (см. рис. 1) по вертикали и пройдем весь диапазон перемещений от 0 до 24 с достаточно мелким шагом (например, 0.001). На каждом шаге решения задачи при фиксированной величине перемещения узла 1 определяется нагрузка Р. График зависимости Р от « представлен на рис. 5. Сравнивая график на рис. 2 (мягкое нагружение) с графиком на рис. 5 (жесткое нагружение), можно обнаружить, что решение, представленное на рис. 5, “обрезает” решение, показанное на рис. 2. Это означает, что не все равновесные состояния системы на рис. 2 реализуются физически, например, на рис. 2 имеются участки, где внешняя сила Р отрицательна. Отметим, что для анализа устойчивости каждой из равновесных форм требуется дополнительное исследование. Другое исследование, которое необходимо провести для данной задачи, — проверка наличия на интегральной кривой точек бифуркации. Критерием бифуркации решения является смена знака диагональных элементов касательной матрицы жесткости на промежутке от £ до £ + Д£. Проверка упомянутого критерия на каждом шаге интегрирования показала, что в данной задаче таких точек нет.

Следует также отметить резкое замедление скорости сходимости итерационного процесса в точках разрыва функции Р = Р(«1), график которой приведен на рис. 5. Численный алгоритм делает в указанных точках 514 и 565 итераций соответственно, для сравнения — во всех остальных точках итерационная процедура Ньютона—Рафсона сходится за 2-4 итерации при относительной точности по смещениям и силам еи = е/ = 10-5, по энергии деформирования ее = 10-10.

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

4. Приложение к задачам наномеханики

Предложенные выше численные алгоритмы могут быть с успехом применены к решению задач о деформировании и потере устойчивости атомных решеток [8, 14-16]. Характер деформирования атомной решетки близок к характеру деформирования стержневой конструкции (фермы): атомы решетки можно рассматривать как узлы фермы, а отрезки, соединяющие атомы, — как стержни с нелинейным поведением. На каждый находящийся в решетке атом действуют внешние силы и силы со стороны ближайших соседей. Поскольку размер конечного элемента согласован с постоянной атомной решетки, решается задача механики деформируемого твердого тела со структурой, причем взаимодействие элементов структуры определяется действующими физическими потенциалами. Уравнения движения атомной решетки следуют из принципа возможных перемещений: виртуальная работа внутренних сил равна виртуальной работе внешних сил и сил инерции. Для интегрирования уравнений движения используется пошаговая процедура (3)-(5), (6), (7): равновесная конфигурация и внутренние силы атомной решетки в момент £ предполагаются известными и решается задача по определению равновесной конфигурации и внутренних сил в момент времени £ + Д£. При достаточно малом шаге интегрирования Д£ решение нелинейной системы уравнений сводится к решению линеаризованной системы уравнений. На каждом шаге по времени это решение уточняется с помощью итерационной процедуры Ньютона—Рафсона, затем происходит переход на следующий шаг.

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

Заключение

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

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

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

[1] Клюшников В.Д. Лекции по устойчивости деформируемых систем. М.: Изд-во МГУ, 1986. 224 с.

[2] Григолюк Э.И. К вопросу устойчивости пологих арок // Инженерный сборник. 1951. Т. 9. С. 178-201.

[3] Вольмир А.С. Гибкие пластины и оболочки. М.: Гостехиздат, 1956. 419 с.

[4] Григолюк Э.И., Шалашилин В.И. Проблемы нелинейного деформирования. М.: Наука, 1988. 232 с.

[5] Шалашилин В.И., Кузнецов Е.Б. Метод продолжения решения по параметру и наилучшая параметризация в прикладной математике и механике. М.: Эдиториал УРСС,

1999. 222 с.

[6] Коробейников С.Н. Нелинейное деформирование твердых тел. Новосибирск: Изд-во СО РАН, 2000. 300 с.

[7] Бате К., Вилсон Е. Численные методы анализа и метод конечных элементов. М.: Строй-издат, 1982. 448 с.

[8] Korobeynikov S.N. Nonlinear equations of deformation of atomic lattices // Archives of Mechanics. 2005. Vol. 57, N 6. P. 435-453.

[9] Bathe K.-J., Dvorkin E.N. On the automatic solution of nonlinear finite element equations // Computers & Structures. 1983. Vol. 17. P. 871-879.

[10] Назаров Д. Обзор современных программ конечно-элементного анализа // САПР и графика. 2000. № 2. С. 52-55.

[11] Данилин А., Зуев Н., СнЕговский Д., Шалашилин В. Об использовании метода конечных элементов при решении геометрически нелинейных задач // САПР и графика.

2000. № 4. С. 26-31.

[12] http://n-t.ru/tp/ts/mka.htm

[13] http://www.cad.dp.ua/obzors/CAE.php

[14] Аннин Б.Д., Коробейников С.Н., Бабичев А.В. Компьютерное моделирование выпучивания нанотрубки при кручении // СибЖИМ. 2008. Т. 11, № 1. С. 3-22.

[15] Корнев В.М., Кургузов В.Д. Моделирование процесса движения краевых дислокаций методом конечных элементов // Физическая мезомеханика. 2008. Т. 11, № 5. С. 25-38.

[16] Kornev V.M., Kurguzov V.D. Multiparametric sufficient criterion of quasi-brittle fracture for complicated stress state // Engineering Fracture Mechanics. 2008. Vol. 75, N 5. P. 1099-1113.

Поступила в редакцию 9 февраля 2009 г.