CC BY
108
УДК 621.031:531.31 + 539.37
О КОРРЕКТНОСТИ РЕШЕНИЯ ГЕОМЕТРИЧЕСКИ НЕЛИНЕЙНЫХ ЗАДАЧ ПАКЕТОМ МКЭ
О. С. Садаков, А.А. Шапиро
Концепция пренебрежимой малости смещений, поворотов и деформаций, которая позволяет ограничиться геометрически линейной постановкой задачи, значительно облегчает расчеты в механике деформируемого тела, но в ряде актуальных задач ведет к серьезным ошибкам. Так, учет геометрической нелинейности необходим для описания деформирования конструкций из высокоэластичных материалов, конструкций из тонких оболочек, для расчета процессов, связанных с большими пластическими деформациями.
Во многих интенсивно развивающихся в последнее время пакетах программ, основанных на методе конечных элементов, заявлена возможность решения геометрически нелинейных задач. Однако известно, что далеко не всегда решение оказывается корректным. В силу значительной сложности проблемы теорию конечных деформаций нельзя назвать законченной; используемые в пакетах теоретические работы, посвященные учету геометрической нелинейности, как правило, корректны лишь в частных случаях, хотя в документации часто декларируются их более широкие возможности.
Ниже предлагается простая тестовая задача, позволяющая оценить корректность соотношений, используемых тем или иным пакетом для геометрически нелинейной задачи. С ее помощью расчетчик может проверить используемый им программный продукт. Приведен иллюстративный пример с пакетом А^УБ.
1. Рассмотрим состояние однородной деформации (дисторсии). При этом любое семейство прямых равноотстоящих параллельных волокон переходит в новое положение, оставаясь семейством прямых равноотстоящих параллельных волокон. Это преобразование линейно и характеризуется тензором дисторсии Б:
(1)
Здесь вектор / определяет текущее положение произвольного волокна, /0 - его начальное положение, точкой здесь и далее обозначено скалярное произведение.
Иногда удобнее работать с другим тензором, также характеризующим дисторсию:
Я = (2)
где I - тензор тождественного преобразования (единичный). Тензор В определяет линейную зависимость «изменения» волокна от /0:
Д/з/-/0 =£>-/0. (3)
Ограничимся для простоты рассмотрением плоского деформированного состояния. Тензор однородной плоской дисторсии однозначно определяется относительным смещением трех точек. Поэтому для анализа кинетики деформированного состояния удобно рассмотреть треугольный
Закрепив узел С, будем задавать перемещения узлов А и В и тем самым деформировать треугольник, сохраняя следующее свойство: обход вершин А, 5, С должен в любой момент времени I осуществляться против часовой стрелки. В конечный момент деформирования возвратим треугольник в исходное положение; конечная деформация будет нулевой.
2. Для определения текущей дисторсии достаточно, кроме начального положения двух волокон (сторон треугольника а0, Ь0\
изопараметрический конечный элемент (рис. 1).
Рис. 1. Эволюция треугольного элемента
знать их текущее положение (a,b). Тензор D линейно связан с изменениями Aa = a~a0 и
Ab = b~b0. Если АЬ = О, то тензор D представляет собой диаду (обозначим ее cd), при этом
D • b0 = cd • b0 ~ 0. (4)
Введем векторы a = Э-а0, Ь' = Э • ¿>0, где Э - тензор поворота на угол л/2; в декартовом базисе с
ортами в\у е2 он равен e2e¡ - Из условия (4) вектор d можно принять равным V . Второй вектор диады с найдется из условия
D • д0 = cd • а0 = Аа . (5)
Получим D - Aab' Ib' • а0 . Из упомянутой линейности найдем, что в общем случае, когда АЬ * О,
А Ъа-АаЬ'
(6)
D =
а0-Э-Ь0
(здесь учтено, что х- Э-у = -у - Э х).
Зная тензор F -1 + Z), нетрудно решить, например, следующую задачу факторизации:
F = RU, (7)
где ¿У - симметричный, а Л - ортогональный тензоры.
В соответствии с выражением (7) предполагается (вполне произвольно), что переход конечного элемента в текущее состояние производился вначале деформированием (тензор U), а затем жестким поворотом, характеризуемым тензором R. Тензор U называют тензором коэффициентов длины, он равен 1+£, где обычный тензор деформации, определяющий Д/ = £-/0 при отсутствии поворота как жесткого целого, то есть при неподвижных главных волокнах. Зная тензор U, нетрудно найти тензор логарифмической деформации (тензор Генки) е:
e = lnU. (8)
Таким образом, задав кинематическое нагружение конечного элемента с помощью двух замкнутых траекторий А А (/) и А В (/), в любой момент истории нагружения и в любом порядке (не обязательно хронологическом) можно найти тензоры R, U, е , е. В рассматриваемом случае плоского деформирования по R нетрудно найти и угол поворота элемента как жесткого целого. Для этого достаточно знать две координаты: =COS <р И /?21 = sin (р .
3. Для примера было произведено тестирование пакета ANSYS 6.1 (разработчик - ANSYS Inc). Данный программный продукт является бесспорным лидером, как по количеству зарегистрированных пользователей, так и по заявленным возможностям. Открытая архитектура ANSYS и достаточно полная документация [1], содержащая теоретическое описание используемых методов [2, 3] делает этот пакет удобным объектом для проведения такого теста.
Процедура расчета кинетики деформирования треугольника Triangle представлена на листинге 1. Описание алгоритма расчета выполнено при помощи так называемого алгоритмического
«языка» - псевдокода, который легко переводится на любой язык программирования. Здесь использована некоторая модификация псевдокода, предложенного Дэвидом Роджерсом [4]. Описанная процедура была реализована авторами статьи в пакетах компьютерной алгебры Maple и MathCAD.
Была рассмотрена замкнутая траектория деформирования треугольника, имеющая три этапа (рис. 2): 0-1, 1-2 и 2-0. Каждый этап разбит на 10...100 шагов, и тензоры деформации вычисляли на каждом шаге. Их сравнивали с данными, полученными из аналогичного расчета в пакете ANSYS (листинг 2). На рис. 3, 4 показаны результаты сравнения ве-
Рис. 2. История деформирования конечного максимальной логарифмической де-
элемента. Конфигурации в моменты времени У ^ ^
о и з совпадают формации и разности е}-е2, в какой-то
мере характеризующей интенсивность сдвига. Видно, что лишь в течение четверти первого этапа результаты совпадают, далее А^УБ сильно завышает деформации и при возвращении элемента в начальное положение (дисторсия становится нулевой) накопленная ошибка не исчезает.
процедуры Triangle и пакета ANSYS процедуры Triangle и пакета ANSYS
В заключение сделаем два замечания. Во-первых, предложенный пример довольно жесткий: элемент испытывает одновременно большие сдвиговые деформации и поворот на 360 градусов. Известно, что такой случай рассчитывается имеющимися на сегодняшний день методами с наибольшей погрешностью. Во-вторых, авторы статьи не готовы предложить корректный алгоритм для решения краевых задач в терминах конечных деформаций и не хотели бы заниматься неконструктивной критикой используемых методов расчета. Надеемся, что предложенная схема может оказаться полезной расчетчикам для оценки достоверности получаемых ими результатов. По-видимому, более детальный анализ по предложенной схеме может помочь определить примерные границы области корректного применения тех или иных пакетов.
Листинг 1
Псевдокод процедуры расчета кинетики деформирования Triangle
N - число моментов времени, в которые производится расчет текущего состояния. {а}ь {b}k - векторы, задающие положение узлов А и В на плоскости в каждый момент времени к. Координаты векторов определяются в предположении, что декартовый базис имеет начало координат в неподвижном узле С.
Формирование векторов {а}к и {b}k, описывающих конфигурацию треугольника во времени, легче запрограммировать, имея несколько (1..10) характерных положений его узлов и задавшись числом подшагов, за которые происходит переход от одной конфигурации к другой. Узел переходит в новое состояние из предыдущего по прямой линии. Таким образом, движение узла характеризуется кусочно-линейной траекторией.
Дисторсия({а}к, {b}k, {а}0, {Ь}0) - функция вычисляет тензор дисторсии по двум волокнам в исходном ({а}0, {Ь}о) и деформированном ({а}ь {Ь}к) состояниях КореньТ([Х]) - операция взятия квадратного корня от тензора [X] ЛогарифмТ([Х]) - операция взятия натурального логарифма от тензора [X]
Поворот({х}) - функция поворачивает вектор на 90 градусов в 2Б-пространстве против часовой стрелки
Eigenvals([X]) - функция возвращает вектор собственных значений тензора [X] Lsolve([A], {b}) - функция возвращает решение {х} СЛАУ вида [А]-{х}={Ь}.
begin
for к - 1 to N
[F]k = Дисторсия({а}к, {b}k, {а}0, {b}0) [и]к-КореньТ( [F]V[F]k) [e]k = [U]k - [I] [e]k = ЛогарифмТ( [U]k) [R^tFMUK1 next к
Визуализация результатов
function Дисторсия ({a}, {b}, {aO}, {bO})
Функция вычисляет тензор дисторсии по двум волокнам в исходном и деформированном состояниях
{а'} = Поворот({аО}) {Ь1} = Поворот({Ь0})
[F] = ({а} - {аО})-{Ь'}Т/ ({Ь'}-{аО}) + ({Ь} - {Ь0})-{ а'}Т/ ({а'}-{Ь0}) + [I] return [F]
х\
function Поворот ({х})
0 -1
1 О return {у}
function КореньТ ([X])
Операция взятия квадратного корня от тензора = Eigenvals([X])
if( А^оЛу ) then
СС(
а,
Lsolve
1 Л) 1 л
Л
[Y] = а0-[1] + а г [X]
else
cy] = Vvm
end if return [Y]
function ЛогарифмТ ([X])
Onepaifm взятия натурального логарифма от тензора П\ = Eigenvals([X])
if ( Aq о Л1 ) then
I - Lsolve
"1 Л, _1 Я,
1п(Л).
[Ц + аг[Х]
else
[Y] = ln(X,o)-[I]
end if return [Y]
* Замечания о нотации:
• Векторные и тензорные величины отображаются координатами в выбранном декартовом базисе.
• Квадратными скобками обозначены квадратные матрицы (2x2). Фигурным скобкам соответствуют векторы-столбцы (2x1).
Листинг 2
Командный файл расчета для пакета ANSYS
/ргер7
len - 2*sqrt(3)
АхО =2 $ АуО = len
ВхО = -2 $ ВуО - len
lsnum = 3 ssnum = 100
*dim,action,array,4, lsnum J Ax Ay Bx By
action ( 1,1 action(1,2 action (1,3. !__________
2, -2, 9, 2 -7, 1, -4, -10 2, len, -2, len
et,1,PLANE2
keyopt,1,3,2 mptemp, mptemp,1,0 mpdata,ex,1,,2e5 mpdata,prxy,1,,0.3 n,1,0,0, n, 2, AxO, AyO, $ e,1,2,3
d,l,UX,0, $ /solu
antype,0 nropt,,,off nlgeom,on
*do,i,1,lsnum t ime, i
*get,ax,node,2,loc,x *get,ay,node,2,loc,y *get,bx,node,3,loc,x *get,by,node,3,loc,y dcum,repl,
d,2,UX,action(1,i)-ax d, 2,UY,action(2, i)-ay d,3,UX,action(3,i)-bx d,3,UY,action(4,i)-by nsubst,ssnum,,ssnum, off
n,3,BxO,ByO, d,1,UY,0,
I включили препроцессор
! Начальное положение узлов А и В;
! Число шагов нагружения;
! Число подшагов/
I История натружения:
! координаты узлов в каждый
! момент времени;
тестировано и для элементов УХвСОЮб; плоское деформированное состояние; свойства материала (не имеет значения)
Е - 2е5 МПа V - 0.3
элемент по 3 узлам начальные граничные условия включили решатель
тип задачи - статический прочностной расчет
расчет с учетом геометрической нелинейности
цикл по моментам времени от 1 до 1зпит время 1
узнаем текущее положение узлов
задаем новые граничные условия
число подшагов = ssnum
outres, all, all, solve *enddo
/post26
xvar, 1
esol, 2, 1, 2, EPEL,x,ex esol,3, 1,2,EPEL,y,ey esol,4 , 1,2,EPEL,xy,exy /axlab, x, time /axlab, y,strain /GRID,1 p]var, 2,3,4 finish
сохранять всю информацию на каждом подшаге решаем задачу конец цикла
включить time-history постпроцессор по оси абсцисс время
переменная 2 - упругая деформация ех переменная 3 - упругая деформация еу переменная 4 - упругая деформация еху надписи по осям
сетка на графике
отобразить переменные 2,3,4 на графике
Литература
1. ANSYS, Inc. Theory Reference, Chapter 3. Structures with Geometric Nonlinearities.
2. Huges T.J.R. Numerical implementation of constitutive models: rate-independent deviatoric plasticity. // Theoretical foundation for large-scale computations for nonlinear material behavior. -Martinus Nijhoff Publishers, Dordrecht, Netherlands. - 1984.
3. Weber G.G., Lush A.M., Zavaliangos A. An objective time-integration procedure for isotropic rate-independent elastic-plastic constitutive equations // International journal of plasticity. - Vol. 6, 1990. - P. 701-749.
4. Роджерс Д. Алгоритмические основы машинной графики. - М.: Мир, 1989. - 512 с.
Поступила в редакцию 27 апреля 2003 г.
