Диссертация (786428), страница 7
Текст из файла (страница 7)
Поэтому для вычисленияих конечных значений используется квадратурные формулы, при этом кинтегралам с сингулярными особенностями применяется каноническаярегуляризация [148].46§ 2.3. Алгоритм решения задачи на сверхзвуковом этапевзаимодействияВ основе алгоритма решения разрешающей системы уравнений (2.11),(2.15), (1.22) лежит модифицированный метод Рунге-Кутта четвертогопорядка и принцип редукции.Согласно принципу редукции заменяем ряды (2.8) конечнымисуммами, ограничиваясь в разложениях первыми членами с номерамиn = 0,1,..., N :w ( θ, τ )=p (θ, τ)Nwn ( τ )∑ p ( τ ) P ( cos θ ) ,nn =0n(2.20)u ( θ, ττ) N un ( ) dPn ( cos θ )=∑.dθc ( θ, τ ) n=0 c n ( τ )Приведем систему уравнений (2.11), (2.15) к нормальной форме Коши,т.е. приведем к системе первого порядка. Затем, используя представления(2.20),получаемдифференциальныхсистему,содержащуюуравнений,которая6 ( N + 1) + 2обыкновенныхдополняетсяалгебраическимуравнением (1.22).
Первые 6 ( N + 1) уравнений можно представить вматричной форме, имеющей блочную структуру: MW + Q,=W(2.21) ,W = W0 , W1 ,..., WN , , Wn = un , wn , χχn , un , wn , nTTun =un , w n =w n , χ n =χ n ,T=Q31Q=,Q,...,Q,Q0,0,0,0,pin ,0 ,∑01Nn 2hγγi =1T47M0M1.=M=, Mn.000=Mijn6× 6.MN1γ2−L11nL21nL31n−L12 nL22 nL32 n| γ2| 0|||−L13nL23n |L33n |00γ0−020−γ2−0К ним добавляется система уравнений движения оболочки какабсолютно твердого тела:uc = uc ,=ucp 3 N)I n (b( )),∑∑ pin (ττm0 =i 1 =n 0(2.22)уравнение для определения радиуса области контакта (1.22) и начальныеусловия: 0,0,0,0,0,0 T , n ≠ 1Wn ( 0 ) ==uc ( 0 ) V=, uc ( 0 ) 0,=0.0 , b ( 0)T 0,0,0,V0 ,V0 ,0 , n = 1(2.23)Пятое уравнение системы (2.21) в каждом блоке с номером n являетсяинтегро-дифференциальным, т.к.
его правая часть содержит слагаемые3∑pini =1, представленные в виде интегралов от неизвестных функций w n ( τ=) w n ( τ ) :N2n + 1p1n (τ) =−H (τττ)∑ w k ( ) I nk (b( ))2k =0tp2 n (t)p2n + 1 N∑ w k (t ) Pk [cos b(t )] dt ∫ υ[θ, b(t ), t − t ]Pn (cos θ)sin θd θ (2.24)2 k =0 ∫0048t2n + 1 Np3n (t) = −∑ w k (t )dt2 k =0 ∫0arcsin[b ( t )]∫0pdPk (cos θ* )d θ* ∫ υ(θ, θ* , t − t )Pn (cos θ)sin θd θ.d θ*0При этом в правую часть входят все N + 1 функции w 0n ( τ ) , n = 0, N , поэтомусистема решается совместно для всех 6 ( N + 1) + 3 уравнений.Отличие модифицированного метода Рунге-Кутты от классическогозаключается в том, что наряду с применением классической схемы [149]необходимо построение и использование квадратурных формул длявычисления интегралов в интегро-дифференциальных уравнениях системы(2.21).Построим дискретный аналог системы (2.21), (2.22), (1.22) с учетом(2.23), (2.24).
Численное решение задачи ищется в узлах равномерной сетки сшагом δm по временной координате τm =δm m . Искомые коэффициентырядов перемещений и их скоростей, контактного давления, радиус границыобласти взаимодействия, глубина погружения ударника как абсолютнотвердого тела и ее скорость заменяются их значениями в дискретныемоменты времени:=wn ( τm ) , χ nm =χτu0 n ( τm ) , w=u=un ( τm ) , w=w n ( τm ) ,nmn ( m), unmnmnm n ( m )=, pnmχ nm =χτ3∑ p ( τ ) , b = b ( τ ) , u=i =1inmmmcmuc ( τm ) .uc ( τm ) , u=cmОснованная на методе Кунге-Кутта разностная схема для системы(2.21), (2.22), (1.22) имеет следующий вид:Wm+1 = Wm +δm( R + 2K + 2Z + S ) ,6(2.25)49TWm+1 = W0,m+1 , W1,m+1 ,..., WN,m+1 ,T,Wn,m+1 =un ,m+1 , wn ,m+1 , χχn , m +1 , un ,m +1 , wn , m +1 , n , m +1Wm = W0m , W1m ,..., WNm ,T ,Wnm =unm , wnm , χχnm , unm , wnm , nmTR=R=, Rn0 , R 1 ,..., R NTK=Z=Tr1n , r2 n , r3n , r4 n , r5 n , r6 n ,K=, Kn0 , K 1 ,..., K NTZ=, Zn0 , Z1 ,..., Z NT=STk1n , k2 n , k3n , k4 n , k5 n , k6 n ,Tz1n , z2 n , z3n , z4 n , z5 n , z6 n ,S=, Sn0 , S1 ,..., S NTTs1n , s2 n , s3n , s4 n , s5 n , s6 n ,где1 1 R MWm + Q(R) , K= M Wm + R + Q(K) , Z= M Wm + K + Q(Z) ,=2 2 (R)(R)(K)(K), Q(K) = Q(K),S M ( Wm + Z ) + Q(S) , Q(R) = Q(R)=0 , Q1 ,..., Q N0 , Q1 ,..., Q NTTTT(S)(S)(Z)(Z),, Q(S) = Q(S)Q(Z) = Q(Z)0 , Q1 ,..., Q N0 , Q1 ,..., Q NTQ(R)nT3311(R)(K) (K) ),0 ,,w=p0,0,0,0,(),00,0,0,0,pinm (wQ=∑∑inmn022hγγhγγi =1i =1TQ(Z)n31 (Z)0,0,0,0,=pinm (w∑0 ),0 ,2hγγi =1Q(S)n31 (S)=0,0,0,0, ∑ pinm (w0 ),0 ,2hγγi =1TU c,m+1 = U cm +δm ∗R + 2K ∗ + 2Z∗ + S∗ ) ,(6T(2.26)TU c,m+1 = uc ,m+1 , uc ,m+1 , U cm = uc ,m , uc ,m ,∗∗1∗ T2R = r ,r∗∗1K = k ,k∗ T2∗1=, r2∗, r = ucmN 31 (R)Rp()I+pw∑∑einmnm ,m0 n 0=i 1=N 3r2∗ ∗ 1 (K) )I nm ,R e + p∑∑ pinm (w, k2, k= ucm + =2m0 =n 0=i 1∗150∗∗1Z = z ,z∗∗1∗ T2∗ T2S = s ,s∗1z ucm +,=k2∗ ∗ p N 3 (Z) )I nm ,pinm (w, z2 =∑∑2m0 =n 0=i 1p N 3 (S) )I nm ,pinm (w,=s ucm + z , s =∑∑m0 =n 0=i 1∗1∗2∗2I nm = I n (bm ) (см.(2.19));=bm+1 (R) = ( w ij )Здесь w(S)=w( wij+ z5i )uc ,m+1 ( 2 − uc ,m+1 ) .(2.27)r k (K)(Z)==wij + 5i , wwij + 5i ,,wN ×m2 N ×m2 N ×mN ×m.Заметим, что на сверхзвуковом этапе взаимодействия справедливоусловие=bm+1(1 − uc ,m +1) uc ,m +1uc ,m+1 ( 2 − uc ,m+1 )≥ 1.Квадратурные формулы для интегралов, входящих в выражения для), p2 n ( ), p3n ( ) (см.(2.24)), строятся с использованиемопределения p1n (τττметода весовых коэффициентов и канонической регуляризации [148] дляполучения конечных значений сингулярных интегралов.Сеточные функции p1nm , p2nm , p3nm с учетом (2.24) запишутся так:p1nm2n + 1 N≈−∑ w km I nkm , I nkm = I nk (bm ) (см.(2.18)),2 k =051p2 nm ≈l2n + 1 N m−1(cos)wPb∑∑ ki ki ∑ Pn (cos θ j )sin θ j ×2 =k 0=i 0=j 1× δm {w1 j (bi ) f1 (θ j , bi , tm − ti ) + w3 j (bi ) f 2 (θ j , bi , tm − ti )} +(2.28) i {w1 j (bi ) f3 (θ j , bi , tm − ti ) + δl f 4 (θ j , bi , tm − ti )} ,+wp3nml2n + 1 N m−1 m1≈−Pn (cos θ j )sin θ j ×∑∑ w ∑[ Pk (cos θ*i1 ) − Pk (cos θ*i1−1 )]=∑2 =k 0=i 0 ki =i1 1j 1× δm {w1 j (θ*i1 ) f1 (θ j , θ*i1 , tm − ti ) + w3 j (θ*i1 ) f 2 (θ j , θ*i1 , tm − ti )} +} i {w1 j (θ*i1 ) f3 (θ j , θ*i1 , tm − ti ) + δl f 4 (θ j , θ*i1 , tm − ti ) ) ,+wπbгде δl = , δm1 = i , θ j = jδl , ti = iδm , θ*i1 = i1δ m1 ;lm1f1 (q j , x, tm − ti ) =f 2 (q j , x, tm − ti ) =f3 (q j , x, tm − ti ) =2∑ϑq =1rq1(q j , x, tm − ti ) ,2∑ϑq =1rq 2(q j , x, tm − ti ) ,rq 3(q j , x, tm − ti ) ,2∑ϑq =1f 4 (θ j , x, tm − ti ) =ϑs (θ j , x, tm − ti ) .Весовые коэффициенты имеют вид: tm − ti, t ≠tdt ln i ∫ = tm − ti +1 i +1 m=ωtm − t tiln tm − tm−1 , ti +1 =tmti +1 F ( x, θ j ) − F ( x, θ j −1 ), x ≠ sin θ j , x ≠ sin θ j −1dθωlj ( x ) =∫= F ( x, θ j ), x =sin θ j −1( l =1,3) ,lsinxθ−) θ j −1 (− F ( x, θ j −1 ), x =sin θ jθj52где F ( x, z ) - первообразная подынтегральной функции.Значения полных эллиптических интегралов 1-го и 2-го родавычисляются с помощью аппроксимации многочленами [150]:K (m ( x)) = a0 + a1m 1 ( x) + ...
+ a4 m 14 ( x) ++ b0 + b1m 1 ( x) + ... + b4 m 14 ( x) ln(1 m 1 ( x)) + ε(m ( x)),ε(m ( x)) ≤ 2 ⋅ 10−8 ,a0 1.38629436112,a1 0.09666344259,a2 0.03590092383,===a3 0.03742563713,a4 0.01451196212,==b0 0.5,b1 0.12498593597,b2 0.06880248576,===b3 0.03328355346,b4 0.00441787012;==E (m ( x)) = 1 + a1m 1 ( x) + ... + a4 m 14 ( x) ++ b1m 1 ( x) + ... + b4 m 14 ( x) ln(1 m 1 ( x)) + ε(m ( x)),ε(m ( x)) ≤ 2 ⋅ 10−8 ,=a1 0.44325141463,=a2 0.06260601220,=a3 0.04757383546,=a4 0.01736506451,=b1 0.24998368310,=b2 0.09200180037,=b3 0.04069697526,=b4 0.00526449639,где m ( x) =4θ j x(θj+ x)20 ≤ m ( x) < 1, m 1 ( x) =1 − m ( x) .Неполные эллиптические интегралы 1-го и 2-го рода вычисляются спомощью метода Симпсона [149] ( δ эл. , n – шаг в эллиптическом интеграле ичисло шагов соответственно):531 δα=Fxmxy(),(),где()q 1 − m ( x)sin 2 α ,I = E δ ( x), m ( x) , где y(α)= 1 − m ( x)sin 2 α qnn −1δ эл.
I≈y0 + yn + 4∑ y2i −1 + 2∑ y2i ,3 =i 1 =i 1 δ q ( x)δ эл. =, y0 = y (0), yn = y δ q ( x) , y2i −1 = y [ (2i − 1)δ эл. ] ,2n (q + x )η t 2 − (q − x ) 2 qkj , q =1,2.αrcsin j=y2i y (2iδ эл. ), δq ( x) = 2(tm − ti )qjxК системе уравнений (2.25) – (2.28) добавляются начальные условия:T 0,0,0,0,0,0 , n ≠ 1Wn 0 ==uc ,0 V=, uc ,0 0,=0.0 , b0T=VVn0,0,0,,,0,100(2.29)Предложенный алгоритм расчета реализован в среде Delphi [151].54§ 2.4. Примеры расчетовВ качестве первого примера рассмотрим вертикальный удар стальнойсферическойоболочкипополупространству,заполненномусталью.Соответствующие безразмерные параметры имеют следующие значения:R 1, =h 0.05,=V 0.1.=γ 1, =На рис.
2.1. и 2.2. представлены зависимости положения границыобласти контакта и ее скорости расширения в зависимости от времени.Рис. 2.1.55Рис. 2.2.Зависимости контактного давления от времени в разных точкахобласти контакта представлены на рис. 2.3. Сплошная кривая соответствуетлобовой точке оболочки, штриховая – точке с координатой θ =0.04 ,штриховая пунктирная –θ =0.08 .
Положения скачков на графикахсоответствуют моментам времени, при которых радиус расширяющейсяобласти контакта b ( τ ) достигает соответствующей точки.56Рис. 2.3.На рис. 2.4. представлены распределения контактного давления поугловой координате θ в различные моменты времени. Сплошная криваясоответствует моменту времени τ =0.01 , штриховая – τ =0.03 , штриховаяпунктирная – τ =0.06 . Скачки на графиках соответствуют положениюграницы области контакта в соответствующий момент времени. Как видно,на подвижной границе области взаимодействия контактное давление имеетразрыв первого рода.