Чайнов Н.Д. - Конструирование двигателей внутреннего сгорания (1037884), страница 12
Текст из файла (страница 12)
Искомая температураТ(x, y, z, t) представляется какрения соответствующего интегралавзвешенной невязки по области V,занимаемой теломòW p RdV(2.100)где R – невязка или погрешность,возникающая при подстановкеприближенного решения T в уравнение (2.2); р = 1, 2, …, n – числоузловых параметров; Wр – определяемые последовательно по области V интерполянты.При применении МКЭ в качестве Wр выбираются функции формыNр и уравнение (2.100) принимаетвидò N p RdV= 0.(2.101)VПодставляя выражение (2.100) вуравнение (2.2) с последующим определением невязки R и используяинтегральную формулировку выражения (2.101), приходят к системеобыкновенных дифференциальныхуравнений первого порядка относительно узловых значений температур Tp(t), которая в матричнойформе имеет вид[c]T ( x , y , z, t) = [N]{T e (t)}. (2.99)Наряду с примененной выше, вслучае стационарной задачи, вариационной формулировкой широко используется более универсальный невариационный подход.Весьма распространенным здесьявляется известный метод Бубнова–Галеркина, являющийся частным случаем метода невязок.
Приближенное решение уравнения(2.2) нестационарной теплопроводности получается на базе рассмот= 0,V¶{T }+[H]{T } = { f }, (2.102)¶tгде [c] – матрица теплоемкости,m[c] = å[c e ]; [c e ] = ò cr[N]T [N]dV .e =1VeИногда по аналогии с задачамидинамики матрицу [c] называютматрицей демпфирования.С помощью Lкоординат интеграл в выражении [ce] определяется по формулам (2.73) и (2.74).Для плоской и осесимметричнойзадач55é2 1 1 ùc rDt ê[c ] =1 2 1ú;ú12 êêë1 1 2 úû(2.103)é2 1 1 ùpr Dcr ê[c e ] =1 2 1ú.ú6 êêë1 1 2 úû(2.106) находят узловые температуры в момент времени (t + Dt). Данное решение предполагает посто¶Tв течениеянство производной¶tинтервала Dt.
Существуют и другиеболее совершенные способы решения. Так, если предположить, что впределах Dt вектор {T } температуринтерполируется по некоторым егозначениям, в частности, линейнопо двум значениям {T }0 и {T }1, соответствующим началу и концу интервала Dt, тоeНаиболее простым приближенным способом решения системыдифференциальных уравнений является переход от производнойтемпературы по времени к ее конечноразностному выражению¶T (T1 -T 0 )=,¶tDtT = [N 0(2.104)где Т0 и Т1 – соответственно начальные и конечные значения температуры в пределах временногоинтервала Dt.При применении центральнойразностной схемы величины {T } и{ f } определяются как среднее значение для интервала Dt, т.е.{T } 0 + {T }1 ü;ïï2ý{ f } 0 + { f }1 ï{f} =.ïþ2(Dt - t)t; N1 =; 0 £ t £ Dt.DtDtТакая постановка задачи предусматривает распространение метода конечных элементов и для координаты времени.
В этом случаегде N 0 =ì{T } 0 ü¶T æ 1 ö=ç÷ [-1 1] íý. (2.108)¶t è D t øî {T }1 þ{T } =(2.105)Подставив зависимости (2.107)и (2.108) в выражение (2.102) и умножив его на N1, после интегрирования в пределах Dt получим систему алгебраических уравнений относительно неизвестных значенийвектора температур {Т}1, значениявектора температур {Т}0 в началеинтервала при этом известныС учетом выражения (2.105) система уравнений (2.102) принимаетвид2 üìí[H] + [c]ý{T1 } =Dt þîæ 2ö= ç [c] -[H] ÷{T } 0 + 2{ f }.Dtèøì{T } üN 1 ]í 0 ý, (2.107)î {T }1 þ{T }1 = -{2 3[H] +[c] Dt} -1 ´(2.106)ì´í[[H]1 3 -[c] Dt]{T } 0 îDtüï2-[2 (Dt) ] ò { f }tdt ý.ïþ0Рассмотренная расчетная схемаявляется явной по времени.
Значения узловых температур {T}0 в момент времени t считаются известными и из системы уравнений56(2.109)Как и при решении задачи теплопроводности вклад в общую суммудадут лишь элементы, включающиеузел, по перемещению которого вычисляются частные производные.В результате получают системулинейных алгебраических уравнений относительно компонентов перемещений узлов. Для двухмерныхзадач порядок системы равен 2n,так как перемещение di каждого узла характеризуется двумя компонентами ui и vi, а для трехмерныхзадач порядок системы равен 3n,так как {d i }T = [u i v i w i ].Для плоской и осесимметричной задач напряженнодеформированное состояние в точке характеризуется соответственно следующими компонентами2.4.2.
Расчет напряженнодеформированного состояниядеталей двигателяПри расчете теплонапряженныхдеталей данному этапу предшествует определение теплового состояния, поэтому все, связанное с разбивкой детали на элементы и ихописанием, берется в этом случае вготовом виде. Если температура детали заведомо близка к постояннойпо объему (например, шатуны,подвеска коленчатого вала и др.),то тепловое состояние не рассчитывается, но требуется выполнитьвсе операции, связанные с разбивкой детали на элементы и их описанием.При решении задачи в перемещениях вектор перемещений точекв пределах элемента d выражаетсячерез компоненты перемещений узлов {de} элементов соотношениями(2.60–2.62).
Компоненты перемещений n узлов являются неизвестными задачи. Их определяют минимизацией функционала (2.28), выражающего полную потенциальнуюэнергию системы. При этом используется зависимость (2.87), в которой величины Ф и {T } следует заменить соответственно на П и {d}.Для минимизации функционалаП необходимо вычислить частныепроизводные от функционалов Пеотдельных элементов по перемещениям узлов, затем просуммироватьодноименные производные, т.е.производные для различных элементов, но по перемещениям одного узла, и результат приравнятьнулю.{s }T = [s xsy{e} = [e xeyTt xy ];g xy ];{d}T = [u v]и{s }T = [s rsz{e}T = [e rezsqeqt r z ];g r z ];{d} = [u v].TОставшиеся компоненты объемного напряженнодеформированного состояния в данном случаеравны нулю.
Так, для плоского напряженнодеформированного состояния sz = tyz = tzx = 0; gyz = gzx = 0.Так как задача решается в перемещениях, с помощью формул(2.21) компоненты деформацииэлемента {ee} выражают через компоненты перемещения, используя57é¶N iê ¶xêê 0êê 0ê[B i ] = ê¶Nê iê ¶yêê 0êê¶N iêë ¶zдля последних интерполяционнуюзависимость (2.60){e e } = [B]{d e },(2.110)где [В] – матрица деформации,получаемая дифференцированиемматрицы [ N] с учетом зависимостей (2.21).
Для линейного треугольного элемента [В] = [В1 В2 В3];для линейного тетраэдра [В] == [В1 В2 В3 В4].В случае плоской и осесимметричной задачé¶N iêê ¶x[B i ] = ê 0êê¶N iêë ¶yébi1 ê=02D êêëc ié¶N iê ¶rêê 0[B i ] = êê Niê rê¶N iêë ¶z0¶N i¶y¶N i¶xébiê0ê1 ê0=6V êc iêê0êdë iùúúú=úúúû0ùc i ú;úbi úû¶N i¶y0¶N i¶x¶N i¶z00ci0bidi00ù0úúdi ú,0úúci úbi úû(2.113)где i принимает значения 1, 2, 3, 4.Температурная деформация, рассматриваемая как разновидность начальной деформации, для плоской,осесимметричной и трехмерной задач может быть представлена в матричной форме соответственно следующим образом:{e e0 }T = a T T e [1 1 0];{e e0 }T = a T T e [1 1 1 0];{e e0 }T = a T T e [1 1 1 0 0 0].Для изотропного материаласвязь напряжений с деформациямивыражается зависимостями (2.24).Разрешая последние относительнонапряжений, окончательно получим(2.111)ù0 ú¶N i úú¶z ú =0 úú¶N i úú¶r ûbi0ùéê0ci ú1 êú,=aci zúêb0++2Diúêrrcibi úûêëù0 úú0 úú¶N i ú¶z úú=0 úú¶N i ú¶y úú¶N i ú¶x úû0(2.112){s e } = [D]({e e } - {e 0e }).(2.114)Матрица упругости [D] для плоской, осесимметричной и трехмерной задач соответственно имеетследующий вид:где i принимает значения 1, 2, 3.В случае трехмерной задачи58полной потенциальной энергиисистемы П в матричной формедля плоской задачиéù0 úê1 mE ê0 ú;[D] =m 11-m 2 ê1-m úê0 0ú2 ûëm æП = åç ò 0,5{d e }T [B]T [D][B]{d e }dV e =1 ç V eè(2.115)- ò {d e }T [B]T [D]{e 0e }dV Ve- ò {d e }T [N]T {P }dV -для осесимметричной задачиE (1 - m)[D] =´(1 + m)(1 - 2m)mméê 11-m 1-mêmmê11-m1-m´êê mm1êê1 - m 1 - mê 000êëVeö- ò {d e }T [N]T { pn }dF ÷ ÷Feø-{d e }T {R},ùúú0 úú ; (2.116)ú0 úú(1 - 2m) ú2(1 - m) úû0где {P} – вектор объемных сил;{pn} – вектор поверхностной распределенной нагрузки; {R} – вектор сосредоточенной нагрузки.При описании элемента индексе сохраняется только при неизвестных узловых перемещениях, хотякомпоненты матрицы деформации[B], векторы объемных сил и подля трехмерной задачиéê 1êê mê1 - mê mE (1 - m) ê1 - mê[D] =(1 + m)(1 - 2m) ê 0êêê 0êê 0êëm1-m1m1-mm1-mm1-m1000000(2.118)ùúú000 úúú000 úú.(1 - 2m)00 úú2(1 - m)ú(1 - 2m)00 ú2(1 - m)ú(1 - 2m) ú002(1- m) úû0Подставив в формулу (2.28)соотношения (2.110) и (2.114) ивыполнив суммирование по всемэлементам, получим выражение00(2.117)верхностной нагрузки вычисляются для каждого элемента.
При расчете теплонапряженных деталей,если требуется учесть зависимость59Е, m и aТ от температуры, компоненты вектора температурной деформации {e0} и матрицы упругости[D] также вычисляются для каждого элемента.Дифференцирование полнойпотенциальной энергии по вектору узловых перемещений производится аналогично дифференцированию функционала Ф(Т ) повектору {Т } с использованием зависимости (2.88). Минимизацияполной потенциальной энергиисистемы П приводит к соотношениюВ результате суммирования поотдельным элементам получимглобальную матрицу жесткости [k]и глобальный вектор нагрузки {G }üïïe =1ý (2.122)me ï{G } = {R} - å{G }.ïþe =1m[k] = å[k e ];Условие стационарности полной потенциальной энергии системы с учетом всех приведенных соотношений представляет системулинейных алгебраических уравнений относительно составляющихперемещений узлов, которые являются неизвестными задачи.
В матричной форме эта система имеетвид¶П m æç= å [B]T [D][B]dV {d e } ¶{d} e =1 ç Vòeè- ò[B]T [D]{e 0 }dV Ve[k]{d} = {G },- ò[N]T {P }dV Veгде {d} = [d1, …, dn] – вектор узловых перемещений.Заданные по условиям задачиперемещения некоторых узловыхточек (на осях симметрии, в местахзакрепления конструкции) непосредственно подставляются в систему (2.123) в качестве известныхкомпонентов вектора {d}.При использовании треугольного симплексного элемента в случаеплоской задачи интегрирование,связанное с определением матрицыжесткости [ke] и вектора нагрузки{G e }, осуществляется в замкнутомвиде, поскольку матрицы [B] и [D]содержат лишь постоянные величины. Интегрирование выражения,содержащего вектор объемных сил{P}T = [X Y], производится по объему элемента, а выражения, содержащие вектор распределенной нагрузки {pn} = [px py], по грани элемента, по которой эта нагрузкаТö- ò[N]T { pn }dF ÷ - {R} = 0.÷Veø(2.119)Вектор узловых перемещений{de} при интегрировании по объемуэлемента в локальной системе координат является постоянным.При этом важное значение имеетинтеграл[k e ] = ò[B]T [D][B]dV , (2.120)Veназываемый матрицей жесткостиэлемента.Остальные интегралы в соотношении (2.119) составляют векторнагрузки {G e } элемента, т.е.{G e } = - ò[B]T [D]{e 0 }dV Ve- ò [N] {P }dV - ò[N]T { pn }dF .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
