Автоматихация производства ЭВА (1075779), страница 23
Текст из файла (страница 23)
Эти узловые условия приводят к системе уравнений, решив которую получим:
1= Фi; 2 =(4Фi-3Фj-3Фk)/L; 3 =2(Фi-2Фj+Фk)/L2.
|
|
| Рис. 16.1 |
Подставляя выраженные через узловые значения коэффициенты в (16.1) и перегруппировывая члены, получим интерполяционный полином для квадратичного комплекс – элемента в матричном виде:
= [Ni Nj Nk] { Фi Фj Фk }T (16.2)
здесь: Ni = (1-2x/L) (1-x/L); Nj = (4x/L) (1-x/L); Nk = (x/L) (1-2x/L).
Существуют и широко применяются на практике формальные способы вычисления функций формы, использующие их известные свойства.
Применение квадратичного элемента. Элементы высокого порядка применяются так же, как симплекс – элементы, поскольку выбор интерполяционного полинома не связан с исходными дифференциальными уравнениями. В качестве примера рассмотрим одномерную задачу переноса тепла. Ее решение было рассмотрено в разделе 12 с использованием симплекс – элементов. Задача касалась определения распределения температуры по длине стержня, подверженного конвективному теплообмену. Исходные уравнения для произвольного элемента, выведенные ранее, имеют вид:
[K(e)] {T} = {F(e)},
где:
[K(e)] = V [B(e)]T [B(e)]dV + S2 h [N(e)]T [N(e)] dS
[F(e)] = - s1 q [N(e)]T dS + S2 h TOC[N(e)]T dS
Все интегралы должны быть вычислены заново, если мы используем квадратичный (вместо линейного) элемент. С этой целью запишем интерполяционный полином, аппроксимирующий температуру во внутренних точках комплекс – элемента:
T = [Ni Nj Nk] { Ti Tj Tk }T
Матрица градиентов с учетом выражения для функций формы примет вид:
| B(e)] = | dNi | dNi | dNi | = | (4x–3L) | (4L–8x) | (4x–L) | |||||
| dx | dx | dx | L2 | L2 | L2 |
Вычисляем объемную часть матрицы теплопроводности, полагая dV=S(1)dx:
| S(e) | (4x–3L) 2 | (4x–3L) (4L–8x) | (4x–3L) (4x–L) | dx | ||||
| L4 | L4 | L4 | ||||||
| (4L–8x) 2 | (4L–8x) (4x–L) | |||||||
| L4 | L4 | |||||||
| (4x–L) 2 | ||||||||
| Симметрично | L4 | |||||||
Вычисляем подробно интеграл для элемента а11 матрицы под интегралом:
| L | |||||||||||||||
| | 16x2-24xL+9L2 | = | 16x3 | L | + | 9x | L | - | 24x2 | L | = | 14 | |||
| L4 | 3L4 | 0 | L2 | 0 | 2L3 | 0 | 6L | ||||||||
| 0 |
Вычисляя аналогично остальные интегралы, приходим к выражению для объемной части матрицы теплопроводности элемента:
| 14 | -16 | 2 | |||||
| [KV(e)] = | (S(e)/6L(e)) | -16 | 32 | -16 | (16.3) | ||
| 2 | -16 | 14 |
Конвективная составляющая матрицы К(е) вычисляется по формуле:
| hР(e)L | NiNi | NiNj | NiNk | dx | ||
| NjNi | NjNj | NjNk | ||||
| NkNi | NkNj | NkNk |
Здесь Р(e) – периметр элемента. Подставляя выражения для функций форм и проводя интегрирование, получим:
| 4 | 2 | -1 | |||||
| [KS2(e)] = | (hLP(e)/30) | 2 | 16 | 2 | (16.4) | ||
| -1 | 2 | 4 |
Конвективная составляющая матрицы Fh(е) вычисляем по формуле:
| Fh(е) =hTOCP(e)L | 1 – (3x/L) + (2x2/L2) | dx = | 1 | |||
| (4x/L) + (4x2/L2) | (hLTOCP(e)/6) | 4 | ||||
| – (x/L) +(2x2/L2) | 1 |
Если конвективный теплообмен наблюдается на конце элемента, например, в узле i, то Ni=1, Nj=Nk=0 и поверхностный интеграл принимает вид:
hTOC Ai(e) [1 0 0]T
где: Ai(e) – площадь поверхности элемента в узле i. Наличие теплообмена в узле i сказывается и на матрице теплопроводности [K(e)], благодаря поверхностному интегралу по S2:
| hS2 | NiNi | NiNj | NiNk | 1 | 0 | 0 | |||
| NjNi | NjNj | NjNk | dx = h Ai(e) | 0 | 0 | 0 | |||
| NkNi | NkNj | NkNk | 0 | 0 | 0 |
Вычисление составляющей вектора нагрузки, обусловленной действием в 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] /63,75 [см] = 3,2 [Вт/oК]
(hLP(e)/30) = 10 [Вт/(см2 oC)] 3,75 [см] 2 [см] /30 = 2,5 [Вт/oК]
(hLTOCP(e)/6)= 10 [Вт/(см2oC)] 3,75 [см] 40[oC]2 [см] /6 = 500 [Вт]
(h S(e)TOC)= 10 [Вт/(см2oC)] 2 [см2]40[oC] = 400 [Вт]
Матрица теплопроводности для 1-го элемента определится суммой:
| 14 | -16 | 2 | 4 | 2 | -1 | |||||||
| [K(1)] = | 3,2 | -16 | 32 | -16 | + 2,5 | 2 | 16 | 2 | ||||
| 2 | -16 | 14 | -1 | 2 | 4 |
или:
| 54,8 | -46,2 | 3,9 | |||||
| [K(1)] = | | -46,2 | 142,4 | -46,2 | |||
| 3,9 | -46,2 | 54,8 |
Матрица теплопроводности для 2-го элемента определится суммой:
| 0 | 0 | 0 | |||||
| [K(2)] =[K(1)] + | Sk(2) h | 0 | 0 | 0 | |||
| 0 | 0 | 1 |
Матрица [K(2)] содержит дополнительное слагаемое, так как на свободном конце 2-го элемента (в k-ом узле) также происходит теплообмен :
Для вектор – столбца {F(1)} имеем:
| 1 | 500 | |||
| [F(1)] = | 500 | 4 | = | 2000 |
| 1 | 500 |
Для вектор – столбца {F(2)} имеем:
| 0 | 500 | |||
| [F(2)] =[F(1)] + | 400 | 0 | = | 2000 |
| 1 | 900 |
Объединяя полученные матрицы по методу прямой жесткости, получаем следующую систему уравнений:
| 54,8 | -46,2 | 3,9 | 0 | 0 | | Т1 | = | 500 | ||
| -46,2 | 142,4 | -46,2 | 0 | 0 | Т2 | 2000 | ||||
| 3,9 | -46,2 | 109,6 | -46,2 | 3,9 | Т3 | 1000 | ||||
| 0 | 0 | -46,2 | 142,4 | -46,2 | Т4 | 2000 | ||||
| 0 | 0 | 3,9 | -46,2 | 64,8 | Т5 | 900 |
Так как Т1=150 оС задана, то как и ранее получим модифицированную систему:
| 54,8 | 0 | 0 | 0 | 0 | | 150 | = | 8220 | ||
| 142,4 | -46,2 | 0 | 0 | Т2 | 8930 | |||||
| 109,6 | -46,2 | 3,9 | Т3 | 415 | ||||||
| 142,4 | -46,2 | Т4 | 2000 | |||||||
| Симметрично | 64,8 | Т5 | 900 | |||||||
Решая данную систему, получим следующие значения температуры в узловых точках: {T}T=[150 80,8 55,8 46,3 43,5] (oC). Эти значения хорошо согласуются с вектором: [150 80,9 55,4 46,2 43,3], представляющим аналитическое решение исходной задачи.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
