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

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

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

Для его вычисления вспомним правило перемножения двух матриц на примере:a bax+buay+bvaz+bwx y zc d cy+dvcz+dw= cx+duu v we fex+fuey+fvez+fwИскомое произведение матриц примет вид:-4400000-40400-404-400040000=16-160000-16320-16000000000-1601600000000000000Принимаем толщину элемента, равной единице, и выносим произведение матриц за знак интеграла (13.32). В результате dV становиться равным 1/32. Следовательно, МЖЭ для первого конечного элемента:К(1) = ½1-10000-120-1000000000-10100000000000000Для вычисления объемного интеграла (13.33) в матрице нагрузки воспользуемся интегральной формулой вычисления площади треугольника с применением системы L –координат. Последняя представляет собой совокупность трех относительных координат точки R внутри треугольника L1, L2 и L3, каждая из которых является отношением расстоянияот точки А до одной из сторон треугольника.

Если A – площадь треугольника, то L – координата точки R относительно стороны (i, j) равна L1=А1/A. Можно показать, что переменныеL1, L2 и L3 представляют собой ФФ для треугольного симплекс – элемента. Это обстоятельство вкладывает дополнительный смысл в понятие ФФ треугольного элемента. Например,для рисунка 13.9 имеем: L1= Nk.69Рис. 13.9Преимущество L – координат проявляется при необходимости вычисления интегралов вдоль сторон конечного элемента и по его площади. Так, если A, B, и С – целые числа, то справедлива следующая интегральная формула:A! B! C!ABCLLLdA=2A(13.38)123(A+B+C+2)!АЧтобы воспользоваться формулой (13.38) необходимо подынтегральное выражение выразить через L – координаты. Вычислим в качестве примера следующий интеграл по площади А: Ni Nj dA =  (L11 L21 L30) dA = { 1!1!0!/(1+1+0+2)! }2A = 2A/4!=A/12.В текущем примере надо вычислить интеграл (13.33).

Вычислим вначале: Ni dA =  (L11 L20 L30) dA = { 1!0!0!/(1+0+0+2)! }2A = 2A/3!=A/3.Запишем интеграл (13.33) в развернутом виде для первого элемента:N1(1)N2(1)f(1) =0dV2G(1)(13.39)(1)N4V00(1)(1)Заменяя ФФ L – координатами: L1=N1 ,L2=N2 , L3=N4 (1), имеем:(1)f(1) = 2G A [ 1 1 0 1 0 0 ]T3Подставляя численные значения, имеем: f(1) = [29 29 0 29 0 0]TТаким образом, результирующая система уравнений для первого конечного элементапримет вид:[K(1)] {Ф} = {f(1)}½1-10000-120-1000000000-10100000000000000Ф1Ф2Ф3Ф4Ф5Ф6=292902900(13.40)Системы уравнений для второго, третьего и четвертого элемента примут вид:70[K(2)] {Ф} = {f(2)}½½½00000001-1000000000010-1000000000000000-120-1000000000-10100000000000000-102-10000-110000000000000000-1-10000-12-10000-11[K(3)] {Ф} = {f(3)}[K(4)] {Ф} = {f(4)}Ф1Ф2Ф3Ф4Ф5Ф6Ф1Ф2Ф3Ф4Ф5Ф6Ф1Ф2Ф3Ф4Ф5Ф6===029290290(13.40-а)029029290(13.40-б)000292929(13.40-в)Окончательная система уравнений примет вид: [K] {Ф} = {f}½1-10000-14-1-2000-120-100-204-2000-1-24-10000-11Ф1Ф2Ф3Ф4Ф5Ф6=298729878729(13.40-г)Величины Ф3, Ф5 и Ф6 равны нулю, так как соответствующие им узлы расположены навнешней границе.

Преобразуя и решая систему, получим:{Ф}Т = {218 160 0 123,6 0 0}ТПредставляет особый интерес поверхность , соответствующая полученномумножеству узловых значений. Эта поверхность представлена на рисунке 13.10. Дело втом, что объем тела, ограниченный этой поверхностью, пропорционален не найденному пока моменту М, закручивающему стержень на один градус на длине 1 метр. Крометого, не определены еще сдвиговые напряжения.1. Вычисляем сдвиговые напряжения, связанные с найденной функцией напряжений формулами (13.27).

