Текст книги "Введение в анализ риска"

Рассмотрим приближенный метод построения функции w(·) для случайного процесса x_i при вычислении показателей риска.

Предположим, что x_i – стационарный случайный процесс. Для вычисления w(t;x_i) строится некоторый вспомогательный векторный гауссовский марковский процесс y(t), размерность которого не больше двух.

Решение задачи разобьем на два этапа: построение одномерной плотности распределения w(t;x_i) и корреляционной функции для процесса x_i(t); построение процесса y(t) = (у₁,…,у_n) с использованием w(t;x_i), .

Наиболее подходящим аналитическим выражением для указанных плотностей являются достаточно компактные выражения плотности S-распределений Джонсона, допускающие автоматизацию процесса их выбора [25, 26].

Класс плотностей Джонсона состоит из следующих трех семейств:

– семейства плотностей S_L-распределения Джонсона

где Z > ε; δ > 0; –∞ < γ^* < ∞; γ^* = γ – δLnλ; –∞ < ε < ∞;

– семейства плотностей S_в-распределения Джонсона

где δ > 0; –∞ < γ < ∞; λ > 0; –∞ < ε < ∞; ε < Z < ε + λ;

– семейства плотностей S_U-распределения Джонсона

где –∞ < γ < ∞, δ > 0; –∞ < γ < ∞; λ > 0; –∞ < ε < ∞. В выражениях (5.15)–(5.17) ε, δ, γ, λ – параметры плотностей; Z обозначает y_i в момент времени t_k.

Принадлежность y_i(x) к одному из семейств плотностей Джонсона выбирается на основе анализа начальных или центральных моментов первых четырех порядков. Вычисление указанных моментов осуществим с помощью модифицированного моментно-семивариантного метода [41]. Корреляционную функцию для x_i(t) найдем с помощью аппарата кумулянтного анализа. При этом получим систему дифференциальных уравнений относительно некоторой совокупности кумулянтных функций как исследуемой i-й, так и остальных компонент векторного процесса X(t). Одна из кумулянтных функций приближенно описывает искомую корреляционную функцию, которая представляется в виде таблицы.

В случае когда система (5.7) допускает линеаризацию, искомую корреляционную функцию K(τ) получим, решив матричные дифференциальные уравнения

где K(τ) = (k_ij); k_ij – элементы корреляционной матрицы K размера (n х n); A = (a_ij) – элементы матрицы A коэффициентов системы возмущенного движения ЛА; C = (c_ij) – матрица коэффициентов уравнения Колмогорова.

Непосредственное решение приведенного уравнения найти довольно сложно. Поэтому рассмотрим иной подход к его решению. Предполагая процесс x(t) стационарным, перейдем от уравнения (5.18) к уравнению KA + BK = –C, а затем к эквивалентному матричному уравнению

NK_o = α,           (5.19)

где В = А^Т;

К_o = (k₁₁,…,k_1n,k₂₁,…,k_n1,…,k_nn)^T; α = (c₁₁,…,c_1n,c₂₁,…,c_n1,…,c_nn)^T– векторы-столбцы. Используя кронекеровские суммы и произведения [6], исходную систему представим в виде

N = A x E + E x B^T = B x E + E x A^T,

где A х E, E х B^T – кронекеровские произведения матриц

Решив (5.19), получим искомую корреляционную функцию для стационарных уравнений после окончания переходного процесса, т. е. при t → ∞.

Если порядок матриц A и В мал, то уравнения (5.19) решаются известными аналитическими методами, как, например, К_о = N^–1α, где N^–1 – обратная матрица для N.

На втором этапе построим вспомогательный марковский процесс, используя полученный на первом этапе вид плотности ^-распределения Джонсона для y_i. Если имеем S_L-распределение Джонсона (5.15), то

если плотность y_i представлена в виде плотности S_в-распределения Джонсона (5.16), то

и, наконец, если плотность y_i представлена в виде плотности S_U-распределения Джонсона (5.17), то

В последних выражениях  – некоторый нормированный гауссовский процесс; λ, γ, δ, ε – параметры плотности распределения Джонсона.

Таким образом, имеем

где 1 < j < 3 в зависимости от конкретного вида плотности S-распределения Джонсона.

Подберем для процесса корреляционную функцию ρ_y(τ) таким образом, чтобы корреляционная функция процесса y_i(t), формируемого в соответствии с (5.20), совпадала с корреляционной функцией, найденной на предыдущем этапе. Такой подбор, как следует из работы [41], удобно осуществить в соответствии с формулой

