The effect of approximating functions in the construction of the stiffness matrix of the finite element on the convergen.pdf Одним из наиболее эффективных численных методов, позволяющих выполнять анализ напряженно-деформированного состояния конструкций из эластомеров, является метод конечных элементов (МКЭ). Применение МКЭ накладывает определенные требования на системы автоматизированного проектирования (САПР) конструкций из эластомеров. Для расчетов с использованием МКЭ проектируемый объект должен быть представлен с помощью конечного числа элементарных объемов, называемых конечными элементами (КЭ) ["-4]. Задача дискретизации объекта на КЭ является весьма сложной и трудоемкой, особенно в трехмерном случае, что приводит к необходимости введения в САПР модуля, позволяющего автоматизировать данный процесс. МКЭ основан на физической идеализации: одно из двух условий статического или кинематического равновесия выполняется точно, а другое приближенно. Это привело к тому, что МКЭ применяется в двух видах: в виде метода сил и метода перемещений. На практике чаще применяется метод перемещений, так как при его использовании легко представляются кинематические связи для произвольных конструкций. Однако в некоторых случаях более предпочтительным является метод сил, что приводит к необходимости реализации в САПР некоторого механизма выбора требуемого метода расчета ["]. Влияние аппроксимирующих функций при построении матрицы жесткости 27 При применении МКЭ возникает система линейных алгебраических уравнений (СЛАУ), которая в общем случае может иметь весьма высокий порядок. Для ее решения приходится использовать специальные модификации методов решения СЛАУ, как прямые, так и итерационные. После решения СЛАУ необходимо вычислять дополнительные параметры напряженно-деформированного состояния, например по полученным перемещениям вычислять деформации и напряжения. Это приводит к появлению больших массивов численных данных, которые, в первую очередь, приходится исследовать на точность, достоверность и соответствие физическому смыслу задачи [2J. Данный процесс по времени и трудоемкости иногда превосходит все предыдущие этапы проектирования и расчетов. Доказательства сходимости и скорости сходимости схем МКЭ рассмотрены в работах О.К. Зенкевича [3J, М.В. Урева [5J, А.С. Сахарова [6J, А.А. Гусева [7J, В. Nedjar [8J и др. Таким образом, структуру САПР, базирующихся на применении МКЭ, условно делят на три взаимосвязанных блока: 1) подготовка начальных данных для расчетов, включая описание объекта конструирования, дискретизации его на КЭ, описание нагрузок, физических и геометрических параметров, уравнения состояния и краевых условий; 2) расчет с использованием МКЭ; 3) анализ полученных результатов, включая синтез дополнительных параметров напряженно-деформированного состояния на основе ранее вычисленных. Покажем, как работает первый блок, ведь он является важным с точки зрения теоретического обоснования математического обеспечения САПР рассматриваемых объектов. Остановимся на вопросе дискретизации КЭ, а именно на вопросах выбора типа КЭ и функций формы, а также условий, которым они должны удовлетворять. Обычно при построении матрицы жесткости трехмерных шестигранных криволинейных конечных элементов удобным является при произвольном порядке аппроксимации перемещений использование произведений одномерных интерполяционных полиномов Лагранжа [9J. Однако с целью упрощения вычислений интегралов, входящих в выражения для коэффициентов матрицы жесткости (МЖ), иногда целесообразным становится использование степенных координатных функций [7J. При этом построенная матрица коэффициентов для системы степенных функций (либо другой системы координатных функций) с помощью тензоров преобразования приводится к матрице жесткости КЭ. Это справедливо для любого принимаемого закона аппроксимации перемещений [7-9J. Рассмотрим следующие законы. Вариацию энергии деформации конечного элемента представим в известной форме: δW=∫∫∫ σiJδεiJdv . v Компоненты тензоров напряжений и деформаций примем согласно следующим зависимостям: где σiJ =CiJklεkl; εkl = 2 +jkl +ξlk ) , ξkl =Ckn u,ln ; (1) (2) (3) 28 Р.В. Киричевский, А.В. Скринникова Ck" - компоненты тензора преобразования базисной z k и местной Xk систем координат. Штрих в индексах указывает на принадлежность компоненты к базисной декартовой системе координат (рис. 1). Fig. 1. Basic and local FE systems Учитывая выражения (1) - (3), уравнение для вариации энергии деформации перепишем в виде где δW=∫∫∫Hm i" lu,l" δ(u,im )dv, v Hm i" l = C jikl Cm C" . jk Введение нового тензора H m i" l как свертки тензора упругих постоянных с (4) (5) тензором преобразования координатных систем в несколько раз сокращает объем вычислений при построении МЖ КЭ [6]. Влияние аппроксимирующих функций при построении матрицы жесткости 29 Примем произвольный закон аппроксимации перемещений внутри КЭ: M ∑bn' ψ(μ)ψ(ν)ψ(η) (μνη) (1) (2) (3) μ,ν,η=0 с использованием системы координатных функций вида ψμ, = (i. ψV) = (≤r ; ψ,,) = (^3r ψ(1) = μ! . ψ(2 ) = ν! . ψ(3) = un (6) η! (7) μ! MMMM ∑M= ∑M∑M∑M, μ,ν,η=0 μ=0ν=0η=0 M =1 - для линейной аппроксимации. M = 2 - для квадратичной аппроксимации. M = 3 - для кубической аппроксимации и так далее. На основании выражения (6) получим следующее выражение для распределения производной функции перемещений внутри КЭ: un = ∑ b,n' ψ'μ-δ1 )ψ'ν-δ2 )ψ'η-δ3) ρ u,l b(μνη) (1) (2) (3) , μ,ν,η=0 где ρ = [δ1 (μ-1) +1] [δ2 (ν -1) +1] [δ3 (η -1) +1] . Подставляя (9) в (4), получаем следующее выражение вариации энергии деформации: δw=∫∫∫Hmn ∑ b'μνη, ∑ δ(b,mβγ,)ψjμ,+α-δ1 -δ1 >ψ(2'+β-δ2-•*iψ(3+γ-δ3-δ3>× -1-1-1 μ,ν,η=0 α,β,γ =0 где (8) (9) где (μ+α-δ1l-δ1i)!(ν+β-δl2-δi2)!(η+γ-δl3-δi3)! =∑∑∑bi∙μ.,, ∑∑∑δ( bmαβγ) )fiyα*'=ς∑F(j φ)δu''j >u^φ), μνη αβγ f φ f =μ+(M +1)ν+(M +1)2 η+1, φ=α+(M +1)β+(M +1)2 γ+1. (10) (11) Выражение (10) позволяет строить МЖ относительно коэффициентов b(μνη) и δb(αβγ) при произвольном законе аппроксимации с использованием численного интегрирования. В случае принятия допущения о постоянстве тензора упругих постоянных и компонент метрического тензора внутри КЭ становится возможным точное интегрирование в выражении (10). В этом случае будем иметь l(μα)I(Vl)I(ηγ)= J J J иП u,m'dv = li(αβγ) li(α) li(β) li(γ) li -1-1-1 г0, при S ∨ t ∨ g - четные, [23 ω, при S ∧ t ∧ g - нечетные. (12) 30 Р.В. Киричевский, А.В. Скринникова Здесь ω= . . . , stg (μ-δ" )!(α-δ" )!(ν-δl2)!(β-δi2)!(η-δ3 )!(γ-δ3)! s =α+μ-δ"ι-δ".+", t =β+ν-δι2-δ.2+", g =γ+η-δι3-δ.3+", I (μ) =P,/ li(a) [γ s (μ-δ" )!(α-δ")!; J ν) P, li(β) ∖2∣t (ν-δ2 )!(β-δ2)!; I (η) =J0', Z li(Y) [У g (η-δ3 )!(γ-δ3)!. ", 2, ("3) ("4) ("5) ("6) ("7) ("8) ("9) ". " Переходя к пределам интегрирования - ≤ x ≤ - чения интегралов ("7)-("9): 7(μ) =ґ " Y" α !μ! =ґ " Y" ∣^δ"(^^-") + "]У"(а-") + "] , li'α'^ І2 J (μ-δ" )!(α-δ")! I2J s ’ τ(ν) = Ґ"Y-" ν!β! = Ґ"Jt-" [δl2 (v-") + "][δ2 (в-") + "] ; li(e) 12 j (ν-δl2 )!(β-δi2)! 12 J I~(n) =r "Jg-" n!Y! = r" J li(Y) 12 J (η-δ3 )!(γ-δ3)! 12 J g Матрица коэффициентов на этом этапе вычисляется из зависимости MM F an' V am XHminl I (μνη) F(fφ)=∑a(μνη)∑δa(αβγ)H I.l(αβγ), μ,ν,η=0 α,β,γ=0 .(μνη) r "^^s+t+g-3^[δl (μ-") + "][δi (α-") + "] ≈l(αβY)=L 2 J s × [δ2 (ν - ") + "] [δi2 (β - ") + "] [δ3 (η - ") + "] [δ3 (γ - ") + "] получаем следующие зна- (20) li(γ) где ν!β! η!γ! (2") " g-" [δ3 (η-")+"][δ3 (γ-")+"] (22) (23) (24) tg Для того чтобы перейти от МЖ при коэффициентах b(nμνη) и δb(mαβγ) к обычной МЖ относительно U(p^qr^)δu(stg), что равносильно переходу от одной системы координатных функций к другой, воспользуемся тензорами преобразования координатных функций. Для этой цели выпишем отдельно аппроксимации полей перемещений внутри КЭ по степенным функциям и по функциям с использованием Влияние аппроксимирующих функций при построении матрицы жесткости 31 одномерных полиномов Лагранжа: MNL Un = ∑∑∑ anm1„, m1=1m2=1m3=1 (25) MNL un = ∑∑∑ U',1P, P. )p1Lp2 L"∙, p1=1 p2 =1 p3 =1 (26) где Lp1 , Lp2 , Lp3 - одномерные полиномы Лагранжа, определяемые по формуле δΠ=1(x(j)-x(jδ) LPk = rm ^(k) f x( j ^-Xi )1∏f X( )-Хп+5((^)"| I (Pk)Jr=1l (Pk) (r) (Pk)J На основании (25) и (26) имеем следующую зависимость для тензора преобразования координатных функций: (27) ∂(mk -1)LPk ∂ k mk-1 Cmpk (mk -1)!= m^ k ' (∂χk} x1=x2=x3=0 а связь между коэффициентами разложения имеет вид an' = CP1P2P3 un' (m1m2m3) m1m2m3 p1p2p3 (28) (29) где Cmp1mp2mp3 представляет собой диадное тензорное произведение тензоров согласно формуле Cp1p2p3=Cp1∙Cp2∙Cp3=Cp. (30) m1m2m3 m1 m2 m3 m При этом имеет место очевидное равенство Cp1=Cp2=Cp3 . (31) m1m2m3 Коэффициенты МЖ Fmn, при a(μνη) и δ(a(αβγ)) вычисляются через коэффициенты МЖ Pi j при u(i pqr)δu(js tg) путем двойного умножения на тензоры преобразования P = cm' cn'^ , (32) i j i j m п Выведем коэффициенты тензора преобразования для различных законов аппроксимации. Для линейного закона аппроксимации индексы, входящие в (25) и (27), приобретают следующие значения m = n = l =2, mk =1,2, Pk =1,2, χ(1) = - 2, χ(2) = 2 (33) Так как тензоры преобразования одинаковы для всех направлений х1, х2, х3, то достаточно вычислять их для одного направления, например, Cmp1 при j = 1 для х1. В этом случае выражение (27) для полинома Лагранжа с учетом (33) при- 32 Р.В. Киричевский, А.В. Скринникова нимает вид -x(11))(x1-x(12)) Л 2 (χ' (х' - χ1p1 ))⅛1)- x('1) +δ(1P1 ))(⅛1)- χ12)+δ(P1)) . (+U - ¾ . (34) (х' -x(P1))(^1P1) + 2 + δ(1P1) Jχ1p1) -2 + δ(P1)} Используя формулу (28), вычисляем значения компонент тензора преобразования C11 = L1 x1=0 1 x1 = 0 x1 =0 1, 2, C1=∂L1 C2 =∂x1 x1 = 0 C2=L2 C1 =Lx1 x1 =0 =-1, 0 =1, =2, (35) x1 =1. 0 Таким образом, имеем следующую матрицу тензора преобразования: [Cm ]= C1 C12 ] c2 C22 _ 1 2 -1 (36) Для квадратичного закона аппроксимации (27), имеют такие значения: M=N=L=3,mk=1,2,3,pk=1,2,3, индексы, входящие в (25) и x(1) = 2, x(2) = 0, x('3) = 2, а одномерный полином Лагранжа суть: (x1-x(11))(x1-x(12))(x1-x(13)) (37) (х' - χ1p1 ))(χ1p1)- χ11) +δ(1P1 ))(χ1p1) - xl2) +δ(P^))(χ∙1p1) - χ∙(3) +δ(Pї)) ■ +δ(3) ) 2 +δ(τ^1) (38) (Х' - χ1p1 ))(⅛) + ⅛ + δ(1P1 ))⅛1) +δ(2p1))dP1) -І По формуле (28) получаем аналогично (35) следующую матрицу коэффициентов тензора преобразования: [Cm1]= Г C^ 1 C2 1 C3] 1 Г0 1 0] C^ C2 2 C *^2 = -1 0 1 C^ [3 C2 3 C *^3 -4 2_ (39) Влияние аппроксимирующих функций при построении матрицы жесткости 33 Для кубического закона аппроксимации значения индексов, входящих в (25) и (27), следующие: M=N=L=4,mk=1,2,3,4,pk=1,2,3,4; x1=-1,x1=-1,x1=1,x1=1, (40) ^(1)= 2, x(2 )= 6, x(3)= 6, x(4 ) = 2, а одномерный полином Лагранжа имеет вид Lp1 =R4 = L1 = R1(p1) = - x(11) )(x1 - x(12)) × (χ--і](χχ і: (χp.6+δ⅛>I 41-2 ÷δffi) Аналогично (35) получаем следующий вид матрицы коэффициентов тензора преобразования в случае кубического закона аппроксимации: 9 16 27 8 -9 -4 27 2 [Cm: ]= C1 1 C1 2 C1 3 C1 4 C2 1 C2 2 C2 3 C2 4 C3 1 C3 2 C3 3 C3 4 C4] 1 C 4 ’-"2 C 4 ’-"3 C 4 ’-"4 1 16 1 8 9 4 9 2 9 16 27 8 -9 -4 27 2 1 16 1 8 9 4 9 2 (41) Зная коэффициенты тензора преобразования для линейного (36), квадратичного (39) и кубического (42) законов аппроксимации перемещений по формуле (30), получаем общий тензор преобразования для коэффициентов МЖ КЭ. Такой подход к формированию коэффициентов МЖ при законах аппроксимации высоких порядков позволяет значительно сократить количество вычислений и легко строить МЖ для любых других видов аппроксимирующих полиномов. Рассмотрим влияние типа аппроксимирующих функций для построенных МЖ КЭ на примере расчета резинового элемента конструкции в виде полого толстостенного цилиндра, находящегося под внутренним давлением, на скорость сходимости полученного приближенного решения МКЭ. 34 Р.В. Киричевский, А.В. Скринникова Толстостенные цилиндры широко используются в химической, нефтяной, военной промышленности, на атомных электростанциях. Они подвергаются высоким давлениям и температурам, которые имеют постоянный или циклический характер. Анализ конечных радиальных и ортогональных перемещений при внутреннем давлении в толстостенных цилиндрах проводят с целью обнаружения предела прочности материала, поскольку промышленные цилиндры часто подвергаются давлению выше предела текучести материала. Следовательно, необходим точный анализ напряженно-деформированного состояния цилиндров, чтобы полностью использовать их несущую способность и обеспечить безопасность эксплуатации. Исходные данные: внутренний радиус R1 = 0.25 м, внешний радиус цилиндра R2 = 0,5 м, модуль сдвига μ=0.7 МПа, коэффициент Пуассона ν=0.49999, давление Р = 0.1 МПа. Расчетная схема представлена на рис. 2, где r ∈ [Λ1; Λ2]. Перемещения внутренней и внешней поверхностей определялись по формуле [10]: (-ν)+(1+ν)R22J (1 -ν)r . Рис. 2. Расчетная схема полого цилиндра под внутренним давлением Fig. 2. Calculation scheme for a hollow cylinder under internal pressure Расчет внутренних радиальных, u1 (при r = 0.25), и внешних радиальных, u2 (при r = 0.5), перемещений выполнен в рамках вычислительного комплекса «МІРЕЛА+» [1]. Результаты представлены в таблице. Здесь М1 - дискретизация по дуге, М2 - вдоль радиуса, М3 - по толщине цилиндра. Как видно из таблицы, увеличение степени аппроксимации перемещений в радиальном направлении способствует улучшению численных результатов при одинаковых сетках разбиения (точное решение: u1∕ u2 = 0,07738’0,04762). Скорости сходимости МКЭ для линейного, квадратичного и кубического законов аппроксимации перемещений проиллюстрирована на рис. 3. Как видно, увеличение степени дискретизации влечет более высокую точность решения. Кубическая аппроксимация дает более быструю сходимость МКЭ и на сетке 2×13×2 дает более точный результат. Тот факт, что аппроксимирующие функции более высокого порядка дают более точные результаты и быструю сходимость МКЭ, также подтверждается результатами, приведенными в работе [6] для неэластомерных конструкций. Влияние аппроксимирующих функций при построении матрицы жесткости 35 Рис. 3. Влияние аппроксимирующей функции на скорость сходимости МКЭ Fig. 3. Effect of the approximating function on the FEM convergence rate Результаты расчета полого цилиндра под внутренним давлением Закон аппроксимации Размеры сетки M1×M2×M3 Количество уравнений Перемещения u1∣ u2 Абсолютная погрешность решения Линейный 2×3×2 36 0.0317/0.0524 - 2×4×2 48 0.0481/0.0678 - 2×5×2 60 0.0612/0.0728 - 2×7×2 84 0.0684/0.0915 0.00898/0.0335 2×9×2 108 0.0711/0.0982 0.00628/0.0268 2×10×2 120 0.0752/0.1049 - 2×11×2 132 0.0742/0.1093 - 2×13×2 156 0.0756/0.1154 0.00178/0.0096 Квадратичный 2×3×2 36 0.0428/0.0591 - '2×5×2 60 0.0623/0.0731 - '2×7×2 84 0.0722/0.0955 0.00518/0.0295 2×9×2 108 0.0728/0.1042 0.00458/0.0208 2×11×2 132 0.0760/0.1124 - 2×13×2 156 0.0791/0.1263 0.00172/0.0013 Кубический 2×4×2 48 0.0492/0.0684 - 2×7×2 84 0.0712/0.0963 0.00618/0.0287 2×9×2 120 0.0753/0.1215 0.00208/0.0035 2×13×2 156 0.0775/0.1244 0.00012/0.0006 Заключение На примере полого толстостенного цилиндра показано функционирование трех условно выделенных блоков САПР конструкций из эластомеров, основанной на МКЭ. Полученные в первом блоке МЖ КЭ для линейного, квадратичного и кубического законов аппроксимации перемещений были использованы во втором блоке при расчете радиальных перемещений МКЭ. Анализ полученных результатов в третьем блоке показал, что использование КЭ с аппроксимирующей кубической функцией позволяет ускорить сходимость МКЭ и получить более точные ре- 36 Р.В. Киричевский, А.В. Скринникова зультаты в сравнении с линейной в 10 раз и с квадратичной функцией в 2 раза по абсолютной погрешности. Таким образом, порядок аппроксимирующей функции при построении матрицы жесткости конечного элемента существенно влияет на скорость сходимости МКЭ, что доказывает перспективность использования аппроксимирующих функций более высоких порядков для определения параметров напряженно-деформированного состояния изделий из эластомеров, в частности толстостенных цилиндров. При выборе аппроксимирующих функций для вычисления МЖ КЭ необходимо учитывать, что количество опорных точек интегрирования должно быть выбрано рациональным, так как от них зависит объем вычислений при интегрирова-3 нии, который пропорционален в трехмерных задачах n .

Скачать электронную версию публикации