Логотип сайта поддержки пользователей САПРО сайте поддержки пользователей САПР Translate to:
Hosting Ukraine

РУССКАЯ РУЛЕТКА или анализ основных ошибок теории геометрически нелинейного метода конечных элементов

канд. техн. наук Назаров Д.И., КузГТУ

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


Рис. 1. Задача №1, конструкция из двух элементов

На рис. 1-4 сплошная линия - исходная схема, пунктирная - правильная деформированная схема конструкции.


Рис. 2. Задача №2, конструкция из трех элементов

Рис. 3. Задача №3, конструкция из четырех элементов

Рис. 4. Задача №4, конструкция из пяти элементов

В первой задаче правый элемент конструкции должен быть повернутым право, а левый повернутым вниз. До вертикального смещения Dy<6,964 правый элемент находится в левом положении, после (т.е. Dy>6,964) правый элемент переходит в правое положение.
Во второй задаче узел с неподвижной опорой после деформирования конструкции должен оказаться правее узла с силой. Для проверки в FEA достаточно поменять местами горизонтальные проекции элементов и помнить, что Rx=P (т.е. реакция правой опоры равна приложенной силе).
Правильный результат третьей задачи так же очевиден, верхние два элемента вытянутся в линию, а узел приложения силы получит горизонтальное смещение вправо. Для проверки достаточно убедится, что верхние два элемента обязательно окажутся растянутыми, и, соответственно, потянут узел вправо.
Четвертая задача предназначена для специалистов, заявляющих, что автор говорит о некой закритической работе конструкций, проявляемой в «прощелкивании» элементов. Уникальность задачи на пять элементов именно в отсутствии «прощелкивания»! Правильный результат - узел приложения силы не перескакивает ниже растягиваемого элемента, а нижние два элемента вытягиваются в линию.

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

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

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

Итак, мы уже убедились, что существующие FEA программы с этими простейшими задачами не справляются. Но в чем проблема? Может быть это лишь случайности? И если взять десять тысяч элементов (а лучше миллион), то, к примеру, ANSYS или MSC/NASTRAN получат правильный результат? К сожалению нет! Но почему? Ответ прост, все FEA программы решают только квазилинейные задачи, а не нелинейные. Дополнительно отметим, что физические квазилинейные задачи описываются, с точки зрения математики, гладкими функциями.

Выше сказанное проще прокомментировать на примере графика функции f(x)=ln(x) изображенного на Рис. 5.


а)б)
Рис. 5. График логарифмической функции: а) с десятичной шкалой осей; б) с логарифмической шкалой абцисс и десятичной шкалой ординат

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

Но с чем связанна такая особенность существующей теории геометрически нелинейного конечно-элементного анализа? А с тем, что существующие методы поиска решений в МКЭ никогда не предназначались для решения нелинейных задач! Продавцы FEA программ назовут подобное заявление абсурдным, компании разработчики FEA программ будут отмалчиваться или заявлять, что это бред, но мы ученые, соответственно совместно постараемся провести беспристрастный анализ теории МКЭ, а для этого рассмотрим систему разрешающих уравнений:

|K|*{u}={F},
где |K| - матрица жесткости конструкции, {u} - вектор перемещений узлов конструкции, {F} - вектор узловых сил. Обратим внимание на особенности, важные в нелинейной постановке задачи:

  • матрица жесткости |K| обязана соответствовать деформированному равновесному состоянию конструкции;
  • вектор {u} является нелинейной функцией от вектора {F}.

Наиболее мощным (на сегодняшний момент) методом для решения подобных задач считается метод Ньютона-Рафсона (full Newton-Raphson procedure). Суть метода описывается итерационными вычислениями уравнения:

|Ki|*{Dui}={Fa}-{Fri}........(1),

где |Ki| - матрица жесткости системы, {Dui} - вектор приращений перемещений узлов, {Fa} - вектор приложенных нагрузок, {Fri} - вектор невязки узловых сил, соответствующий вектору перемещений узлов конструкции {ui}, индекс i -соответственно номер итерации.

