С.Ю. Никитин, С.С. Чесноков - Механика (1111872), страница 29
Текст из файла (страница 29)
Численный анализ в механике239§16. Численный анализ в механикеКраткие теоретические сведенияЗадача о движении системы из N взаимодействующих материальных точек описывается N векторными уравнениями, вытекающими из законов Ньютона (см. § 4):mi&&ri =∑ f ij + F ,j ≠iii = 1, 2, ⋅⋅⋅, N ,(16.1)с начальными условиямиri (0) = ri0 , &ri (o) = vi0 , i = 1, 2, ⋅⋅⋅, N .(16.2)Здесь mi – масса i–й точки, ri – ее радиус-вектор, fij – сила взаимодействиямежду i–й и j–й точками, входящими в систему (внутренняя сила), Fi – сумма сил, действующих на i–ю точку со стороны тел, не входящих в систему(внешняя сила). Задача, описываемая уравнениями (16.1) с начальными условиями (16.2), представляет собой хорошо известную в математике задачуКоши для системы обыкновенных дифференциальных уравнений второгопорядка.
В общем виде такая задача при N > 2 решения в квадратурах неимеет. Поэтому в задачах механики многих тел широко применяются методы численного анализа с использованием современной вычислительнойтехники.Проиллюстрируем вначале применение численных методов напримере одного уравнения динамикиmd2rdt 2dr⎛ dr ⎞= F⎜ r, , t ⎟ , r(0) = r0 ,(0) = v0 ,⎝ dt ⎠dt(16.3)а затем обобщим на случай произвольного числа уравнений.
При постановке задачи для численного решения обычно осуществляется переход от ре-§16. Численный анализ в механике240альных переменных, имеющих конкретную размерность (метр, секунда,килограмм и т.п.) к так называемым безразмерным переменным, которыеобразуются из соответствующих размерных переменных путем деления ихна заранее выбранные масштабы. Необходимость такого перехода вызвана,во-первых, ограниченностью диапазона представления чисел в памяти компьютера и неизбежной потерей точности при выполнении процессоромарифметических операций над числами, отличающимися друг от друга намного порядков, а во-вторых, возможностью выявления весьма полезныхсоотношений между параметрами подобия задачи, существенно способствующих ее более глубокому пониманию.
Вопрос о конкретном выборемасштабов для переменных решается в каждой задаче отдельно.Рассмотрим приведение к безразмерному виду задачи (16.3). Пустьхарактерная область изменения радиус-вектора частицы имеет размер L, ахарактерное время задачи равно T. Введем безразмерные радиус-вектор ρ ивремя τ в соответствии с выражениямиρ=rt, τ= .LT(16.4)При этом уравнение (16.3) преобразуется к видуmL d 2ρT2L dρ⎛⎞, T τ⎟ .= F ⎜ L ρ,⎝⎠τTddτ2(16.5)Вводя безразмерную силу Φ по формулеL dρ⎞⎛ dρ ⎞ T 2 ⎛Φ ⎜ ρ, , τ⎟ =F⎜ L ρ,, T τ⎟ ,⎠⎝ dτ ⎠ L m ⎝T dτ(16.6)получаем окончательноd 2ρdτ 2dρ⎛ dρ ⎞= Φ ⎜ ρ, , τ⎟ , ρ(0) = ρ0 ,(0) = ϕ 0 ,⎝ dτ ⎠dτ(16.7)где ρ0 = r0/L, ϕ0 = v0T/L – начальные условия для безразмерных перемен-§16. Численный анализ в механике241ных. Векторная задача (16.7) распадается на три скалярных задачи в проекциях на декартовы оси координат:d 2ξ⎛ dρ ⎞= Φ ξ ⎜ ρ, , τ⎟ , ξ(0) = ξ 0 ,⎝ dτ ⎠dτdξ(0) = ϕ ξ0 ,dτd2ηdη(0) = ϕ η0 ,dτ2⎛ dρ ⎞= Φ η ⎜ ρ, , τ⎟ , η(0) = η0 ,⎝ dτ ⎠dτ2d 2ς⎛ dρ ⎞= Φ ς ⎜ ρ, , τ⎟ , ς(0) = ς 0 ,⎝ dτ ⎠dτ2(16.8)dς(0) = ϕ ς0 ,dτгде {ξ, η, ς} = ρ .
Рассмотрим в качестве примера одно скалярное уравнение:d 2ξ⎛ dξ ⎞= Φ⎜ ξ, , τ⎟ , ξ(0) = ξ 0 ,⎝ dτ ⎠dτ2dξ(0) = ϕ 0 .dτ(16.9)Запишем дифференциальное уравнение второго порядка (16.9) каксистему двух уравнений первого порядкаdξ= ϕ,dτdϕ= Φ( ξ, ϕ, τ),dτξ(0) = ξ 0 ;(16.10)ϕ(0) = ϕ 0 .Основная идея численного интегрирования систем обыкновенных дифференциальных уравнений вида (16.10) состоит в дискретизации переменнойτ, т.е. разбиении оси Oτ на отрезки длиной Δτ и введении сетки τ i = iΔτ ,i = 0, 1, 2, ⋅ ⋅⋅ . Особенности, связанные с такой дискретизацией, удобно рассмотреть вначале на примере одного уравнения первого порядка, а затемобобщить на случай системы уравнений.
Пустьdu= f (u, τ) , u(0) = u0 .dτПроинтегрируем (16.11) по интервалу Δτ между точками τi и τ i+1 :(16.11)§16. Численный анализ в механике242τ i+1ui +1 = ui +∫ f (u, τ)dτ .(16.12)τiРазные численные методы отличаются друг от друга подходами к аппроксимации интеграла в правой части (16.12).Простейший способ аппроксимации состоит в замене подинтегральной функции ее значением в точке ui , τ i . В результате получаем алгоритм, называемый методом Эйлера:ui +1 = ui + f (ui , τ i ) Δτ , i = 0, 1, 2, ⋅⋅⋅(16.13)Расчетная схема (16.13) проста и экономична. Однако надо иметь в виду,что она не всегда устойчива по отношению к малым ошибкам, допущеннымна каком-либо шаге вычислений.
Исследования показывают, что метод Эйлера устойчив, если выполнены следующие условия:∂f2.≤ 0 , Δτ ≤∂u| ∂f / ∂u|(16.14)В противном случае ошибки вычислений на отдельных шагах быстро накапливаются и могут сильно исказить решение, сделав его непригодным дляпрактических целей.
Поскольку в задачах механики условие ∂f / ∂u ≤ 0 выполняется не всегда, применимость метода Эйлера в этих задачах ограничена.Одним из наиболее употребительных на практике является методРунге-Кутта четвертого порядка. В этом методе на каждом шаге вычисляются четыре числа:k1( i) = f (ui , τ i ) ⋅ Δτ,k 2( i) = f (ui + k1( i) / 2, τ i + Δτ / 2) ⋅ Δτ,k 3( i) = f (ui + k 2( i) / 2, τ i + Δτ / 2) ⋅ Δτ,k 4( i) = f (ui + k 3( i) , τ i + Δτ) ⋅ Δτ.(16.15)§16.
Численный анализ в механике243Последовательные значения ui+1 искомой функции определяются так:ui +1 = ui + Δui , Δui =()1 ( i)k + 2k 2( i) + 2k 3( i) + k 4( i) .6 1(16.16)Аналитические оценки устойчивости метода Рунге-Кутта в общем случаепроизвести не удается. Поэтому для определения правильности выбора шага Δτ обычно применяют на каждом этапе расчета, состоящем из двух шагов, двойной пересчет: с выбранным (Δτ) и половинным (Δτ/2) шагами. Если расхождение полученных при этом решений не превышает наперед заданной погрешности, то оставляют выбранный шаг; в противном случаешаг уменьшают вдвое.Метод Рунге-Кутта легко обощается на случай системы уравнений.Например, для двух уравненийdu= f (u, v, τ),dτdv= g(u, v, τ)dτ(16.17)расчеты ведутся по формулам:k1( i) = f (ui , v i , τ i ) ⋅ Δτ,m1( i) = g(ui , v i , τ i ) ⋅ Δτ,k 2( i) = f (ui + k1( i) / 2, v i + m1( i) / 2, τ i + Δτ / 2) ⋅ Δτ,m2( i) = g(ui + k1( i) / 2, v i + m1( i) / 2, τ i + Δτ / 2) ⋅ Δτ,⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅k 4( i) = f (ui + k 3( i) , v i + m3( i) , τ i + Δτ) ⋅ Δτ,m4( i) = g(ui + k 3( i) , v i + m3( i) , τ i + Δτ) ⋅ Δτ,(()1 ( i)k + 2k 2( i) + 2k 3( i) + k 4( i) ,6 11Δv i = m1( i) + 2m2( i) + 2m3( i) + m4( i) ,6ui +1 = ui + Δui , v i +1 = v i + Δv i .Δui =)(16.18)§16.
Численный анализ в механике244Примеры численного анализаПример 16.1. Ракета ведет преследование цели так, что вектор ее скоростиv в каждый момент времени направлен на цель, а модуль скорости v является постоянным. Цель движется по параболической траектории как тело,брошенное со скоростью v0 под углом α к горизонту.
Рассмотреть несколько точек запуска ракеты и провести для них расчет семейства траекторийракеты, запускаемой в разные моменты времени. Для каждого семействавыбрать траекторию с минимальной кривизной, двигаясь по которой ракетауспеет поразить цель над поверхностью Земли.Анализ. По условию задачи движение ракеты и цели происходит в однойплоскости, перпендикулярной поверхности Земли. Совместим с этой плоскостью координатнуюплоскость XOY, начало координат поместим встартовую точку цели. Обозначим через x, y координаты ракеты, через x1,y1 – координаты цели.За начало отсчета времени примем момент пускаРис. 16.1цели, момент запуска ракеты обозначим через ts.В соответствии со стратегией преследования касательная к траектории ракеты в каждый момент времени проходит через точку C, отображающую местонахождение цели (рис.
16.1). Следовательно,dy y1 − yy& y − y, или = 1.=dx x1 − xx& x1 − x(16.19)Воспользовавшись соотношением x& 2 + y& 2 = v 2 , выразим из (16.19) проекции скорости ракеты:v( x1 − x )dx=, t ≥ t s,dt( x1 − x ) 2 + ( y1 − y) 2v( y1 − y)dy=, t ≥ t s.dt( x − x ) 2 + ( y − y) 211(16.20)§16. Численный анализ в механике245Обозначив через x0, y0 координаты точки запуска ракеты, сформулируемначальные условия для ее движения:x (t s ) = x 0 , y(t s ) = y 0 .(16.21)Движение цели по условию задачи описывается кинематическимисоотношениямиx 1 (t ) = v0 cos α ⋅ t,y1 (t ) = v0 sin α ⋅ t −gt 2.2(16.22)В качестве масштабов длины и времени выберем дальность полета цели L ивремя ее полета от начальной до конечной точки T:L=2v022vsin α cos α, T = 0 sin α .gg(16.23)Введем безразмерные координаты и время по формулам:ξ=xyxyt, η = , ξ1 = 1 , η1 = 1 , τ = .LLLLT(16.24)В этих переменных соотношения (16.22), описывающие движение цели,принимают вид:()ξ1 = τ, η1 = tg α τ − τ 2 .(16.25)Переход к безразмерным переменным в системе уравнений (16.20) дает:v( ξ1 − ξ)1 dξ=, τ ≥ τ s , ξ( τ s ) = ξ 0 ,L dτ( ξ1 − ξ) 2 + ( η1 − η) 2v( η1 − η)1 dη, τ ≥ τ s , η( τ s ) = η0 ,=L dτ(ξ − ξ) 2 + ( η − η) 211(16.26)§16.
Численный анализ в механике246где τ s = t s / T , ξ 0 = x 0 / L , η0 = y0 / L . Подставляя в (16.26) координатыцели как явные функции времени, окончательно получаем:dξ=dτdη=dτu( ξ1 − ξ)2( (( τ − ξ) + tg α ⋅ τ − τ2) − η)2( ( ) )( τ − ξ) + (tg α ⋅ ( τ − τ ) − η)(16.27)u tg α ⋅ τ − τ 2 − η22, τ ≥ τ s , ξ( τ s ) = ξ 0 ,2, τ ≥ τ s , η( τ s ) = η0 ,где u = v / (v0 cos α) – безразмерная скорость ракеты.Для интегрирования системы уравнений (16.27) воспользуемся методом Рунге-Кутта. Шаг интегрирования Δτ здесь и далее будем подбиратьтаким образом, чтобы его изменение вдвое вызывало относительное изменение решения в конечной точке не более, чем на δ = 10 −5 . Кривизну траектории K будем вычислять по формуле22⎛ d 2ξ ⎞⎛ d2η⎞d 2ρ1K = = 2 = ⎜⎜ 2 ⎟⎟ + ⎜⎜ 2 ⎟⎟ ,Rds⎝ ds ⎠⎝ ds ⎠(16.28)где R – радиус кривизны, ds = (dξ) 2 + (dη) 2 – дифференциал длины дугитраектории.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