где  – корреляционная функция процесса y_i(t), значения которой найдены на предыдущем этапе;

H_n – полином Эрмита n-й степени.

Для вычисления неизвестной функции воспользуемся методами степенных рядов и вычислим , где ρⁱ, ρ⁽ⁱ⁺¹⁾ – функция ρ_y(τ), полученная из (5.21) в случае, когда ряд содержит i и (i+1) членов. В качестве погрешности используем . Если эта величина меньше заданной малой величины ε, то принимается в качестве расчетной. При этом k есть то значение индекса, при котором вычислительный процесс прекращается. Полученные числовые значения аппроксимируем одной из стандартных корреляционных функций, например [62]:

Последней функции соответствует двумерный марковский процесс

где Y = (Y₁,Y₂)^T; A = (a_ij)_2×2; B = (b₁,b₂)^T; n – белый шум с нулевым средним и ковариационной функцией cov[n(t₁), n(t₁)] = δ(t₁ – t₂); δ(τ) – дельта-функция; a₁₁=0; a₁₂=1; a₂₁= –(b²+ω²); a₂₂= –2b; b₁ = [2(b – βω)]^1/2; b₂ = [2(b + βω)(b + ω)]^1/2 – 2b[2(b – βω)]^1/2.

5.3. Параметры, допускающие выбросы в критическую область

Для ряда параметров движения, например угла атаки или скорости полета, допустимы кратковременные выбросы в критическую область. В этом случае вероятности P_пр и P_лс определяются из (2.3) и для их вычисления необходимо знать условную плотность вероятности W(θ/x_i) времени θ выброса параметра x_i и допустимое θ₀ время выброса в критическую область ограничиваемого параметра x_i.

Пусть самолет совершает полет с углом атаки α₀ < α_СВ, где α₀, α_СВ – соответственно угол атаки в исходном режиме и угол атаки, при котором происходит сваливание самолета на крыло. Определим величину допустимой длительности θ₀ пребывания текущего значения угла атаки α = α₀ + Δα в области критических значений α < α_СВ, при которой самолет за счет инерционности изменения аэродинамических сил и моментов не входит в режим сваливания. Все дело в том, что не угол α определяет критический режим, а поле сил аэродинамического давления, время его образования плюс инерция всей системы.

Процесс сваливания ЛА протекает следующим образом. Внешнее возмущение, вызванное приращением угла атаки Δα в начальный момент времени, на каждой половине крыла приводит к изменению циркуляции Γ. В результате происходит перераспределение давления на крыле, что приводит к изменению подъемной силы Y и моментов M_x, M_y, M_z. При достижении углом атаки критического значения происходит срыв потока на левой и правой полуплоскостях крыла различным образом. Образующиеся при этом угловые скорости и ускорения приводят в соответствующие моменты времени к изменению начальных углов и соответствующему изменению на величину ΔΓ_i циркуляции Γ. С учетом сказанного и результатов работы [67], построим уравнения, описывающие динамику сваливания.

Изменение циркуляции сечения крыла на больших углах атаки в нестационарных условиях описывается интегродифференциальным уравнением

где C_упр = C_упр(α) – коэффициент подъемной силы прямого крыла; l/2 — полуразмах крыла;

Здесь X₁, Z – координаты точек, в которых вычисляется индуктивная скорость;

θ_y – нормальная составляющая скорости набегающего потока; = 2X / b; X₀(Z) – абсцисса средней линии хорд; – область, занятая вихревым следом. Величина θ_y записывается так:

Подставив θ_y в (5.25), после интегрирования получим

где θ_∞ = V_L(1 + βtgχ); V_L – путевая скорость невозмущенного движения.

Вблизи критического угла атаки α_кр изменение функции C_у = C(α) зависит от направления изменения угла атаки α. При этом, когда α уменьшается, график функции C_у = C_у(α) располагается несколько ниже, чем при увеличении α, т. е. имеет место явление гистерезиса. Поэтому при вычислении подъемной силы Y и моментов M_x, M_y, M_z вблизи критического угла атаки необходимо учитывать как нелинейность зависимости C_у = C_у(α), так и явление гистерезиса. При этих условиях возмущенное движение ЛА описывается системой нелинейных дифференциальных уравнений

q = ρV²_∞ / 2, S – площадь.

Величины C_y, M_x, M_y, m_z, входящие в (5.28), определяются по формулам

В этих уравнениях учиты ваются нестационарные эффекты образования аэродинамических сил, скос потока, нелинейные эффекты обтекания, влияние инерционных и демпфирующих сил.