Согласно выражениям (13.27) и (13.30), запишем выражение для матрицы–столбца сдвиговых напряжений первого элемента, выраженное черезопределенную уже матрицу градиентов (13.37): (1) =xy= [B(1)]{Ф} =zy(1) = -x-4 4 0 0 0 00 –4 0 4 0 0= 232,6(н/см2); zy(1)=y2181600123,600=-232,6-145,4= -145,4(н/см2);71Компоненты тензора напряжений для других элементов вычислим аналогично:ZY(2) =640,0(н/см2); ZX(2) =0,0(н/см2);ZY(3) =494,0(н/см2); ZX(3) = -145,4(н/см2);ZY(4) =494,0(н/см2); ZX(4) =0,0(н/см2);Полученные компоненты тензора напряжений схематически показаны на рисунке.13.11. Внутри каждого конечного элемента сдвиговые напряжения постоянны, так какинтерполяционные полиномы взяты линейными по Оx и Оy.Рис. 13.10Рис.

13.112. Вычисляем момент, который согласно определению равен:М=2 dSSгде: S – площадь поперечного сечения стержня.Учитывая аддитивный характер момента, можно записать выражение для частимомента, действующего в пределах рассматриваемого участка стержня:МО= М (е) =   2 dА(13.41)(е)е АеПоскольку площади всех четырех конечных элементов равны, обозначим:А(1)=А(2)=А(3)=А(4)=АВеличины (е) при известных функциях формы определяются по (13.29).Для первого элемента имеем:(е)(е)М (1) = 2   dААОткуда:(1)(1)(1)= 2  [N1 N2 0 N4 0 0] {Ф}dAА72(1) ТМ (1) = 2 {Ф}Т  [N ] dААИнтеграл в полученном выражении идентичен вычисленному ранее интегралу(13.39).

Поэтому можно сразу записать:22М (1) =А (Ф1 + Ф2 + Ф4)А {Ф}Т [ 1 1 0 1 0 0] =33Подставляя значения , найденные в узловых точках, имеем:М (1) = 0,67 А (218 + 160 +123,63) = 0,67 А (501,8)Аналогично находим для остальных элементов:М (2) = 0,67 А (Ф2 + Ф3 + Ф5) = 0,67 А (160)М (3) = 0,67 А (Ф2 + Ф5 + Ф4) = 0,67 А (283,6)М (4) = 0,67 А (Ф4 + Ф5 + Ф6) = 0,67 А (123,6)Суммируя эти соотношения по формуле (13.40), получим:МО= 0,67 А (501,8 + 160 + 283,6 + 123,6) = 0,67 А (1069)Поскольку на элементы разбивалась только 1/8 часть поперечного сечениястержня и А=1/32, полный крутящий момент M составит:М = 8 МО = 8 (0,67)  (1/32)  (1069) = 178,16 (нсм)Это означает, что крутящий момент величиной 178 Хнсм) вызовет закручивание на 1о стального стержня длиной 1 м. Теоретическое значение связи между приложенным крутящим моментом и углом закручивания квадратного стержня со стороной,равной L (м), дается формулой:М = 0,1406GL4 =0,14068106[н/м2]1(рад/м)1(м)=196,3(нсм)Полученное расхождение в 9,5% объясняется выбором грубой сетки разбиения.13.4 Метод прямой жесткостиПрименение метода конечных элементов, как показано на ряде примеров, приводит к системе алгебраических уравнений.

Порядок системы уравнений совпадает счислом неизвестных узловых значений и может достигать сотен тысяч. Метод построения глобальной матрицы жесткости, рассмотренный в разделе 13.3, весьма неэффективен, если в таком виде пытаться реализовать его ЭВМ. Дело в том, что матрицажесткости отдельного конечного элемента [k(e)] (назовем ее локальной матрицей жесткости – ЛМ), имеет такое же число строк и столбцов, что и глобальная матрица жесткости (ГМ) [K], однако, большинство коэффициентов ЛМ равно нулю. Действительно, на рисунке 13.12 показана область, разбитая на 16 элементов и имеющая 15 узловых точек. Матрица жесткости для первого элемента этой области, показанная на ри-73сунке 13.13, содержит 15 х 15 = 225 коэффициентов. Видно, что только девять из них ненулевые, что составляет 4% от общего числа ячеек, занимаемых в ОЗУ матрицейЛМ.Рис.13.12Рис.13.13С ростом числа узлов процент нулевых элементов резко увеличивается.

