Курс лекций - Математическое моделирование технических объектов (1075784), страница 17

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 17)

С этой целью запишем интерполяционный полином, аппроксимирующий температуру во внутренних точках комплекс – элемента:T = [Ni Nj Nk] { Ti Tj Tk }TМатрица градиентов с учетом выражения для функций формы примет вид:B(e)] =dNidxdNi dNidx dx(4x–3L)L2=(4L–8x)L2(4x–L)L2Вычисляем объемную часть матрицы теплопроводности, полагая dV=S(1)dx:(4x–3L) 2L4S(e)(4x–3L) (4L–8x)L4(4x–3L) (4x–L)L4(4L–8x) 2L4(4L–8x) (4x–L)L4dx(4x–L) 2L4СимметричноВычисляем подробно интеграл для элемента а11 матрицы под интегралом:L16x2-24xL+9L2=L416x33L4L+0L9xL2024x22L3L=0146L0Вычисляя аналогично остальные интегралы, приходим к выражению для объемной части матрицы теплопроводности элемента:142-16[KV(e)] =(S(e)/6L(e))-1632-162-1614(16.3)Конвективная составляющая матрицы К(е) вычисляется по формуле:hР(e)LNiNiNiNjNiNkNjNiNjNjNjNkdxNkNi NkNjNkNkЗдесь Р – периметр элемента.

Подставляя выражения для функций форм и проводяинтегрирование, получим:(e)89[KS2(e)] = (hLP(e)/30)42-12162-12(16.4)4Конвективная составляющая матрицы Fh(е) вычисляем по формуле:1 – (3x/L) + (2x2/L2)Fh(е) =hTOCP(e)L(4x/L) + (4x2/L2)1dx =(hLTOCP(e)/6)41– (x/L) +(2x2/L2)Если конвективный теплообмен наблюдается на конце элемента, например, в узле i, то Ni=1, Nj=Nk=0 и поверхностный интеграл принимает вид:hTOC Ai(e) [1 0 0]Tгде: Ai(e) – площадь поверхности элемента в узле i. Наличие теплообмена в узле i сказывается и на матрице теплопроводности [K(e)], благодаря поверхностному интегралупо S2:NiNiNiNjNiNk100hS2NjNiNjNjNjNkdx = h Ai(e)000NkNi NkNjNkNk000Вычисление составляющей вектора нагрузки, обусловленной действием в i-м узле теплового потока q (составляющая Fq(е)), аналогично вычислению конвективной составляющей вектора нагрузки Fh(е), поэтому можно сразу записать:Fq(е) = S1 q[N(e)]TdS = q Ai(e) [1 0 0]TПример 16.1.

Определить распределение температуры в стержне кругового сечения, изображенном на рисунке 16.2.Рис. 16.2Разбиение на конечные элементы показано в нижней части рисунка 16.2.1.Запишем матрицы теплопроводности для 1-го конечного элемента, для чего первоначально вычислим коэффициенты в матрицах (16.3) и (16.4):(S(e)/6L(e)) = 72 [Вт/(см oC)]  [см2] /63,75 [см] = 3,2  [Вт/oК]90(hLP(e)/30) = 10 [Вт/(см2 oC)] 3,75 [см] 2 [см] /30 = 2,5  [Вт/oК](hLTOCP(e)/6)= 10 [Вт/(см2oC)] 3,75 [см] 40[oC]2 [см] /6 = 500  [Вт](h S(e)TOC)= 10 [Вт/(см2oC)] 2 [см2]40[oC] = 400  [Вт]Матрица теплопроводности для 1-го элемента определится суммой:[K(1)] =3,214-162-1632-162-1614+ 2,542-12162-124или:[K(1)] =54,8-46,23,9-46,2142,4-46,23,9-46,2 54,8Матрица теплопроводности для 2-го элемента определится суммой:0 0 0[K(2)] =[K(1)] + Sk(2) h0 0 00 0 1Матрица [K(2)] содержит дополнительное слагаемое, так как на свободном конце 2-го элемента (в k-ом узле) также происходит теплообмен :Для вектор – столбца {F(1)} имеем:1[F(1)] =500  41(2)Для вектор – столбца {F } имеем:0(2)(1)[F ] =[F ] +400  0= 5002000500= 15002000900Объединяя полученные матрицы по методу прямой жесткости, получаем следующуюсистему уравнений:54,8-46,23,900-46,2142,4-46,2003,9-46,2109,6-46,23,900-46,2142,4-46,2003,9-46,264,8Т1Т2Т3Т4Т5=500200010002000900Так как Т1=150 оС задана, то как и ранее получим модифицированную систему:54,80142,4Симметрично0-46,2109,600-46,2142,4003,9-46,264,8150Т2Т3 =Т4Т582208930415200090091Решая данную систему, получим следующие значения температуры в узловых точках:{T}T=[150 80,8 55,8 46,3 43,5] (oC).

Эти значения хорошо согласуются с вектором: [15080,9 55,4 46,2 43,3], представляющим аналитическое решение исходной задачи.Интересно отметить одну особенность, касающуюся полученных матриц квадратичных элементов: поверхностный интеграл в матрице теплопроводности (16.4) содержит отрицательные коэффициенты, чего не было в случае линейного элемента и что является обычным делом при использовании элементов высокого порядка.