Для вычисления допустимой продолжительности θ₀ выброса в опасную область возмущенного угла атаки зададим ступенчатую функцию α = α₀ + α_w, где α_w – ступенчатая функция возмущения, амплитуда которой такова, что α > α_CB. Процесс интегрирования системы (5.22)–(5.29) по времени продолжается до тех пор, пока α не станет меньше α_св. Продолжительность этого процесса даст продолжительность пребывания α в закритической области, т. е. длительность выброса θ₀.

Расчет величины θ₀ проведен на примере ЛА МИГ-25 при следующих начальных условиях: α₀ = 18^o, V₀ = 135 м/с. Зависимость коэффициента подъемной силы от угла атаки зададим в виде

где K – коэффициент, принимающий одно из значений: 1.338; 1.115; 1.561, с помощью которого моделируется гистерезис зависимости C_y(α); A, B – постоянные коэффициенты.

При разработке алгоритма расчета θ₀ коэффициенты уравнений (5.22)–(5.28) определяются по формулам (5.29), что позволяет рассчитывать приращение параметров в каждый последующий момент времени. Решение интегродифференциальных уравнений (5.22)–(5.28) проводится по методу, изложенному в работе [67], путем сведения их к системе алгебраических уравнений. Порядок расчета следующий.

1. Задается время t₀ начала интегрирования.

2. Задаются исходные данные (для конкретного ЛА) и начальные условия: k – число сечений крыла, в которых вычисляется Г; b(z) – хорда крыла; L – размах крыльев; S_кр – площадь крыла; b_A – средняя аэродинамическая хорда; C^φ_y, , m^φ_z,  – аэродинамические коэффициенты; J_x, J_y, J_z – моменты инерции; m – масса объекта; α₀ – начальный угол атаки; C_y(α₀) – коэффициент подъемной силы при α = α₀; M_y(α₀); M_z(α₀) – моменты относительно осей OY, OZ соответственно; q₀ – скоростной напор;  – коэффициенты.

3. Определяются Γ⁰_i, .

4. Для полученных Γ⁰_i при α = α₀ определяются C_y, M_x, M_y.

5. Вычисляются .

6. Приближенным методом решается система (5.28) при соответствующих начальных условиях.

7. Полученные значения α, ω_z, β в момент времени t₀ + nΔt используются для определения Γ⁰_i.

8. Пункты 3–7 повторяются для моментов времени t₀ + nΔt до тех пор, пока не достигается α = α₀.

9. Если α → ∞, то θ′_w = θ⁰_w + Δθ.

10. Максимальное значение θ_w, при котором , и дает θ₀.

Автором проведена оценка влияния величины Δα в (5.22), характеризующей влияние скоса потока на крыле от вихревого следа, на изменение циркуляции. При этом предполагалось, что начальные значения угловых скоростей ω_x, ω_y, ω_z равны нулю. Величины Γ_ν, где v – номера рассматриваемого сечения крыла, вычисленные с учетом и без учета вихревого следа, отличаются для различных сечений v по размаху крыла не более чем на 0.1 %. Поэтому в расчетах третьим членом в (5.22) пренебрегали, т. е. исключали его из системы. Заметим, что при больших начальных угловых скоростях ω_x, ω_y, ω_z значение углов скоса потока на крыле от вихревого следа может быть существенным. Тогда расчетная схема должна включать полную модель.

При интегрировании системы за счет Δα происходит сначала рост подъемной силы на отрезке времени θ₀, а затем ее уменьшение. Значение θ₀, начиная с которого система не возвращается в исходный режим полета, т. е. в область α < α_CB, определяется для различных значений k.

В результате расчетов получено:

– при увеличении коэффициента подъемной силы C_y(α), когда k =1.561, и его уменьшении, когда k =1.338, время θ₀ = 5.4 c;

– при увеличении коэффициента подъемной силы C_y(α), когда k =1.338, и его уменьшении, когда k =1.115, время θ₀ = 6.2 c;

– если увеличение и уменьшение коэффициента подъемной силы C_y(α) происходят при постоянном k =1.115, то время θ₀ = 7,3 с.

Таким образом, уменьшение коэффициента k в функции f(α) приводит к увеличению времени θ₀ пребывания угла атаки в области критических значений без потери устойчивости ЛА. При этом уменьшение k приводит к уменьшению C_y1 = f (α_кр) и уменьшает крутизну графика f(α) в области α_кр = α_св.

Отметим, что этот раздел предназначен для специалистов в области нестационарной аэродинамики.

Разработанная ниже методика расчета позволяет определить допустимую продолжительность θ₀ пребывания параметров движения в критической области по заданным конструктивно-аэродинамическим параметрам ЛА.