Например, для75 узлов он составляет (75 х 75 –9)/(75 х 75) 100% = 99,84%. Кроме того, матрица элемента[k(e)] должна быть вычислена отдельно от глобальной матрицы [K] и затем прибавлена к последней. Это, в свою очередь, требует запоминания в ОЗУ обеих матриц: и [k(e)] и [k(e)].В эффективных программах процедура построения ГМ использует сокращенную форму локальных матриц элементов [N(e)] при получении уравнений для элементов. Такой методизвестен как метод прямой жесткости (МПЖ).Рассмотрим процедуру кодирования ЛМ в методе МПЖ. С этой целью перепишемсистему интерполяционных полиномов (13.29), исключив из неё все глобальные степени свободы , не относящиеся к конкретному элементу:(1)(2)(3)(4)====N1(1)Ф1N2(2)Ф2N2(3)Ф2N4(4)Ф4++++N2(1)Ф2N3(2)Ф3N4(3)Ф4N5(4)Ф5++++N4(1)Ф4N5(2)Ф5N5(3)Ф5N6(4)Ф6(13.42)ФФ записываются здесь в соответствии с порядком следования индексов i,j и kпротив часовой стрелки, начиная с узла, отмеченного звездочкой на рисунке (13.7).Будем называть такой порядок установленным.На основании системы (13.42) для каждого конечного элемента составляем локальные матрицы градиентов (13.30) (матрицы сдвиговых напряжений).Подробно для первого конечного элемента имеем:1) Матрица градиентов: (1) =xy(1)=1/(2A)(1)(1)b1 b2 b4c1(1) c2(1) c4(1)Ф1Ф2Ф4742) Локальная матрица жесткости для первого элемента [k(1)] записывается сфиксацией номеров строк и столбцов в установленном порядке:[k(1)] =110,52-0,5402-0,51-0,50-0,5 0,54Числа, отмечающие строки и столбцы этой матрицы представляют собой номераглобальных степеней свободы.

Применив подобную процедуру к интегралу (13.39),имеем:291f (1) =292294Таким образом, система локальных уравнений, описывающих первый конечныйэлемент, примет вид:1240,5 -0,5011Ф1292-0,5-0,51=Ф22920-0,5 0,544Ф429Полученные локальные матрицы [K(1)] и [f (1)] содержат 12 ячеек вместо 42,требуемых ранее.

Ее необходимо прибавить к ГМ ( предполагается, что элементы ГМпредварительно сброшены в 0):½1-10000-120-1000-10000000100000000000000Ф1Ф2Ф3Ф4Ф5Ф6=292902900При сложении ЛМ с ГМ устанавливаются строгие формальные правила, продиктованные сущностью метода прямой жесткости. Изложим эти правила для элемента Нлокальной матрицы жесткости, стоящего на пересечении ее i-й строки и j-го столбца,который прибавляется к ГМ:1) Прочитать в ячейку iG номер глобальной степени свободы, отмечающей i-юстроку ЛМ.2) Прочитать в ячейку jG номер глобальной степени свободы, отмечающей j-мустолбцу ЛМ.3) Прибавить к содержимому ячейки ГМ, расположенной на пересечении iG-йстроки и jG-го столбца, элемент Н.Сложив по этому правилу ЛМ [k(1)] с ГМ [K], приходим к полученной ранеесистеме (13.39).75Вычислив аналогично ЛМ для второго конечного элемента [k(2)] и [f(2)], получимпромежуточный вариант локальной системы:2350,5 -0,5022Ф2293-0,51-0,550-0,50,5=Ф3Ф5293295Сложив результат с ГМ [K], приходим к глобальной системе:½1-10000-13-1-1000-120-100-1010000-1010000000Ф1Ф2Ф3Ф4Ф5Ф6=29582929290Далее проводим аналогичные действия для третьего и четвертого элементов.Последовательно получаем:1) Локальная система уравнений для третьего элемента:2540,50-0,522Ф229500,5-0,54-0,5-0,51=Ф5Ф42952942) Глобальная система после добавления сюда локальной примет вид:½1-10000-14-1-2000-120-100-203-1000-1-120000000Ф1Ф2Ф3Ф4Ф5Ф6=298729585803) Локальная система уравнений для четвертого элемента:440,55-0,5605-0,51-0,560-0,50,5Ф4=Ф5Ф62942952964) Глобальная система после добавления сюда локальных матриц для четвертогопоследнего конечного элемента приобретает окончательный вид:½1-1000-14-1-200-120-10-204-200-1-240000-1Ф1Ф2Ф3Ф4Ф5=2987298787760000-11Ф629Как и следовало ожидать, мы получили ту же самую систему (13.40-г).13.5 Задача изгиба опертой балкиЗадача изгиба опертой балки с точки зрения метода конечных элементов представляет частный случай рассмотренной выше задачи о кручении бруса.Условие задачи.

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

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