Об использовании метода конечных элементов при решении геометрически нелинейных задач

Александр Данилин, Николай Зуев, Дмитрий Снеговский, Владимир Шалашилин

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

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

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

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

Уравнения равновесия сил, действующих на узлы A и B (рис. 2), соответственно, имеют вид P=N₁ cos a и N₁ cos b=N₂ sin g, где P — сила внешнего воздействия на упругий механизм, а N₁ и N₂ — внутренние силы, действующие в поперечных сечениях первого (левого) и второго (правого) стержней. Эти уравнения необходимо дополнить соотношениями, определяющими геометрию системы: sin a=x/l₁ , cos b=(b–y)/l₁ , sin g=(a–x)/l₂. Здесь , — длины деформированных стержней. Наконец, необходимо установить физические соотношения между силами N₁, N₂ и относительными деформациями стержней e₁= (l₀₁ –l₁)/l₀₁, e₂= (l₀₂ –l₂)/l₀₂, где l₀₁ и l₀₂ — длины стержней в начальном (недеформированном) состоянии. Физические соотношения имеют вид: N₁ = C₁e₁, N₂ = C₂e₂. Эти уравнения в совокупности однозначно определяют некоторое равновесное состояние системы, зависящее от силы P. Они и были подвергнуты нами дальнейшему анализу и решению.

Несмотря на то что работа данного механизма представляется элементарной, математическая формулировка задачи таковой не является. Система уравнений состояния системы является нелинейной, и ее решение в данном случае необходимо строить с применением численных процедур. Для нахождения всего множества равновесных состояний системы, как устойчивых, так и неустойчивых, которое включало бы и все предельные точки¹ , мы воспользовались методом наилучшей параметризации — новым вариантом известного метода продолжения решения по параметру. В книге Шалашилина В.И. и Кузнецова Е.Б. «Метод продолжения по параметру и наилучшая параметризация в прикладной математике и механике» разработан алгоритм перехода к наилучшему параметру продолжения; показано, что такое преобразование насыщает исследование многими новыми полезными качествами, включая как постановку задачи, так и методы построения решения.

Применительно к данной задаче метод наилучшей параметризации сводится к введению нового параметра продолжения l, который определяется условием dl²=dx²+dy²+dP². Это означает, что в качестве параметра продолжения используется длина дуги интегральной кривой в пространстве R³:{P,x,y}. Теперь для получения решения необходимо положить, что все линейные и угловые координаты и нагрузка являются функциями параметра l, и продифференцировать уравнения равновесия по l. В результате решение исходной нелинейной алгебраической системы сводится к интегрированию задачи Коши для системы трех обыкновенных дифференциальных уравнений. Последнее было выполнено численно с использованием абсолютно устойчивого BDF-метода с допустимой локальной погрешностью 10^-12.

Для нахождения всех возможных равновесных состояний системы задача решалась для уровня нагружения 600 Н.

Результаты строгого решения рассматриваемой задачи представлены на рис. 3-7.

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

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

Из рис. 5 и 6 видно, что используемый метод построения интегральной кривой не накладывает ограничений на характер изменения внешней нагрузки: нам известны только начальный (P=0 Н) и искомый уровни нагрузки. Каждой точке интегральной кривой решения задачи соответствует некоторое равновесное состояние системы. Таким образом, вся совокупность возможных равновесных состояний образует бесконечномерное множество, упорядоченное способом введения параметра продолжения. В процессе построения решения как упорядоченной последовательности равновесных состояний задачи нагрузка даже меняет знак (до P=–174 Н), адекватно отображая характер изменения интегральной кривой.

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

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

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

Следующий этап — расчет этой же задачи по методу конечных элементов.

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

Для решения этой задачи мы использовали программный комплекс UAI/NASTRAN 20.0, в котором реализован «классический» вариант метода наилучшей параметризации (Arc Length Method). Здесь также вводится новый параметр продолжения, аналогичный описанному выше, однако решение системы нелинейных алгебраических уравнений не сводится к интегрированию задачи Коши. Аналогичные методы реализованы в большинстве коммерческих программных комплексов (MSC/NASTRAN, MARC, ABAQUS и др.).

Для читателей, не знакомых с NASTRAN, поясним последовательность решения задач с применением этого комплекса на примере схемы, показанной на рис. 9. Вначале (обычно с помощью пре/постпроцессора: FEMAP, MSC/PATRAN и др.) формируется конечно-элементная модель конструкции. Вся информация о модели, равно как и о нагрузках, внешних связях и т.п., записывается в так называемый Исходный файл задачи. Далее NASTRAN выполняет преобразование данных о модели в матричную форму, решает систему уравнений и выводит искомые результаты (напряжения, перемещения, деформации и т.д.) в Файл результатов. Далее этот файл может анализироваться опять же через пре/постпроцессор или иными средствами (например, результаты в табличной форме могут быть загружены в Excel).

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

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

1. Выбор метода решения (метод длины дуги или Arc Length Method).
2. Метод итерации (выбирается автоматически: в зависимости от скорости сходимости выбирается либо метод касательных к кривой s-e, либо метод итерации неуравновешенных нагрузок).
3. Уровень приращения нагрузки (выбирается автоматически на каждой итерации, в зависимости от сходимости решения).
4. Критерии сходимости (по ошибке в уровнях энергии, нагрузки и перемещений).
5. Исходное число приращений нагрузки (при автоматическом уровне приращения нагрузки рассматривается как число, на которое делится полная нагрузка для определения первого приращения).
6. Максимально допустимый поворот в узле модели на каждом приращении.

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

В данном случае наибольшее влияние на ход решения оказали:

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

Счет задачи на UAI/NASTRAN 20.0 занял 13 секунд на первом этапе и столько же — на втором. Решение велось на ПК Pentium MMX-200, ОЗУ 64 Mбайт.

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

SUBCASE 1
 STEP 1
   NLTYPE = GEOM
   LOAD = 1
   NLSOLVE = 4
 STEP 2
   NLTYPE = GEOM
   LOAD = 2
BEGIN BULK
NLSOLVE,4,ARC,AUTO,AUTO,EPU,,71.,,,,,15,,,,,,,,,,15.0
NLSOLVE,5,ARC,AUTO,AUTO,EPU,,71.,,,,,25,,,,,,,,,,15.0

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

Несмотря на указанные проблемы, результаты, полученные на UAI/NASTRAN, очень близки к строгому решению (рис. 10).

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

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

«САПР и графика» 4'2000