При некоторых допущениях методика расчета величины θ₀, полученная для перегрузки n_y и угла атаки α, применима для таких параметров состояния ЛА, как число Маха M, минимальная и максимальная приборные скорости полета, координата X_T положения центра масс самолета.

Пусть текущее значение числа М превышает его критическое значение M_кр на величину ΔM на отрезке времени θ. При этом происходит потеря устойчивости. Это обусловливает практически неконтролируемые летчиком движения, при которых превышаются или максимально допустимые по прочности конструкции перегрузки (n_ymax, n_zmax), или предельно допустимые по сваливанию в штопор углы атаки α_CB [47, 58]. Так как в большинстве случаев достижение n_ymax или α_CB случайно, то при расчетах системы предупреждения критических режимов следует принять значение θ₀, полученное для перегрузки n_y.

Рассмотрим действие силы F, которая внезапно прикладывается к крылу самолета в момент времени t = 0 и остается постоянной в течение некоторого отрезки времени Δt, а затем также внезапно исчезает. При этом даже весьма значительная нагрузка может оказаться безопасной, если длительность Δt мала по сравнению с периодом свободных колебаний крыла. За время Δt крыло ЛА испытывает дополнительные относительные удлинения ΔE₁ и ΔE₂ от изгиба и кручения соответственно. При высоких скоростях упругой деформации крыла допустимая величина относительного удлинения статического нагружения E_доп может быть увеличена на ΔE ≤ kE_доп [76]. В работе [76] предложено выбирать k = 0.1. Будем иметь в виду, что значение k взято по нижнему пределу. В результате переходим к следующей задаче.

Пусть самолет совершает установившееся движение при n_y = n_y доп, при котором относительное удлинение E = E_доп. Начиная с момента t = 0, на крыло действует дополнительная сила Δn_y, которая приводит к приращению относительного удлинения на величину ΔE. Требуется определить время θ₀ действия силы Δn_y, в течение которого неравенство ΔE ≤ kE не нарушается.

Установим связь между относительной деформацией E и упругими деформациями крыла большого удлинения. С этой целью деформации изгиба Y(Z,t) и угол кручения φ(Z,t) представим в виде [45]:

где {f_i(Z)}, {φ_i(Z)} – система линейно независимых функций; g_i(t) – обобщенные координаты, являющиеся функциями времени;  – индексы форм упругих колебаний крыла. При этом приращение относительной деформации при изгибе крыла запишется в форме

где ; Δl(Z_k) – удлинение крыла в сечении Z_k; m = m(Z) – погонная масса крыла; EJ = EJ(Z) – жесткость крыла на изгиб; h(Z) – высота профиля по размаху крыла. Для деформации кручения получим выражение

где B = (B₁, B₂, …, B_N);  – жесткость крыла на кручение;  – погонный массовый момент инерции. Теперь расчет допустимого времени θ₀ действия силы Δn_y сводится к отысканию таких удлинений ΔE₁ и ΔE₂ [74], при которых выполняются равенства

где E₁, E₂ – относительные удлинения от деформаций изгиба и кручения крыла соответственно. Наименьшее значение θ¹₀ или θ²₀, определенное из этих равенств, является искомой величиной θ₀.

При заданных значениях Γ_i, B_i определим обобщенные координаты g_i и их вторые производные . Рассмотрим продольное движение ЛА. Предполагая малость упругих деформаций, используя гипотезу стационарности, а также учитывая, что углы атаки докритические, получим:

где

θ,  – угол тангажа и скорость его изменения соответственно. В последних соотношениях величины угла кручения φ(Z,t) и изгиба Y(Z,t) крыла определяются из системы уравнений [47]

где OZ – ось жесткости крыла; m(z) – погонная масса крыла; n(Z) – положение центра тяжести каждого сечения крыла, отсчитываемого от оси жесткости; , M(Z) – погонные внешние силы и моменты;

– сила и момент конструктивного демпфирования соответственно; , F_BH = b(Z)Δn_y(Z)G – m(Z)g – момент и сила внешних воздействий; Δn_y(Z) – приращение перегрузки в сечении Z.

Функции Y(Z,t), φ(Z,t) должны удовлетворять граничным условиям вида

где l – полуразмах крыла. Кроме того, должны удовлетворяться условия сопряжения в местах скачкообразного изменения жесткостей, в местах сочленения крыла и фюзеляжа, условиям скачков сил и моментов в точках крепления сосредоточенного груза на крыле.

