Конечно-элементная реализация метода геометрического погружения на основе вариационного принципа Кастильяно для плоской задачи теории упругости

Метод конечных элементов (МКЭ) представляет собой широко известный численный метод решения задач механики деформируемого твердого тела. Наибольшее распространение получил вариант МКЭ для задач теории упругости в перемещениях, формулировка которого может быть получена на основе вариационного принципа Лагранжа (принципа минимума общей потенциальной энергии системы). Метод конечных элементов в напряжениях, основанный на вариационном принципе Кастильяно, также получил свое развитие, но в связи с рядом проблем имеет ограниченное применение [1]. В частности, имеются сложности в построении элементов с криволинейными границами при решении задач для областей произвольной, неканонической конфигурации. В данном случае может оказаться полезным использование идеологии метода геометрического погружения [2], позволяющего свести краевую задачу для тела сложной формы к итерационной последо- вательности задач для тела канонической конфигурации, при этом возможно применение конечных элементов в напряжениях простой формы [3].

1. Рассмотрим вариационную постановку задачи теории упругости в рамках вариационного принципа Кастильяно. Функционал дополнительной энергии упругой системы имеет вид

^П C ^ ) = ¹ К' ^£ ( ^У ) ^d V "  I t i^ ) U i dS u ’              ⁽¹⁾

2 у              S

u где V - область, в которой разыскивается решение; су , £ - тензоры напряжений и деформаций с компонентами oj (F), Ej (г) соответственно; F=(x,y,z)gV; U - заданные перемещения на поверхности Su; t - усилия на поверхности, где заданы кинематические граничные условия.

Функционал (1) определен на множестве статически допустимых полей напряжений yj (F), удовлетворяющих в области V уравнениям равновесия и статическим граничным условиям на части поверхности тела S (S uS = S). 0X0 и y jj a j = T ’ F G Sy, где T - заданные усилия на поверхности Sy; aj - компоненты вектора внешней единичной нормали к поверхности тела.

Условие минимума функционала (1)

5П C = 0

эквивалентно выполнению в области V уравнений совместности в напряжениях (уравнений Бельтрами-Митчела) и граничных условий в перемещениях на границе S u

U = Ui , ^{F G S} u •

Введем в рассмотрение некоторую каноническую область V₀, включающую исходную область V произвольной формы, так что V₀ = V U V , где V \ - дополнение области V до области V₀.

Выполним тождественное преобразование функционала (1) дополнительной энергии линейно-упругого тела:

^п с (^а ) = 1 J^a j г _ч (^а) d V + ¹ J ^ i,. г, ( а ) d V . ^{- 1} j ^ г, (^а) d V . ^- j t_i ( ^a U d S u ,

2 V                2 V.                2 V.

П C (a ) = 1 Jaij 8 j ^) d V -1 Jaj" 8j (a )d V.- j ti(a) UidSu •

• 2 Vo                2 V.

Можно показать, что поля напряжений, минимизирующие функционал (2), в области V доставляют минимальное значение функционалу (1).

Условие минимума функционала (2) примет вид

8П с (а) = jsaj) d V - js^a) d V.- jti (8a)Ui dSu = 0.(3)

V                   V.

Вариационное уравнение (3) будем решать, используя следующую итерационную процедуру:

j8ai,

aj⁰⁾ = 0, k = 1,2,3...,n, где k - номер итерации.

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

• 2.    Рассмотрим процедуру построения дискретного аналога вариационного уравнения (4) с помощью метода конечных элементов в напряжениях для плоской задачи теории упругости в декартовой системе координат. Для дискретизации выберем прямоугольный конечный элемент с четырьмя узлами в углах (рис. 1).

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

NT—к

^Лдф^Л

^дф) ^

^2 Л д²ф

I dx Д dу )\дx ду ,

. Подробно использование данного эле

мента описано в статье [4].

Рис. 1. Прямоугольный конечный элемент в локальной системе координат

Уделим особое внимание статическим граничным условиям, являющимся главными для данной вариационной формулировки. Для выполнения главных граничных условий можно использовать традиционный прием, заключающийся в модификации глобальной матрицы податливости, задавая значения функции напряжений и ее первых производных на границах. Такой подход вызывает затруднения, связанные с определением узловых неизвестных на границах, поэтому для удовлетворения статических граничных условий широко используется метод множителей Лагранжа [5, 6]. Применение данного метода иллюстрирует следующий пример: пусть на стороне элемента, параллельной оси Ox, действуют поверхностные силы Ty , как изображено на рис. 2.

Связь функции напряжений Эри и компонент тензора напряжений имеет вид д2ф       д2ф         д2ф ст — ——, ст — ——, т —--.

^xxду^{2    уу}дx^{2    xy}    дxду

y

a

◄---------►

i

->

x

Ty (x)

Рис. 2. Элемент, на сторону которого действуют поверхностные силы

Дважды проинтегрировав по переменной j на стороне элемента с узлами i и j выражение

д²ф 5x²

Ty (x) и определив константы интегриро

вания через ф и

(5Ф)

V dx у

в узловых точках

( ф

V

Гдф^

’ V dx )J

Фj