Последовательность вычислений следующая:

  1. вычисляется матрица жесткости |Ki| и вектор невязки узловых сил {Fri}, соответствующие вектору перемещений {ui} (для начальной итерации {ui}={0}, {Fri}={0});
  2. вычисляется (см. форм. 1) вектор приращений перемещений узлов {Dui};
  3. определяется вектор перемещений {ui+1}={ui}+{Dui};
  4. повторяются предыдущие шаги итерации до сходимости уравнения 1 (условием сходимости является или |{Fa}-{Fri}|<eF, или |{Dui}|<eu, где e - точность вычислений).

Рис. 6. Графическое представление метода Ньютона-Рафсона

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

Данный метод сегодня применяется для решения большинства геометрически нелинейных задач и признан наиболее мощным и правильным методом. Однако в решении истинно нелинейных задач (например, см. Рис. 1-4) этот метод показывает свою несостоятельность (впрочем как и arc-length метод).

Отметим основные недочеты (если их можно так назвать) в методе конечных элементов, некритичные для квазилинейных задач, но критичные для нелинейных задач:

  1. в корректной нелинейной постановке уравн. 1 не является равенством;
  2. матрица жесткости |Ki| не может быть нулевой или отрицательной;
  3. при выполнении итерации решения уравн. 1 не учитывается взаимное влияние компонентов вектора перемещений {u}.
Рассмотрим более детально описанные недочеты. Для первого недочета достаточно более внимательно рассмотреть график реализации метода Ньютона-Рафсона (см. Рис. 6.). На графике можно увидеть, что на i-ой итерации при определении {Fri} и {Dui} используется |Ki| соответствующий деформированному состоянию предыдущей итерации. Графически (пока) вроде бы верно, но вот физически есть некоторые сомнения, уточнить которые помогут уравнения, формирующие уравнение 1:

|A|*{N}={F} - уравнение статического равновесия конструкции,

{D}=|A|t*{u} - уравнение совместности деформаций,

{N}=|G|*{D} - физические уравнения (соответствие деформаций элементов усилиям),

где {N} - вектор усилий в элементах, {D} - вектор деформаций (удлинений) элементов, |G| - матрица жесткостей элементов конструкции, |A| - матрица уравнений статики, |A|t - матрица коэффициентов уравнений совместности деформаций (транспонированная матрица уравнений статики), которые после преобразования дают

|A|*|G|*|A|t*{u}={F},

и после обозначения |K|=|A|*|G|*|A|t, окончательно получаем |K|*{u}={F}.

Преобразуем формирующие уравнения в итерационные, соответствующие уравнению 1:

|Ai|*{Ni}={Fa}-{Fri} - уравнение статического равновесия конструкции,

{Di}=|Ai|t*{Dui} - уравнение совместности деформаций,

{Ni}=|G|*{Di} - физические уравнения,

учитывая, что в методе Ньютона-Рафсона {Fri} и {Dui} (и соответственно {Di} и {Ni}) относится к текущему состоянию конструкции, а |Ki| (и соответственно |Ai| и |Ai|t) относится к предыдущему состоянию конструкции, переписываем формирующие уравнения в соответствии с физическим смыслом:

|Ai|*{Ni}{Fa}-{Fri},

{Di}»|Ai|t*{Dui} (при значительных относительных перемещениях - неравенство, т.е. {Di}|Ai|t*{Dui}),

{Ni}=|G|*{Di} (закон Гука не зависит от особенностей геометрической постановки задачи),

Прокомментируем уравнения:

|Ai|*{Ni}{Fa}-{Fri} - усилия в элементах текущего состояния, спроецированные в предыдущее состояние конструкции не равны разности приложенной нагрузки и невязки узловых сил в текущем состоянии, несколько странновато звучит по простой причине - тяжело красиво сформулировать абсурдные зависимости;

{Di}|Ai|t*{Dui} - фактически не требует комментариев, т.к. в точной (нелинейной) постановке удлинение элемента не равно сумме проекций на его ось перемещений узлов, элементарная геометрия господа.

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