Для ЛА, обладающих большой массой и большими моментами инерции, ограничимся только упругими колебаниями, предполагая, что влияние угловых и линейных перемещений на общий характер движения пренебрежимо мало. В этом случае решение задачи сводится к исследованию системы (5.31) с соответствующими начальными и краевыми условиями.

Как и выше, упругие перемещения представим в виде следующих разложений:

где функции f_i(Z), φ_i(Z) образуют некоторую полную систему координатных функций, которая удовлетворяет всем граничным условиям, условиям сопряжения и скачков рассматриваемой задачи. В результате применяется метод Бубнова-Галеркина по пространственным координатам, т. е. решение представляется в виде линейной комбинации базисных функций с неопределенными коэффициентами q_i(t), зависящими от времени t. Эти неопределенные коэффициенты выбираются так, чтобы соответствующая невязка была ортогональна всем базисным функциям. Следуя такому требованию, приходим к системе обыкновенных дифференциальных уравнений относительно неопределенных коэффициентов следующего вида:

где q(t) = (q_–1, q₀, q₁, …, q_N); A – матрица коэффициентов инерции; C – матрица коэффициентов жесткости конструкции; G – матрица коэффициентов, зависящих от реактивных и инерционных сил, обусловленных продольным ускорением; B – матрица коэффициентов аэродинамической жесткости; D – матрица коэффициентов аэродинамического демпфирования; D^* — матрица коэффициентов конструктивного демпфирования; E – матрица коэффициентов демпфирования, обусловленного силами Кориолиса, F – вектор внешних сил и моментов согласно (5.31). Эту систему запишем в операторном виде (с помощью оператора Лапласа S):

[AS² + (D + D^* + E)S + (B + C + G)]g(S) = F(S),

и сведем к одному дифференциальному уравнению n-го порядка. Предполагая, что внешняя сила постоянной амплитуды воздействует на ограниченном интервале времени θ₀, получим значения обобщенных координат g_i для единичной возмущающей силы в виде

где

Так как исходный режим является установившимся, то в начальный момент времени t = 0 все производные для g_i(t) равны нулю. Это позволяет определить все постоянные A, B, φ_k1, φ_k2, φ_k в выражении (5.32). При этом начало действия дополнительной силы совпадает в момент времени t = 0.

Допустимые величины (E₁)_доп и (E₂)_доп относительного удлинения крыла при статическом нагружении, входящие в уравнение (5.30), предполагаются известными из расчетов или натурных испытаний. Тогда уравнения (5.30) запишутся в виде

где

индекс i означает, что коэффициент определен для i-й компоненты вектора g = (g₁, g₂, …, g_i, …, g_n).

Уравнения (5.33) являются трансцендентными относительно θ₀, точное решение их в аналитическом виде установить затруднительно. При приближенном решении этих уравнений будем отбрасывать высокие тона упругих колебаний. В частном случае рассмотрим первый тон изгибных колебаний крыла Y = f(z)g(t). Силы и моменты от внешних возмущений, входящие в уравнения (5.31), запишутся так:

В результате представим исходное уравнение следующим образом:

где

Вычислим f(z) приближенно из уравнений колебаний крыла постоянного сечения в пустоте. Опуская промежуточные выкладки, получим

Подставив выражения для коэффициентов из (5.34), найдем

При этом решение уравнения (5.34) запишется в виде

где ; . Для расчета θ₀ определим время t^*, при котором g(t) достигает максимума. При этом t^* величина ΔE достигает максимального значения. Приравняв к нулю производную по времени от g(t) и проделав необходимые преобразования, получим

Осталось определить

При постоянном сечении крыла, когда m = const и EJ = const, имеем .

Определив вторую производную и подставив ее в (5.33), придем к уравнению для определения θ₀

Решение поставленной задачи упрощается, если θ₀ мало, а t – θ₀  0. Предполагая, что

получим

где t^* определяется из соотношения (5.35), которое с учетом принятых допущений принимает вид .

В результате численных расчетов допустимой продолжительности θ₀ выброса перегрузки получены графики зависимостей θ₀ = f₁(F_b,EJ), θ₀ = f(F_b,V), приведенные на рис. 5.4 и 5.5. Изменения скорости полета и жесткости крыла различным образом влияют на изменение θ₀. При приближенных расчетах изменением скорости полета на θ₀ можно пренебречь.

Таким образом, установлена связь между воздействием быстроисчезающих сил и допустимой длительностью выброса перегрузки для крыла ЛА с учетом процессов динамики полета, аэродинамики обтекания и упругих деформаций изгиба и кручения.

Рис 5.4

Рис 5.5