Есть и другие особенности. Этоговорит о том, что бессмысленно предугадывать результаты интегрирования, когда мы имеемдело с элементами высокого порядка.17.1 Вывод уравнений элементов методом ГалеркинаМетод Галеркина является приближенным методом решения дифференциальных уравнений первого и второго порядка. Преимуществом его использования является то, что он исключает необходимость вариационной формулировки задачи, поскольку оперирует непосредственно с самим дифференциальным уравнением  () = 0. Решение задачи методом конечных элементов в контексте использования метода Галеркина начинается непосредственнос записи общего вида искомой системы уравнений вида:ne=1 L[N(e)] T  () dx = 0(17.1)где: n – общее число конечных элементов; L – верхний предел интегрирования, равный длине одномерной области, в которой производится поиск решения.Обязательным условием построения системы (17.1) является то, что в неё могут включаться производные порядка не выше первого.

Таким образом, если исходное дифференциальное уравнение  () = 0 имеет первый порядок, то никаких дополнительных выкладок привыводе системы уравнений для конечных элементов делать не нужно. Рассмотрим вначалеименно такое уравнение.Осевое нагружение стержняИзвестно, что перемещение (u [см]) стержня с площадью поперечного сечения S[см2],подверженного осевому нагружению (силой F [H]), описывается дифференциальным уравнением 1-го порядка вида:duF= 0(17.2)dxAEгде: Е – модуль упругости [H/см2].Рассмотрим методику составления матриц элементов и системы линейных уравнений для (17.2), взяв исходные данные из примера 13.1.1.

Разбиваем стержень на 2 элемента, как показано на рисунке 13.3.2. Запишем ДУ (17.2) для e–го (е=1,2) элемента:du(е)dx-F(е)A E= 0(17.3)3. Подставляя (17.3) в (17.1), получим систему уравнений 1-го элемента:92du(1)F(17.4)) dx = 0[N]((1)LdxS E4. Заменив в (17.4) функцию перемещений u(x) ее интерполяционным полиномом первого порядка, получим:d([N(1)] { U(1)})F(1) T(17.5)) dx = 0((1)L [N ]dxS EЗаписываем выражения для матриц ФФ и градиентов:(1) TN (1) = [ (1-x/L) (x/L)] ;B (1) = 1/L [ -1 1]Подставляем найденные матрицы в (17.5) и выполняем перемножение:1L21L2 (L-x)(1)dx {U } -xL ([ -1 1]L-x(L-x)х-х)dx {U(1)} -FS EL(1)dx =0xL(x-L)LF(1)S ELL-xdx =0 (17.6)xLВычисляем промежуточные интегралы:LLx dx =0L-(L-x) dx = Lx0L=0L220LL-x dx=0-11(x-L) dx = -L220Подставляем результат в (17.6):11U1-12x22{U2FL2S(1)E}-1{ 1}=0Величина FL/2S(1)E согласно (13.5) равна: 1,25 мм, следовательно, искомая система уравнений для первого элемента такова:1U11,25-1==01U21,25-1{}{}Аналогичные вычисления для 2-го элемента приводят к системе:-1-111{U2U3}={0,6250,625}=0Объединяя обе системы по методу прямой жесткости, приходим к системе:10U11,25-101U21,875-1=01U30,625-1{}93Поскольку U1 = 0, из первого уравнения получаем: -U1 + U2 = 1,25.

То есть U2 = 1,25мм. Тогда из второго уравнения имеем: : -U1 + 0 U2 + U3 = 1,875. То есть U3 = 1,875 мм. Видим, что результат совпадает с перемещениями, полученными ранее в примере (13.1).Изгиб консолиРешим теперь дифференциальное уравнение второго порядка (13.47), описывающееупругую линию консоли, неподвижно закрепленной на одном конце и подверженной действию перпендикулярной к ее оси силы – на свободном конце. Исходные данные для этой задачи подробно описаны в разделе 13.6. Поэтому можно сразу приступить к ее решению методом Галеркина, приняв уравнение (13.47) в качестве исходного дифференциального уравнения.

То есть, в соотношении (17.1) имеем:2yx2My(0) = 0;=  (y);EJ1. Записываем уравнения метода в общем виде:2yMn(e) Tdх = 02e=1 L [N ]xEJ-()(17.7)Здесь n – число элементов; L – длина отдельного элемента.Прежде, чем начинать вычисления, необходимо: (1) выбрать ФФ и (2) преобразоватьполученный интеграл так, чтобы он содержал производные порядка не выше первого (!?только в этом случае мы получим систему линейных уравнений для решения ?!).1.

Чтобы иметь возможность сравнить результаты расчетов с результатами, полученными в разделе 13.6, выберем то же разбиение консоли на конечные симплекс – элементы,представленное на рисунке 13.15.2. Запишем интерполяционную модель упругой линии:y = Ni(e) Yi + Nj(e) Yj = [(1 –x/L) (x/L)] [Yi Yj]T = N (e) [Y]3. Кривизна консоли =M/EJ - функция длины элемента. Аппроксимируем ее с помощью линейной модели: = N (e) [i j] T(17.8)4. Избавиться от производной второго порядка в уравнении (17.7) можно, взяв по частям интеграл: [N(e) TL](2yx2) dх(17.9)Обозначим v=(dy/dx), u=[N(e)] T и по формуле интегрирования по частям:XjXju dv = uv -XiXjv du =dydx-[N(e)] TXiXjXidydx[B(e)] TdxXiВ теории доказывается, что первое слагаемое здесь учитывается только в том случае,если на концах элемента определены производные.

Характеристики

Список файлов лекций