|Ki|*{Dui}{Fa}-{Fri}........(2)
Разумеется подобное уравнение будет справидливым во всех случаях!!!

Хотелось бы написать «что и требовалось доказать!», но, учитывая, что миллионы пользователей FEA программ и тысячи ученых в области МКЭ «успешно» применяют и развивают FEA, проведем дальнейший анализ.
А может быть, несмотря на отсутствие физического смысла в уравнении 1, применяемый математический аппарат на столько гибок и универсален, что все же позволяет получить правильный результат?
Однако в применяемом математическом аппарате для численного решения системы нелинейных уравнений, а именно в методе Ньютона для векторного уравнения (в отечественной литературе более известного как метод касательных) есть свои проблемы, обычно не обсуждаемые в математической физике и механике.
С математической точки зрения это не проблемы, это всего лишь ограничения области применимости этого метода. А именно - требование отсутствия перегибов функции, т.е. df(x)/dx0 и требование непрерывности и возможности двойного дифференцирования функции на всем рассматриваемом отрезке (т.е. отсутствие вертикальных асимптот, горизонтальных и вертикальных линий, точек излома как для f(x), так и для df(x)/dx).
Первое требование свойственно только методу Ньютона, второе - всем численным методам простых итераций. Однако почти все геометрически нелинейные задачи имеют именно такие особенности, достаточно посмотреть на графики представленные в документации ANSYS (см. Рис. 7.1. главы №7 «Buckling Analysis», Рис. 8.1. главы №8 «Nonlinear Structural Analysis» раздела «Analysis Guides»).

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

Рассмотрим второй указанный недочет существующей теории нелинейного МКЭ - «матрица жесткости |Ki| не может быть нулевой или отрицательной».
На первый взгляд, это высказывание можно и не комментировать, т.к. оно кажется абсолютно справедливым, более того, это высказывание может показаться, с точки зрения уместности, абсурдным, но все же…


Рис. 7. График геометрически нелинейного поведения конструкции

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

Классически поведение конструкции описывается следующим образом: процесс нагружения - кривая ABCD, процесс разгружения - кривая DEFA. Участки кривой BC и EF - скачки конструкции из одного «стабильного» состояния в другое, т.е. участки с отсутствующим промежуточным положением конструкции. Обычно в этих теориях не акцентируется внимание на особенность - конструкция не оказывается на участке BE, хотя он явно представлен на графике F(u), а вот прямые BC и EF введены дополнительно.

Обратимся снова к задаче представленной на Рис. 1. Проведя несложный анализ приходим к выводу, что для определения правильного положения правого стержня нам необходимо оказаться на участке BE (см. Рис. 7). В общей постановке это можно сформулировать следующим образом - «Переход любого элемента конструкции из одного стабильного состояния в другое может привести к переходу любого количества других элементов конструкции из одного стабильного состояния в другое (эффект домино)».
Однако на участке BE (см. Рис. 7) K (жесткость) оказывается отрицательной! Напомним, что K=dF/du, или в графической интерпретации - тангенс угла наклона функции F(u). Очевидно, что с подобным препятствием ни одна программа не справится, т.к. это противоречит существующей теории МКЭ и механике в целом.

Хотелось бы верить, что существующие FEA программы (да и теория МКЭ) каким-то чудом обходят это препятствие. Однако, формулируя четко, необходимо сказать, что и из-за этого «недочета» все геометрически нелинейные FEA программы не способны получить правильного результата!

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

Проанализируем последовательность определения промежуточных положений конструкции с использованием уравнения |Ki|*{Dui}={Fa}-{Fri}.
В отличие от предыдущего «недочета» в этот раз мы разберем не конкретную итерацию метода, а взаимосвязь последовательных итераций.


Рис. 8. Исходная модель к анализу взаимосвязи итераций в нелинейном МКЭ

За основу примем конструкцию (см. Рис. 8) из одного элемента (x и y - горизонтальная и вертикальная проекции) с одной опорой и одной силой (Px и Py - горизонтальная и вертикальная составляющие), разумно предполагая, что нелинейное поведение нескольких элементов окажется более сложным, чем поведение одного.