(SI

V dx / j у

, получим

(^

V dx у

+

(дф^Л

V dx у j

xj

= T (x) dx, xi

-ф,+фj

(^

V dx у

xj x a = jT (x) dx dx.

xixi

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

[ g >e }={ xe}.

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

Согласно процедуре метода конечных элементов в напряжениях выражения, входящие в уравнение (4), можно представить для отдельного конечного элемента в матричном виде:

|§Ч/⁽k- (^k⁺¹⁾)dV = ⁵{ф^е(k⁺¹⁾}[k^eJ₀{ф^e(k⁺¹⁾},

V

J5aj^kj^k )dV .^-'^k+1>}[k^eJ_a{₉^еk>},

Ve

Jti (sd<^k+1V)uds„ =з{ф^^k+1)}{ pu},

Su где [ke J0 - матрица податливости элемента, принадлежащего канонической области; [ke J^ - матрица податливости элемента, принадлежащего области дополнения; {Р^} - вектор обобщенных узловых перемещений элемента.

Проведя теперь типовую процедуру ансамблирования конечных элементов и учет статических граничных условий методом множителей Лагранжа, получим итерационную последовательность СЛАУ:

Г[к1, [gГ>k■ :

,

<

_ [G] [0] J /л'^k*»}_j

=[K>* ^k 'j+EFF Д

S^}J

где [K]0 =^[keJ0 — глобальная матрица податливости канонической e области; [ к к=яke J^ - глобальная матрица податливости дополне-e ния; {Ри}-^{РЩ} - глобальный вектор обобщенных узловых переме-e щений; {ф} - глобальный вектор узловых неизвестных; {X} - вектор множителей Лагранжа.

Условие остановки итерационной процедуры выбрано в виде

{Фk"-;   ;_1P-k»}

< 0,001,

ф где в качестве нормы вектора {ф} выбрана норма Чебышева.

• 3.    Рассмотрим задачу о плосконапряженном состоянии прямоугольной упругой пластины с прямоугольным вырезом при продольном ее растяжении переменной распределенной нагрузкой и заданными на контуре выреза компонентами вектора перемещений (рис. 3). Так как задача симметрична относительно осей Ox и Оу, будем рас

  сматривать четверть пластины.

  

  Ъ

  Рис. 3. Расчетная схема задачи

  
  
  
  
  
  

Граничные условия задачи: Ux (0, у ) = 0, Uy (x,0) = 0, у е[Ъ/2,Ъ], xе [a/2,a],

U_x (x, b/2) = 0,1E -10 • x,   U_y (x, bl2) = 0,1E -10,   x е [0, a/2],

U_x (a/2,у) = 0,4E -10,   U_y (a/2,у) = 0,05E -10• у,   у е [0,Ъ/2],

^xx ( a, у ) = T ( у ) , ^Тxy ( a, у ) = ⁰у^е[⁰, ^Ъ] ,

^Тxy (^x, ^Ъ) = °уу (^x, ^Ъ) = ^{0 x е}[^0,a] ,

Т xy ( 0 у ) = 0 у е[ Ъ/ 2, Ъ ] , тxy (x,0) = 0, x е[a/2,a].

Для конкретизации зададим a = 8 м, Ъ = 4 м, E = 2Е + 11 Па - модуль Юнга, v = 0,35 - коэффициент Пуассона.

Согласно введенному выше условию V0 = V U VA каноническая область будет представлять собой прямоугольник с координатами вершин {(0,0); (a,0); (a, b); (0, b)} . Следует отметить, что для канонической области граничные условия при x = 0 и у = 0 будут справедливы для всей стороны.

Построим итерационную процедуру для поставленной задачи в форме (4):

ba

J J(5c xx(k ' ;: xx (а< k') + 5а„(k 'г. „ (а(k+1))+5т „ < k ; „ (а< k')) dxdу = b/2 a/2

= J J (^5аxx^(k ^ xx (^{С( k}      уу^{( k}^ уу (^С(k ^^ xy^{( k+1)}Y xy (^{а( k )})) dX^dу +

+ J (tx (5а<^k+1)) Ux (x, b/ 2) +1, (5а<^k' ) U„ (x, bl2 )) dx + 0

b/2

+ J(tx (5а^(k 4Ux (al2,у) +Ц (5а^(k+1))U_y (a/2,у))dу. 0

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

Все вычисления произведены в среде MATLAB 7.0.1. Для выявления оптимальной сетки поставленная задача была решена методом конечных элементов в напряжениях. Исследование метода геометрического погружения, основанного на вариационном принципе Кастильяно, производилось при числе степеней свободы, равном 1284. Итерационная процедура сходится за 20 итераций при остановке по норме (8).

Поставленная задача также решена в программном комплексе ANSYS с использованием формулировки в перемещениях при аналогичном числе узловых неизвестных (рис. 4).

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

Рис. 4. Распределение компонент тензора напряжений, полученных методом конечных элементов (ANSYS) на основе вариационного принципа Лагранжа: а - касательное напряжение тxy ; б - нормальное напряжение ayy ;

в - нормальное напряжение о_хх

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

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

СТyy , ^Па

с xy,^Па

д

Рис. 5. Графики распределения компонент тензора напряжений: а, б, в - при y = b / 2 ; г, д, е - при x = a / 2; * — метод геометрического погружения

е

в напряжениях,        метод конечных элементов в перемещениях

Рис. 6. Сходимость итерационной процедуры с увеличением числа итераций: а - на примере ст_хх, Па на границе у = 0, x е[4,8]; б - на примере стyy, Па на границе у = 0, x е[4,8]; —*— решение, полученное методом конечных элементов в перемещениях; * — решение, полученное методом геометрического погружения на каждой итерации

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