Диссертация (531291), страница 35
Текст из файла (страница 35)
Techn. Paper ASME 97-AA-020, Singapore ASME, New York, USA – 6p.284. Vorobjov Yu.S., Kanilo S.P. Sharing 3D finite element and beam model for turbomachine blades dynamic analysis // Multiple Scale Analysis and260Coupled Physical Systems. Sent–Venant Symposium. – Paris, 1997. – P. 469–473.285. Vorobjov Yu.S., Korsunsky M.L. The modify variational method ofcalculation ofturbomachinesbladingvibrations//Proc. Int.
Conf.Computational Eng. Sc. – 1995. – P. 1614–1618.286. Walker K.P. Vibrations of cambered helicoidal fan blades // J. Soundand Vibration. – 1978. – V. 59. – P. 35–57.287. Wolter I. Experimentelle Untersuchungen des Schwingungsverhaltensvon Turbinen-Laufdchaufeln unter realen Betriebsbedingungen mit und ohneKopplung durch eingelegten Dämpferdraht: Dis. … Dr. – Ing. – Stuttgart, UNI,1980.
– 179s.288. Yanecki S., Vorobiev Yu.S., Kanilo S.P. The numerical analysis of turbomachine blade packet vibration // там же. – P. 291–296.289. Yang B.D., Chen J.J., Meng C.H. Prediction of resonont responce ofshrouded blades with three-demensional shroud constraint // Trans. ASME. J.Eng. Gas Turbines and Power. – 199. – 121, №3. – P.
523–529.290. Yang T., Botkin M. A. A modular approach for threedementional shapeoptimization of structures // AIAA Journal. – 1987. – V.25, №3. – P. 492–497.291. Yang B.–D., Meng C.–H. Modeling of friction contact and it’s application to the design of shroud contact // Trans. ASME. J. Eng. Gas Turbines andPower. – 1997. – 119, №4. – P. 958–963.292. Zhang J., Wen—liang W., Xiang—jun C. Natural mode analysis of Nblades disc coupled system. Modal synthesis of symmetric structure with convcroup // Acta mech. solida sin., 1984. – N4. – P.
469–481.293. Zhu Baotion, Wu Houyu Xi’an Jiaotong daxue xuebao // J. Xi’an Jiaotong Univ. – 1999. – 33,№9. – P. 47–49.294. Косозубая зубчатая передача. / Айрапетов Э.Л., Апархов В.И.,Генкин М.Д., Федосеев Ю.Н., Авторское свидетельство СССР №442329,кл. F16H1/02, опубл. Б. И. № 33, 1974.261295. Косозубая зубчатая передача. / Айрапетов Э.Л., Апархов В.И.,Генкин М.Д., Косарев О.И., Авторское свидетельство СССР №700723, кл.F16H1/02, опубл.
Б. И. № 44, 1979.296. Цилиндрическая передача. / Айрапетов Э.Л.,Гребенников А.С.,Федосеев Ю.Н., Авторское свидетельство СССР №1698530А1, кл.F16H1/08, опубл. Б. И. № 46, 1991.297. Improvements in or relating to epicyclic gearings. / Stoekicht W.G.
Patentof Germany №1. 163. 181, F16 H1/00, 1966.298. Пинежанинов Ф., Пинежанинов П. Свойства базисных функций //Статьипометодуконечныхwww.exponenta.ru/soft/Mathemat/pinega/main.asp15.09.2013)262элементов.(датаобращения–:8 ПРИЛОЖЕНИЯПриложение 1Основные соотношения теории упругости и метода конечныхэлементовВ соответствии с методом конечных элементов, выбрав в качестве степеней свободы узловые перемещения, с помощью интерполяционной формулы (базисных функций) можно аппроксимировать перемещение любой точки внутри элемента. В технической литературебазисные функции принято называть функциями форм и обозначатьматрицу функций форм как [ N ] . Тогда интерполирующая формула дляперемещений {δ } запишется в виде:{δ } = [N ]{δ e },где индекс e (здесь и далее) указывает на принадлежность вектора илиматрицы элементу, т.е.
мы имеем дело с вектором узловых перемещений по всем степеням свободы КЭ.Для трехмерного случая перемещения раскладываются по осямx,y,z и обозначаются соответственно u,v,w, интерполирующая формула для n-узлового КЭ может быть записана в виде: u N1{δ } = v = 0 w 000N100N1N20000N200N2 u1 v 1L 0 w1 L 0 L = [ N ]{δ e } .L N n un vn w nСоотношения, связывающие деформации произвольной внутренней точки конечного элемента и узловые перемещения, записываютсяследующим образом:263 ∂u ∂N1 ∂x ∂x ∂v ε x ∂y ε ∂w y ε z ∂z {ε } = = = ∂Nγ∂u∂vxy + 1γ yz ∂y ∂x ∂y γ zx ∂v + ∂w ∂z ∂y ∂w + ∂u ∂N1 ∂x ∂z ∂zL∂N1∂y∂N n∂x∂N n∂yL∂N1L∂z∂N1∂x∂N1∂zL∂N n∂y∂N1L∂y∂N1∂N nL∂x∂z∂N n∂x∂N n∂z u1 v1 ∂N n w1∂z L u n v∂N n n ∂y wn ∂N n ∂x {ε } = [ B ]{δ e } ,илигде [B ] – матрица связи узловых перемещений и деформации.
В этойматрице не учтены члены высшего порядка малости, т.е. речь идет олинейной постановке задачи. При небольших перемещениях мы получаем вполне корректные результаты.Деформации {ε} и напряжения{σ} связаны соотношением:{σ } = [ D ]{ε } ,где [D ] – матрица упругости в общем случае зависит от свойств мате риала.Для изотропного материала, например, с модулем упругости E икоэффициентом Пуассона ν , матрица упругости примет следующийвид :264 1 ν (1 −ν ) νE(1 −ν ) (1 −ν )[ D] =(1 +ν )(1 − 2ν ) 0 0 0νν(1 −ν )(1 −ν )0000100001 − 2ν2(1 −ν )00001 − 2ν2(1 −ν )00001ν(1 −ν )ν(1 −ν )0 0 0 0 1 − 2ν 2(1 −ν ) 0Таким образом, имея функции форм, соответствующим образомпродифференцировав их , а также зная матрицу упругости материала ,можно записать основные соотношения теории упругости для конечного элемента:{δ } = [ N ]{δ e }{ε } = [ B ]{δ }{σ } = [ D ]{ε }с определенными выше матрицами [ N ] , [ B ] и [ D ] .Уравнения, связывающие узловые силы с узловыми перемещениями, можно записать относительно перемещений (уравнения жесткости) или относительно сил ( уравнения податливости):[ K ]{δ } = {F } или [ f ]{ F } = {δ } , где[ K ] – матрица жесткости,{F } – вектор узловых сил,[ f ] – матрица податливости.При имеющемся векторе перемещений {δ}, размерность которо го равна числу степеней свободы системы, потенциальная и кинетическая энергии запишутся в виде:265Ï∑1{δ } T [ K ] {δ } ,2=K∑ =1 &δ2T{ } [ M ] {δ&} ,причем матрицы жесткости [K] и масс [M] являются квадратными исимметричными.
Размер этих матриц равен числу степеней свободывсей системы, и вычисляются они путем соответствующей процедурысуммирования соответствующих матриц элементов, которые, в своюочередь, определяются интегральными выражениями:[ Kå ] = ∫ [ Bå ] T [ D ][ Bå ] dV ,(П.1)V[ M å ] = ρ ∫ [ N å ] T [ Nå ] dV .(П.2)VК свойству этих матриц можно добавить еще одно : при оптимальной нумерации узлов в системе эти матрицы приобретают ярковыраженный ленточный характер. Все эти свойства эффективно используются в вычислительных алгоритмах [159].Уравнение движения системы с распределенными по узловымточкам степенями свободы имеет вид :[ K ]{δ } + [C ]{δ&} + [ M ]{δ&&} = {F } ,где [C] – матрица демпфирования, а {F} – суммарный векторвсех действующих на систему сил.Для рассматриваемых в данной работе систем демпфирование неоказывает существенного влияния на частоту собственных колебаний,поэтому исключив внешние силы, и демпфирование мы получаемуравнение собственных колебаний:[ K ] { δ } + [ M ] { δ&&} = 0 .Это уравнение будет иметь вещественное периодическое решение{δ } = {δ 0} cos(ωt )266при выполнении условия[K ] − ω 2 [M ]= 0,представляющего собой классическую задачу на собственныезначения.
Для решения таких задач разработано большое количестворазличных алгоритмов [27, 92, 156, 160].При суперэлементном подходе, этот способ был предложен Айронсоном [10] и несколько позднее Гайяном [72], все степени свободы{δ } необходимо разделить на две части:δ1 ,δ 2{δ } = где перемещения{δ1 }.Переменные{δ 2 } однозначно выражаются через перемещения{δ1 } будем называть главными, а {δ 2 } – вспомога-тельными, тогдаI T {δ 2} = [T ]{δ1} и {δ } = T * {δ1} = {δ1} , где(П.3)[I ] – единичная матрица,[Т] – матрица связи главных и второстепенных перемещений.Уравнение равновесия системы можно представить в блочномвиде: K11K 21K12 δ1 F1 = .K 22 δ 2 F2 Если допустить предположение, что по второстепенным степеням свободы отсутствует нагрузка, в том числе и инерционная, то естьF2 =0, то из второго уравнения этой системы выразим {δ 2 }:{δ 2} = − [ K 22 ] −1 [ K12 ] T {δ1} .Таким образом, связь главных и второстепенных степеней свободы можно записать в виде:[T ] = − [ K 22 ] −1 [ K12 ] −T .267(П.4)Для случая собственных колебаний системы можно записать потенциальную и кинетическую энергию в виде:Ï =K=T111TTT{δ } [ K ]{δ } = {δ1} T * [ K ] T * {δ1} = {δ1} K * {δ1} ,2221 &δ2TTT{ } [ M ]{δ&} = 12 {δ&1} T * [ M ] T * {δ&1} = 12 {δ&1} M * {δ&1} ,Tа уравнение собственных колебаний без учета демпфирования:∂2[ Ê ]{δ1} + [ M ] 2 {δ1} = 0 ,∂t∗где∗] = Ò∗ Ò[ Ì ] Ò∗ ,[ Ê ∗ ] = Ò∗ Ò[ Ê ] Ò∗ .[Ì∗Таким образом, мы пришли к системе уравнений с меньшим ко личеством степеней свободы.
Решая характеристическое уравнение,записанное в виде: Ê ∗ + ω 2 M ∗ {δ1} = 0 ,получим значения собственных частот и формы соответствую щих им колебаний, выраженные через главные степени свободы.Для восстановления перемещений по второстепенным степенямсвободы логично было бы использовать соотношение П.3 с учетомП.4, однако соотношение П.4 было выведено в предположении, чтовторостепенные степени свободы не нагружены. Такой подход получил название статической конденсации. В реальности же к второстепенным степеням свободы должны быть приложены инерционныенагрузки, вызванные смещениями по главным степеням свободы.
Длянахождения уточненного решения записывается характеристическоеуравнение для полной системы в блочном виде: Ê11Ê 21Ê12 Ì−ω 2 Ê 22 Ì1121ÌÌ268 δ1 δ = {0} .22 2 12 После блочного перемножения второе уравнение разрешаетсяотносительно { δ2 }:{δ 2} = − [ Ê 22 ] − ω 2 [ M 22 ]−1[ Ê 21 ] − ω 2 [ M 21 ] {δ1} ,в результате чего получаем матрицу связи главных и второстепенныхстепеней свободы:−1[Ò] = − [ Ê 22 ] − ω 2 [ M 22 ] [ Ê 21 ] − ω 2 [ M 21 ].(П.5)Решение, полученное таким способом, естественно, будет зависеть от того , насколько удачно производится разделение степенейсвободы на главные и второстепенные, поэтому не может быть универсальным. Если амплитуды колебаний по второстепенным степенямсвободы минимальны, то влиянием инерционных сил, то есть слагаемыми, содержащими множитель ω2, можно пренебречь и полученноерешение является достаточно точным.
В противном случае полученное решение можно использовать в качестве начального приближенияи, получив с помощью соотношения 2.6 новые матрицы [К*] и [М*],итерационным путем уточнить частоту и форму колебаний на этой частоте.Приложение 2Функции формы используемых элементовДля записи функций форм вводится локальная система координат ξ-η- ς с центром в середине элемента. Координаты нормируютсятаким образом, чтобы диапазон их значений внутри элемента был впределах от –1 до 1.
Вводятся переменные ξ0 =ξξi, η0 =ηηi и ς0 =ςςi. Тогда функции формы можно записать следующим образом.Элемент 3D48:Ni =1(1 + ξ0 ) (1 + η0 ) (1 + ζ 0 ) (ξ0 + ζ 0 − 1)8269– для угловых узлов,Ni =11 − ξ 2 (1 + η0 )(1 + ζ 0 ) – для узлов на ребрах при ξ=0,4Ni =1(1 + ξ0 )(1 + η0 ) 1 + ζ 24()()– для узлов на ребрах при ζ=0.При этом в направлении оси 0η имеет место линейная аппроксимация перемещений, а в направлении других осей – квадратичная.Описание данного элемента и тестовые расчеты консольно закрепленных пластин приведены в [41].Элемент 3D72:Ni =1(1 + ξ0 )(1 + η0 )(1 + ζ 0 )[−10 + 9(ξ 2 + ζ 2 )] ï ðè ξ = ±1; η = ±1; ζ = ±1;64Ni =91(1 + ξ 0 )(1 + η0 )(1 − ζ 2 )(1 + 9ζ 0 ) ï ðè ξ = ±1; η = ±1; ζ = ± ;643Ni =91(1 + η0 )(1 + ζ 0 )(1 − ξ 2 )(1 + 9ξ 0 ) ï ðè ξ = ± ; η = ±1; ζ = ±1 .643В направлении осей 0ξ и 0ζ перемещения интерполируются полиномом третьей степени, а в направлении 0η – линейным.Элемент 3D96:для угловых узлов ( ξ = ±1; η = ±1; ζ = ±1 )Ni =1(1 + ξ0 ) (1 + η0 ) (1 + ζ 0 ) [9 (ξ 2 + η 2 + ζ 2 ) − 19] ,64для узлов на ребрахNi =91(1 + η0 )(1 + ζ 0 )(1 − ξ 2 )(1 + 9ξ 0 ) ï ðè ξ = ± ; η = ±1; ζ = ±1 ;643Ni =91(1 + ξ 0 )(1 + ζ 0 )(1 − η 2 )(1 + 9η0 ) ï ðè ξ = ±1; η = ± ; ζ = ±1 ;643Ni =91(1 + ξ 0 )(1 + ζ 0 )(1 − η 2 )(1 + 9η0 ) ï ðè ξ = ±1; η = ± ; ζ = ±1 .643270Приложение 3Исходный текст блока расчета матриц для элемента со смешанной линейно-кубической аппроксимацией перемещенийSUBROUTINE SAVE_72(NE,FLAG,CORD,NOP)LOGICAL FLAGCOMMON/NTAPE/NT1,NT2,NT3,NT4DIMENSION ESTIFM(72,72),ESMATM(72,72),CORD(1),NOP(1)REWIND NT3REWIND NT4DO 400 N=1,NEPRINT 106, NCALL IQT372(N,ESTIFM,ESMATM,FLAG,CORD,NOP)WRITE(NT3) ESTIFM,ESMATMWRITE(NT4) ESTIFM,ESMATM106 FORMAT(70X,I3)400 CONTINUERETURNENDSUBROUTINE IQT372(N,ESTIFM,ESMATM,FLAG,CORD,NOP)LOGICAL FLAGCOMMON/CONTR/NP,NE,NB,NDF,NCN,NSZF,CMY,XMAS,YMAS,RO,ZMAS,EMATCOMMON/NTAPE/NT1,NT2,NT3,NT4DIMENSIONESTIFM(72,72),ESMATM(72,72),XX(4),SKK(4),YY(2),BB(24),*CC(24),DD(24),ES(6,6),B(72,72),C(72,72),EST(72,72),E(72,72),*CORD(NP,NDF),NOP(NE,NCN)Cвесовые коэффициенты для численного интегрированияDATA SKK/.3478548,.6521451,.6521451,.3478548/,Cвесовые координаты для точек разбиения по оси у = 1Cкоординаты точек интегрирования1 XX/-.8611363,-.3399810,.3399810,.8611363/2 YY/-0.557735,0.557735 /AREAS=0.IX=NOP(N,1)DO 10 J=1,NCNIY=NOP(N,J)BB(J)=(CORD(IY,1)-CORD(IX,1))*XMASCC(J)=(CORD(IY,2)-CORD(IX,2))*YMAS10 DD(J)=(CORD(IY,3)-CORD(IX,3))*ZMASESTIFM=0.ESMATM=0.ES=0.E1=(1.-CMY)/((1.+CMY)*(1.-2.*CMY))ES(1,1)=E1ES(1,2)=E1*CMY/(1.-CMY)ES(1,3)=ES(1,2)ES(2,1)=ES(1,2)ES(2,2)=E1ES(2,3)=ES(1,2)ES(3,1)=ES(1,2)ES(3,2)=ES(1,2)ES(3,3)=E1G=E1*(1-2.*CMY)/(2.*(1-CMY))ES(4,4)=GES(5,5)=G271ES(6,6)=GESS=1.E-20DO 50 IX=1,4X=XX(IX)DO 50 IY=1,2Y=YY(IY)DO 50 IZ=1,4Z=XX(IZ)CALL ISOQ72(X,Y,Z,AREA,BB,CC,DD,EST,E)IF(AREA.GT.ESS) GOTO 230WRITE(NT2,23) AREA,N,X,Y,Z23 FORMAT(5X,' Det = ',E12.5,' Nelem =',I5,' COORD :',*' X=',F7.3,' Y=',F7.3,' Z=',F7.3)230 CONTINUECALL DDTED(C,EST,ES,72,6,72,72,6,B)A1=AREA*SKK(IX)*SKK(IZ)AREAS=AREAS+A1A2=A1*RO/981.DO 48 I=1,72DO 48 J=I,72ESTIFM(I,J)=ESTIFM(I,J)+C(I,J)*A1 * EMATESMATM(I,J)=ESMATM(I,J)+E(I,J)*A248 CONTINUE50 CONTINUEDO 65 I=1,71DO 65 J=I+1,72ESTIFM(J,I)=ESTIFM(I,J)ESMATM(J,I)=ESMATM(I,J)65 CONTINUEwrite(NT2,103)N,AREASprint103,N,AREASCPoljarwinckelnC Преобразование в полярную систему координатIF(FLAG) THENDO 1101 J=1,NCNI=NOP(N,J)1101DD(J)=ATAN(CORD(I,2)/CORD(I,3))110211031104DO 1103 I=1,72DO 1102 J=1,72B(I,J)=0.B(I,I)=1.DO 1104 J=1,NCNI=(J-1)*NDF +2B(I ,I )= COS(DD(J))B(I ,I+1)= SIN(DD(J))B(I+1,I )=-SIN(DD(J))B(I+1,I+1)= COS(DD(J))CONTINUECALL TMT(B,ESTIFM,72,72,C)CALL TMT(B,ESMATM,72,72,C)ENDIFCCRETURN220 WRITE(6,100)N100 FORMAT(5X ,' Nelem = ',I4,'ASRK4 Subr.