Рассмотрим систему разрешающих уравнений:

G*sin2(a)*Dy=Py-DPy
G*cos2(a)*Dx=Px-DPx,

где: G=EF/l - жесткость элемента; cos(a)=(x-Dx)/((x-Dx)2+(y-Dy)2)1/2 - косинус, а sin(a)=(y-Dy)/((x-Dx)2+(y-Dy)2)1/2 - синус угла наклона элемента; Dx - горизонтальное, а Dy - вертикальное перемещение узла приложения силы; DPx - горизонтальная, а DPy - вертикальная реакция элемента.

Построим график зависимости горизонтальной проекции силы (Px-DPx=Px(x+Dx)), в зависимости от горизонтального смещения узла (Dx), при исходных данных: EF=1, l=(5+5)1/2; x=5; y=5; Dy1=-5, Dy2=0, Dy3=5.


Рис. 9. График Px(x+Dx)

По графику (см. Рис. 9) несложно убедится, что при решении системы уравнений независимое определение Dx и Dy приводит к определению нового положения элемента абсолютно не соответствующего предыдущему положению. Фактически определяются два независимых новых положения, соответствующих точкам [(x+Dx);y] и [x;(y+Dy)], но, при определении DPy и DPx необоснованно используется третья точка [(x+Dx);(y+Dy)].
С физической точки зрения - это перемещение из текущего (неравновесного) состояния в новое (случайное), которое может быть как более, так и менее неуравновешенным.
Рассматривая график (см. Рис. 9), убеждаемся, что подобное определение неправильного нового положения, для достаточно большого числового диапазона конкретных значений, приводит к определению возможно принципиально неправильного DPx, например, вместо отрицательного - положительное (и наоборот). Следствием подобной ошибки является возможное увеличение невязки (Px-DPx), по сравнению с предыдущей итерацией.
Соответственно, вместо сходящегося решения мы получаем расходящееся! Отметим, что с использованием аппроксимации следующего положения по дуге интегральной кривой (arc-leght method) новая точка может оказаться вообще где угодно, например [(x-Dx/2);(y-Dy/3)].

Попробуем рассмотреть эту проблему с математической точки зрения, вполне возможно, применяемый математический аппарат имеет скрытые резервы, позволяющие обойти физические ошибки.
Построим графики зависимости горизонтальной (см. Рис. 10) и вертикальной (см. Рис. 11) проекций реакций элементов в зависимости от проекций элемента, в нашем случае при EF=1, l=(5+5)1/2; x=5; y=5.


Рис. 10. График Px(x,y)

Рис. 11. График Py(x,y)

С математической точки зрения, определение правильного деформированного состояния конструкции - это пошаговый переход из точки, соответствующей исходным координатам x и y (лежащей на нулевом круге), через промежуточные точки, в точку, соответствующей значению силы P (с координатам x’ и y’). При этом метод определения положения следующей точки (направление и величину шага) напоминает «метод покоординатного спуска», давно известного в математике.
В математике этот метод используется для определения точки минимума функции. В нашем случае, проекции сил Px и Py имеют различное положение своих минимумов и соответственно будут «притягивать» решение в различных направлениях.
Более того, в моделируемой физической задаче точка, соответствующая силе P, не является минимумом функции, что противоречит математической постановке задачи.
Другой момент, для задачи даже с одним элементом, соответствующих правильному положению конструкции при заданной силе P (с проекциями Px и Py), правильных решений может быть два (см. Рис. 10, 11).

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

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

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

Стоит ли подвергать опасности жизни тысяч людей?

Такое могут позволить себе только «любители» FEA, а профессионалы обязаны не совершать ошибок!

В последнее десятилетие резко увеличилось количество конструкций кораблей, самолетов, ракет, лодок и т.д., полностью запроектированных с использованием современных средств САПР, но в неполной мере экспериментально проверенных (наиболее легендарный - Боинг 747).

Резко увеличилось и колличество катастроф...

Господа - будьте профессионалами, не играйте в русскую рулетку с чужими жизнями!



Copyright © Сайт поддержки пользователей САПР