Numerical method for structural and parametric identification of a mathematical model of incomplete inverse deformation of creep strain

Full Text

Введение

Во многих работах, например [1–4], отмечалось существенное влияние случайных микронеоднородных возмущений механических характеристик материала на поля напряжений и деформаций и, как следствие этого, — построение соответствующих стохастических реологических моделей для расчетов на прочность. Особую актуальность эта проблема приобретает при построении моделей ползучести, где разброс реологической деформации составляет до 50–70 %, и такие результаты приходится рассматривать как приемлемые [1, 4–6]. Как правило, построение стохастических уравнений ползучести базируется на обобщении соответствующих детерминированных соотношений, в которых часть параметров и функций полагаются случайными, а остальные — детерминированными [5–9]. С использованием стохастических уравнений установившейся ползучести решен ряд краевых задач для элементов конструкций из микронеоднородных материалов [10–15]. Анализ экспериментальных данных по ползучести материалов, представленных в работах [1, 16–22], показывает, что поведение реономных материалов отличается большим разнообразием. Описание всех особенностей их деформирования не укладывается в рамки наиболее употребительных определяющих уравнений на основе теорий упрочнения, течения, старения, кинетических уравнений Ю. Н. Работнова и других [23]. Основной недостаток подавляющего большинства теорий ползучести заключается в невозможности описать частичную обратимость деформации ползучести после полной разгрузки. Разновидности теории наследственности не приспособлены для описания отмеченного эффекта, поскольку теоретическая величина возврата оказывается слишком большой по сравнению с экспериментально наблюдаемой [19, 20]. К тому же область применимости большинства теорий ограничена случаем подобия стационарных кривых ползучести в целом. Вариант теории ползучести в пределах первой и второй стадий, учитывающей частичную обратимость реологической деформации, предложен Ю. П. Самариным [22] и обобщен в более поздней работе [20] для учета всех трех стадий ползучести.

Однако все представленные в работах [1, 16–23] уравнения состояния являются детерминированными, процедуры идентификации параметров моделей, как правило, являются эвристическими и многоступенчатыми с использованием некоторого набора характерных точек на кривых стационарной ползучести. Отсутствие алгоритмов статистической обработки экспериментальных данных и статистического анализа результатов вычислений не позволяют строго оценить адекватность построенных моделей ползучести и достоверность оценок их параметров. В связи с этим целью данной работы является разработка метода идентификации параметров реологической модели неполной обратимости деформации ползучести [22], который принципиально устраняет вышеизложенные недостатки известных методов.

1. Численный метод параметрической идентификации аппроксимации кривых ползучести при постоянных напряжениях

В соответствии с теорией неполной обратимости деформации ползучести при переменных напряжениях $\sigma=\sigma(t) \geqslant 0$ при наличии первой и второй стадий имеем следующую математическую модель [22]:
\[ \begin{equation}
\begin{gathered}
p(t)=u(t)+v(t)+w(t), \\
\begin{cases}
\displaystyle
u(t)= \sum\limits_{i=1}^s u_i(t),\\
\dot u_i(t)=\alpha_i \Bigl[ \beta b_i \varphi_i \Bigl(\dfrac{\sigma(t)}{\sigma_*}\Bigr)-u_i(t)\Bigr],\\
\displaystyle
v(t)=\sum\limits_{i=1}^s v_i(t),\\
\dot v_i(t)=
\begin{cases}
\alpha_i \Bigl[ (1-\beta) b_i \varphi_i \Bigl(\dfrac{\sigma(t)}{\sigma_*}\Bigr)-v_i\Bigl(\dfrac{\sigma(t)}{\sigma_*}\Bigr)\Bigr], \\
(1-\beta)b_i \varphi_i \Bigl(\dfrac{\sigma(t)}{\sigma_*}\Bigr) > v_i(t), \\
0, (1-\beta)b_i \varphi_i \Bigl(\dfrac{\sigma(t)}{\sigma_*}\Bigr) \leqslant v_i(t); \\
\end{cases}
\end{cases} \\
\dot w(t)=f \Bigl(\frac{\sigma(t)}{\sigma_*}\Bigr),
\end{gathered}
\end{equation} \tag{1} \]
где $p(t)$ — деформации ползучести; $u(t)$, $v(t)$, $w(t)$ — вязкоупругая, вязкопластическая и вязкая компонента $p(t)$ соответственно; $u_i(t)$, $v_i(t)$ — компоненты вязкоупругой и вязкопластической компонент соответственно; $\sigma_*$, $\alpha_i$, $b_i$, $\beta \in [0{,}1]$ — параметры модели; $\varphi_i ( {\sigma}/{\sigma_*} )$ и $f ( {\sigma}/{\sigma_*} )$ — функции от безразмерных напряжений.

Из (1) следует, что $u(t)+v(t)$ описывает деформацию, накопленную на первой стадии ползучести при постоянном напряжении, при этом при $\beta=1$ эта деформация полностью обратима после разгрузки, а при $\beta=0$ — полностью необратима.

Исходной информацией для определения параметров и функций соотношений (1) являются экспериментально полученные кривые стационарной ползучести при действии постоянного напряжения с последующей разгрузкой:
\[ \begin{equation*}
\sigma (t) =
\begin{cases}
\sigma^j, & 0 < t \leqslant T, \\
0, & t > T, \; j=\overline{1,m}.
\end{cases}
\end{equation*} \]
$m$ — количество кривых ползучести, $\sigma^j=\mathrm{const}$. Величина $\sigma_*$ в (1) полагается равной одному из напряжений $\sigma^j$, т.е. $\sigma_*=\sigma^j$, $1 \leqslant j \leqslant m$.

Модель (1) применима при следующих ограничениях: материал является нестареющим; напряжение $\sigma(t)$ не изменяет знака и не достигает предела текучести, третья стадия отсутствует, температура постоянна.

При постоянном напряжении $\sigma=\sigma^j= \mathrm{const}$ $(0 \leqslant t \leqslant T )$ интегрирование (1) дает аналитическое выражение для кривой ползучести вида
\[ \begin{equation}
p(t)= \sum\limits_{i=1}^s a_i \bigl[1-\exp(-\alpha_it)\bigr]+ct,
\end{equation} \tag{2} \]
где $a_i=b_i \varphi_i ( {\sigma^j}/{\sigma_*})$, $c=f({\sigma^j}/{\sigma_*})$.

Ставится задача достоверной оценки параметров аппроксимации (2) по результатам эксперимента $p(t_k)$, $k=0, 1, 2, \dots, N-1$, где $N$ — объем выборки результатов наблюдений, полученных в ходе промышленного или научно-технического эксперимента. Аппроксимация (2) содержит $2s+1$ параметров (как правило, количество экспоненциальных составляющих $s \leqslant 3$), достоверная оценка которых по результатам наблюдений является основной проблемой параметрической идентификации. При этом одновременно решается задача структурной идентификации, под которой понимается оценка числа экспоненциальных составляющих $s$ в модели (2).

Одним из известных методов решения этой задачи является метод, в основе которого лежит последовательное выделение экспоненциальных составляющих при аппроксимации кривых ползучести [24]. Однако этот метод имеет ряд существенных недостатков: требование выпуклости и монотонности функции, описывающей экспериментальную кривую ползучести, что не всегда выполняется в практике эксперимента и предварительно необходимо выполнять методы непараметрического выравнивания опытных данных; применение интерполяции к предварительно сглаженным экспериментальным данным, что существенно искажает оценки параметров экспоненциальных составляющих при наличии случайной помехи в результатах наблюдений. Однако основным недостатком этого метода является детерминированный подход к решению задачи идентификации.

Принципиально устранить указанный недостаток можно на основе численных методов нелинейного регрессионного анализа [26, 27]. При таком подходе задача параметрической идентификации решается на основе нелинейной регрессионной зависимости, построенной с учетом функциональной зависимости (2):
\[ \begin{equation}
y_k= \sum\limits_{i=1}^s a_i \left(1-\exp(-\alpha_i t_k)\right)+ct_k+\varepsilon_k, \quad
k=0, 1, 2, \dots , N-1,
\end{equation} \tag{3} \]
где $y_k=p(t_k)$ — результаты эксперимента объемом $N$, $\varepsilon_k$ — случайное возмущение в результатах наблюдений, $0 \leqslant k \leqslant N-1$. При этом среднеквадратичные оценки параметров модели (3) находятся из условия минимизации величины отклонения построенной модели $\hat{y}_k=\displaystyle \sum\limits_{i=1}^s \hat{a}_i \bigl(1-\exp(-\hat{\alpha}_i t_k)\bigr)+\hat{c}t_k$ от результатов наблюдения $y_k$, $0 \leqslant k \leqslant N-1$:
\[ \begin{equation}
\|y-\hat{y}\|^2=\|e\|^2 \to \min.
\end{equation} \tag{4} \]

Нелинейный характер регрессионной модели (3) обуславливает одну из основных проблем при оценивании параметров модели — проблему выбора начального приближения вектора оценок, обеспечивающего сходимость итерационной процедуры уточнения среднеквадратичных оценок к параметрам, удовлетворяющим критерию (4). Кроме того, большое число идентифицируемых переменных, а также неравномерность вклада различных экспоненциальных составляющих в наблюдаемый выходной сигнал может привести к проблеме вычислительной устойчивости оценок, без решения которой невозможно не только обеспечить достоверность результатов, но даже реализовать необходимые алгоритмы вычислений.

В работах [28, 29] задача параметрической идентификации аппроксимации кривой ползучести с экспоненциальным ядром решается на основе разностных уравнений. При таком подходе строится модель временного ряда, описывающая последовательность результатов наблюдений деформации ползучести. Известные соотношения между коэффициентами разностного уравнения и параметрами нелинейной регрессии позволяют свести задачу к оценке коэффициентов линейно-параметрической дискретной модели [30]. Это дает возможность не только решить проблему выбора начального приближения вектора оценок, но и непосредственно находить оценки параметров модели кривой ползучести, а также проводить статистический анализ результатов вычислений.

По сравнению с известным методом [22] построенная на основе разностных уравнений модель более адекватна результатам наблюдений, среднеквадратичные оценки параметров более достоверны и для них могут быть найдены доверительные интервалы, однако, несмотря на это, проблема устойчивости вычисления оценок остается по-прежнему актуальной. Кроме того, метод [28, 29] не позволяет решить задачу структурной идентификации — определить число экспоненциальных составляющих в модели на основе статистических методов обработки результатов наблюдений.

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

• последовательное выделение из модели (3) каждой экспоненциальной составляющей $a_i \exp(-\alpha_i t_k)$, $i=1, 2, 3$, как в известном методе [22];
• параметрическая идентификация на основе разностных уравнений каждой отдельной экспоненциальной составляющей, в том числе статистический анализ результатов расчета.

Такой подход не только позволяет существенно повысить устойчивость вычислений, а, следовательно, и эффективность оценок параметров $a_i$ и $\alpha_i$, но и дает возможность определить количество составляющих, в совокупности обеспечивающих заданную адекватность построенной модели $\hat{y}(t)$ результатам эксперимента (например, не менее $R^2=97\,\%$, где $R^2$ — коэффициент детерминации).

Следует отметить, что поскольку в алгоритме метода используются разностные уравнения, описывающие экспоненциальные составляющие, перед его применением следует провести предварительную обработку кривой ползучести, которая заключается в формировании равномерной выборки результатов эксперимента $y_k$, $0 \leqslant k \leqslant N-1$, с шагом $\tau$, равным наименьшему промежутку времени измерений деформации ползучести. Объем выборки $N$, сформированной в результате предварительной обработки, вычисляется по формуле $N= [ {t_{\max}}/{\tau}]$, где $t_{\max}$ — заданное время наблюдений.

Предлагаемый численный метод включает следующие основные этапы.

Этап 1. На первом этапе строится аппроксимация второй стадии ползучести:
\[ \begin{equation}
\hat{y}^{}_{1k}=A+c t_k.
\end{equation} \tag{5} \]
В равенстве (5) параметр $A$ равен сумме коэффициентов в экспоненциальных составляющих: $A=\sum_{i=1}^s a_i$. Параметрическая идентификация линейной модели (5) проводится на множестве точек $t_k$, $n_1\leqslant k\leqslant N-1$, в которых экспоненциальной составляющей можно пренебречь:
\[ \begin{equation}
y_k=A-\sum\limits_{i=1}^sa_i \exp (-\alpha_i t_k ) + ct_k + \varepsilon_k \approx A+ct_k+\varepsilon_k,
\quad k= n_1 , n_1+1 , \dots, N-1.
\end{equation} \tag{6} \]

Момент времени $t_k$ $(k \geqslant n_1)$, начиная с которого экспоненты считаются полностью затухающими и не влияют на вторую стадию, определяется из условия
\[ \begin{equation*}
{\Delta y'_k \over \max |y'_k|} \approx {y_k-2y_{k-1}+y_{k-2} \over y_1}< 0.001, \quad
k= 2, 3, \dots, N-1.
\end{equation*} \]

Среднеквадратичные оценки параметров линейной регрессии (6) находятся из условия
\[ \begin{equation*}
\|y-\hat{y}_1\|^2=\sum \limits_{k=n_1}^{N-1} (y_k-\hat{y}_{1k})^2 =\sum \limits_{k=n_1}^{N-1}(y_k-\hat{A}-\hat{c}t_k)^2 \to \min
\end{equation*} \]
на основе решения нормальной системы уравнений
\[ \begin{equation*}
\begin{cases}
\displaystyle
\hat{A} \cdot (N-n_1 )+\hat{c} \cdot \sum _{k=n_1}^{N-1} t_k = \sum _{k=n_1}^{N-1} y_k; \\
\displaystyle
\hat{A} \cdot \sum _{k=n_1}^{N-1} t_k +\hat{c} \cdot \sum _{k=n_1}^{N-1} t_k^2 = \sum _{k=n_1}^{N-1}t_k y_k
\end{cases}
\end{equation*} \]
по формулам
\[ \begin{equation}
\begin{array}{l}
\displaystyle
\hat{A} = {\sum \limits_{k=n_1}^{N-1}y_k \cdot \sum \limits_{k=n_1}^{N-1}t^2_k - \sum \limits_{k=n_1}^{N-1}t_ky_k \cdot \sum \limits_{k=n_1}^{N-1}t_k\over (N-n_1 )\sum \limits_{k=n_1}^{N-1}t^2_k -
\Bigl(\sum \limits_{k=n_1}^{N-1}t_k\Bigr)^2}, \\
\displaystyle
\hat{c}={ (N-n_1 )\sum \limits_{k=n_1}^{N-1}t_ky_k - \sum \limits_{k=n_1}^{N-1}y_k \cdot \sum \limits_{k=n_1}^{N-1}t_k \over \sum \limits_{k=n_1}^{N-1}t^2_k - \Bigl(\sum \limits_{k=n_1}^{N-1}t_k\Bigr)^2}.
\end{array}
\end{equation} \tag{7} \]

Результатом расчетов, выполненных на первом этапе алгоритма рассматриваемого численного метода, является построенная аппроксимация для второй стадии ползучести:
\[ \begin{equation}
\hat{y}_{1k}=\hat{A}+\hat{c}t_k, \quad k=0, 1, 2, \dots, N-1.
\end{equation} \tag{8} \]

Этап 2. На втором этапе находятся среднеквадратичные оценки параметров первой экспоненциальной составляющей:
\[ \begin{equation}
f_k=-a_1 \exp [-\alpha_1 \tau (k-1)], \quad k=1, 2, \dots, N.
\end{equation} \tag{9} \]

Вначале формируется выборка результатов расчета $y^{(1)}_k=y_k-\hat{y}_{1k}$, $k=0, 1, 2, \dots, N-1$. Среднеквадратичные оценки параметров $a_1$ и $\alpha_1$ вычисляются на основе нелинейной регрессии
\[ \begin{equation*}
y^{(1)}_k=- a_1 \exp (- \alpha_1 \tau k) + \varepsilon_k, \quad k=0, 1, 2, \dots, N-1,
\end{equation*} \]
из условия $\|\hat{\varepsilon}\|^2=\displaystyle \sum \limits_{k=0}^{N-1}\bigl[y^{(1)}_k+\hat{a}_1 \exp (-\hat{\alpha}_1 \tau k)\bigr]^2 \to \min$.

В формате данного метода эта задача решается с использованием разностного уравнения, построенного на основе дискретной функции (9):
\[ \begin{equation}
\begin{cases}
f_1=\lambda_2; &\\
f_k=\lambda_1f_{k-1}, & k= 1, 2, \dots, N,
\end{cases}
\end{equation} \tag{10} \]
где $\lambda_1= \exp(-\alpha_1 \tau)$, $\lambda_2=-a_1$.

Введем обозначения: $z_k=y^{(1)}_k$, $0 \leqslant k\leqslant N-1$. Тогда с учетом равенства $f_{k+1}=z_k-\varepsilon_k$, $0 \leqslant k\leqslant N-1$, из системы (10) получаем
\[ \begin{equation}
\begin{cases}
z_0=\lambda_2+\eta_0; &\\
z_k=\lambda_1 z_{k-1}+\eta_k, & k= 1, 2, \dots, N-1;\\
\eta_0=\varepsilon_0; &\\
\eta_k=\varepsilon_k-\lambda_1\varepsilon_{k-1}, & k= 1, 2, \dots, N-1
\end{cases}
\end{equation} \tag{11} \]
или в матричной форме:
\[ \begin{equation*}
z=F \lambda+\eta,
\end{equation*} \]
где $z=\bigl(z_0,z_1, \dots, z_{N-1}\bigr)^\top$, $\eta=\bigl(\varepsilon_0, \varepsilon_1-\lambda_1 \varepsilon_0, \varepsilon_2-\lambda_1\varepsilon_1, \dots, \varepsilon_{N-1}-\lambda_1\varepsilon_{N-2}\bigr)^\top$, $\lambda=\left(\lambda_1, \lambda_2\right)^\top$, $F=
\begin{bmatrix}
1& 0 & 0 & \dots & 0\\
0& z_0 & z_1 & \dots & z_{N-2}
\end{bmatrix}^\top$, ${}^\top$ — символ транспонирования.

Начальные оценки параметров нелинейной регрессии (11) находятся из условия $\|\hat{\eta}\|^2=\|z-F\lambda\|^2 \to \min$ по формуле $\hat{\lambda}^{(0)}= (F^\top F )^{-1}F^\top z$, откуда получаем оценки
\[ \begin{equation}
\hat{\lambda}_1^{(0)}= {\sum\limits_{k=1}^{N-1}z_kz_{k-1}\over \sum\limits_{k=1}^{N-1}z^2_{k-1}}, \qquad \hat{\lambda}_2^{(0)}=z_0.
\end{equation} \tag{12} \]

Для уточнения среднеквадратичных оценок параметров математической модели (11) с учетом выполнения условия $\|\hat{\varepsilon}\|^2 \to \min$ рассмотрим нелинейную регрессию относительно коэффициентов разностного уравнения
\[ \begin{equation}
z_{k-1}=f_k(\lambda_1, \lambda_2)+ \varepsilon_k, \quad k= 1, 2, \dots, N,
\end{equation} \tag{13} \]
где нелинейная дискретная функция $f_k(\lambda_1, \lambda_2)$ описывается рекуррентной формулой (10).

Итерационная процедура уточнения среднеквадратичных оценок нелинейной регрессии (13) описывается формулой [26, 27, 30]
\[ \begin{multline}
\hat{\lambda}^{(i+1)}=\hat{\lambda}^{(i)}+
\bigl[ \bigl(W (\hat{\lambda}^{(i)} )\bigr)^\top W (\hat{\lambda}^{(i)} )\bigr]^{-1}
\bigl(W (\hat{\lambda}^{(i)} )\bigr)^\top
\bigl[y^{(1)}-f (\hat{\lambda}^{(i)} )\bigr],
\\
i=0, 1, 2, \dots,
\end{multline} \tag{14} \]
где $W (\lambda_1, \lambda_2 )$ — матрица Якоби размера $ [N {\times} 2 ]$, элементы которой имеют вид
\[ \begin{equation*}
\begin{array}{l}
\displaystyle
w_{1k} (\lambda_1, \lambda_2 )={\partial f_k (\lambda_1, \lambda_2 )\over \partial \lambda_1}=
\begin{cases}
\qquad \qquad \qquad 0, & k=1;\\
f_{k-1} (\lambda_1, \lambda_2 )+\lambda_1w_{1,k-1} (\lambda_1, \lambda_2 ), &
2\leqslant k\leqslant N;
\end{cases}
\\
\displaystyle
w_{2k} (\lambda_1, \lambda_2 )={\partial f_k (\lambda_1, \lambda_2 )\over \partial \lambda_2}=
\begin{cases}
\qquad \quad 0, & k=1;\\
\lambda_1w_{2,k-1} (\lambda_1, \lambda_2 ),& 2\leqslant k\leqslant N.
\end{cases}
\end{array}
\end{equation*} \]

Начальное приближение в рекуррентной формуле (14) находится по формулам (12). Процесс уточнения среднеквадратичных оценок $\hat{\lambda}^{(i)}$ заканчивается при выполнении условия
\[ \begin{equation}
\|\hat{\lambda}^{(i+1)}-\hat{\lambda}^{(i)}\|<0.001 \|\hat{\lambda}^{(i)}\|.
\end{equation} \tag{15} \]

С учетом полученных среднеквадратичных оценок коэффициентов разностного уравнения можно найти оценки параметров первой экспоненциальной составляющей:
\[ \begin{equation}
\hat{a}_1=-\hat{\lambda}^{(i)}_2, \quad \hat{\alpha}_1=-\dfrac{1}{\tau} \ln \hat{\lambda}^{(i)}_1.
\end{equation} \tag{16} \]
и записать ее зависимость в виде
\[ \begin{equation}
\hat{y}^{(1)}_k=-\hat{a}_1 \exp (-\hat{\alpha}_1t_k), \quad k= 0, 1, 2, \dots, N-1.
\end{equation} \tag{17} \]

На основе полученных выше результатов можно перейти к решению задачи уточнения среднеквадратичных оценок параметров нелинейной регрессии (3) при использовании одной экспоненциальной составляющей.

На этом шаге с учетом найденных ранее оценок коэффициентов разностного уравнения $\hat{\lambda}^{(i)}_1$ и $\hat{\lambda}^{(i)}_2$ (формула (14)), а также оценки $\hat{c}$ (формула (7)) вычисляется вектор $\hat{\mu}^{(0)}= \bigl(\hat{a}^{(0)}_1, \hat{\alpha}^{(0)}_1, \hat{c}^{(0)}\bigr)^\top$ начальных оценок параметров функциональной зависимости
\[ \begin{equation*}
f_k(\mu)=a_1\bigl[1-\exp (-\alpha_1t_k)\bigr]+ct_k, \quad k=0, 1, 2, \dots, N-1,
\end{equation*} \]
по формулам
\[ \begin{equation*}
\hat{a}^{(0)}_1=-\hat{\lambda}^{(i)}_2, \quad
\hat{\alpha}^{(0)}_1=-\dfrac{1}{\tau} \ln \hat{\lambda}^{(i)}_1, \quad
\hat{c}^{(0)}=\hat{c}.
\end{equation*} \]

Далее находятся элементы матрицы Якоби $W(\mu)$ размера $ [N {\times} 3]$:
\[ \begin{equation*}
\begin{array}{l}
w_{k1}(\mu)=\dfrac{\partial f_k (\mu)}{\partial a_1}=1-\exp(-\alpha_1t_k), \quad k=0, 1, \dots, N-1; \\
w_{k2}(\mu)=\dfrac{\partial f_k (\mu)}{\partial \alpha_1}=a_1t_k\exp(-\alpha_1t_k), \quad k=0, 1, \dots, N-1; \\
w_{k3}(\mu)=\dfrac{\partial f_k (\mu)}{\partial c}=t_k, \quad k=0, 1, \dots, N-1.
\end{array}
\end{equation*} \]

Итерационная процедура уточнения среднеквадратичных оценок нелинейной регрессии (3) при использовании одной экспоненциальной составляющей описывается формулой
\[ \begin{multline}
\hat{\mu}^{(i+1)}=\hat{\mu}^{(i)}+ \bigl[\bigl(W (\hat{\mu}^{(i)} )\bigr)^\top W
(\hat{\mu}^{(i)} )\bigr]^{-1}
\bigl(W (\hat{\mu}^{(i)} ) \bigr)^\top
\bigl[y-f (\hat{\mu}^{(i)} )\bigr],
\\
i=0, 1, 2, \dots.
\end{multline} \tag{18} \]

Процесс уточнения среднеквадратичных оценок $\hat{\mu}^{(i)}$ заканчивается при выполнении условия
\[ \begin{equation}
\|\hat{\mu}^{(i+1)}-\hat{\mu}^{(i)}\|< 0.001\|\hat{\mu}^{(i)}\|.
\end{equation} \tag{19} \]

Результатом расчетов, выполненных на этом шаге алгоритма рассматриваемого численного метода, является построенная аппроксимация кривой ползучести с одной экспоненциальной составляющей:
\[ \begin{equation}
\hat{y}_k=\hat{a}_1\bigl[1-\exp(-\hat{\alpha}_1t_k)\bigr]+\hat{c}t_k, \quad
k=0, 1, 2, \dots, N-1,
\end{equation} \tag{20} \]
в которой оценки параметров соответствуют элементам вектора $\hat{\mu}^{(i+1)}$:
\[ \begin{equation*}
\hat{a}_1=\hat{\mu}_1^{(i+1)},\quad
\hat{\alpha}_1=\hat{\mu}_2^{(i+1)},\quad
\hat{c}=\hat{\mu}^{(i+1)}_3.
\end{equation*} \]

Для вычисления дисперсий величин $a_1$, $\alpha_1$, $c$ находится матрица
\[ \begin{equation*}
V[\mu]= (W^\top W )^{-1}s^2_{\text{ост}},
\end{equation*} \]
где $s^2_{\text{ост}}=\displaystyle \frac{1}{N-3} \sum\limits_{i=0}^{N-1} (y_k-\hat{y}_k )^2$.

Диагональные элементы $V_{11}(\mu)=D(a_1)$, $V_{22}(\mu)=D(\alpha_1)$ и $V_{33}(\mu)=D(c)$ задают дисперсию параметров $a_1$, $\alpha_1$ и $c$.

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

Например, можно использовать относительную величину среднеквадратической нормы отклонения аппроксимации (20) от данных эксперимента:
\[ \begin{equation}
h, \,\% = \Biggl( {\sum \limits_{k=0}^{N-1}\left(y_k-\hat{y}_k\right)^2\over \sum \limits_{k=0}^{N-1}y^2_k } \Biggr)^{1/2} \cdot 100 \,\%.
\end{equation} \tag{21} \]

В качестве другого критерия адекватности модели можно использовать «коэффициент детерминации» $R^2$ %, который определяется через выборочный коэффициент множественной корреляции (коэффициент корреляции между случайными векторами $y$ и $\hat{y}$:
\[ \begin{equation}
R={\sum \limits_{k=0}^{N-1} (y_k-\overline{y} ) (\hat{y}_k-\overline{\hat{y}} )
\over
\biggl({\sum \limits_{k=0}^{N-1} (y_k-\overline{y})^2}\biggr)^{1/2}
\biggl({\sum \limits_{k=0}^{N-1} (\hat{y}_k-\overline{\hat{y}} )^2\biggr)^{1/2}}},
\quad
R^2, \,\% = R^2 \cdot 100 \,\%.
\end{equation} \tag{22} \]

Условием адекватности построенной модели результатам эксперимента можно считать выполнение одного из неравенств
\[ \begin{equation}
h, \% \leqslant 3\,\%, \qquad R^2, \% \geqslant 97\,\%.
\end{equation} \tag{23} \]

Этап 3. Если условия (23) не выполняются, то следует перейти к третьему этапу численного метода — выделению в соотношении для кривой ползучести (3) второй экспоненциальной составляющей: $a_2\bigl[1-\exp(-\alpha_2 t_k)\bigr]$.

При решении этой задачи используются результаты, полученные на первых двух этапах: построенные аппроксимации второй стадии ползучести (8) и первой экспоненциальной составляющей (17). На первом шаге этого этапа формируется выборка результатов расчета $y^{(2)}_k=y_k-\hat{y}_{1k}-\hat{y}^{(1)}_k$, $k =0, 1, 2, \dots, N-1$. Для этой выборки рассматривается аппроксимация второй экспоненциальной составляющей вида
\[ \begin{equation*}
f_k=-a_2 \exp \bigl[-\alpha_2 \tau (k-1)\bigr], \quad k=1, 2, \dots, N.
\end{equation*} \]

Среднеквадратичные оценки параметров $a_2$ и $\alpha_2$ вычисляются с применением нелинейной регрессии
\[ \begin{equation*}
y^{(2)}_k=- a_2 \exp (- \alpha_2 \tau k) + \varepsilon_k, \quad k=0, 1, 2, \dots, N-1,
\end{equation*} \]
из условия $\|\hat{\varepsilon}\|^2= \displaystyle \sum \limits_{k=0}^{N-1}\bigl[y^{(2)}_k+\hat{a}_2 \exp (-\hat{\alpha}_2 \tau k)\bigr]^2 \to \min$.

Эта задача решается на основе алгоритма, описанного выше на втором этапе численного метода, с использованием формул (10)–(15) с учетом обозначений $z_k=y^{(2)}_k$, $0 \leqslant k \leqslant N-1$, $\lambda_1=\exp(-\alpha_2 \tau)$ и $\lambda_2=-a_2$.

Благодаря полученным среднеквадратичным оценкам коэффициентов разностного уравнения можно найти оценки параметров второй экспоненциальной составляющей:
\[ \begin{equation}
\hat{a}_2=-\hat{\lambda}^{(i)}_2, \quad \hat{\alpha}_2=-\dfrac{1}{\tau} \ln \hat{\lambda}^{(i)}_1
\end{equation} \tag{24} \]
и записать ее выражение в виде
\[ \begin{equation*}
\hat{y}^{(2)}_k=- \hat{a}_2 \exp (- \hat{\alpha}_2 \tau k) + \varepsilon_k, \quad k=0, 1, 2, \dots, N-1.
\end{equation*} \]

На основе полученных выше результатов можно перейти к решению задачи уточнения среднеквадратичных оценок параметров нелинейной регрессии (3) при использовании двух экспоненциальных составляющих.

На этом шаге с учетом найденных на предыдущих шагах оценок параметров экспоненциальных составляющих (формулы (16) и (24)), а также оценки $\hat{c}$ (формула (7)) вычисляется вектор $\hat{\mu}^{(0)}=\bigl(\hat{a}^{(0)}_1, \hat{\alpha}^{(0)}_1, \hat{a}^{(0)}_2, \hat{\alpha}^{(0)}_2, \hat{c}^{(0)}\bigr)^\top$ начальных оценок параметров функциональной зависимости
\[ \begin{equation*}
f_k(\mu)=a_1\bigl[1-\exp (-\alpha_1t_k)\bigr]+a_2\bigl[1-\exp (-\alpha_2t_k)\bigr]+ct_k, \quad k=0, 1, 2, \dots, N-1,
\end{equation*} \]
по формулам
\[ \begin{equation*}
\hat{a}^{(0)}_1=\hat{a}_1,\quad \hat{\alpha}^{(0)}_1=\hat{\alpha}_1,\quad
\hat{a}^{(0)}_2=\hat{a}_2,\quad \hat{\alpha}^{(0)}_2=\hat{\alpha}_2,\quad \hat{c}^{(0)}=\hat{c}.
\end{equation*} \]

Далее находятся элементы матрицы Якобы $W(\mu)$ размера $ [N {\times} 5]$:
\[ \begin{equation}
w_{k1}(\mu)={\partial f_k (\mu)\over \partial a_1}=1-\exp(-\alpha_1t_k), \quad k=0, 1, \dots, N-1;
\end{equation} \tag{25} \]
\[ \begin{equation}
w_{k2}(\mu)={\partial f_k (\mu)\over \partial \alpha_1}=a_1t_k\exp(-\alpha_1t_k), \quad k=0, 1, \dots, N-1;
\end{equation} \tag{26} \]
\[ \begin{equation}
w_{k3}(\mu)={\partial f_k (\mu)\over \partial a_2}=1-\exp(-\alpha_2t_k), \quad k=0, 1, \dots, N-1;
\end{equation} \tag{27} \]
\[ \begin{equation}
w_{k4}(\mu)={\partial f_k (\mu)\over \partial \alpha_2}=a_2t_k\exp(-\alpha_2t_k), \quad k=0, 1, \dots, N-1;
\end{equation} \tag{28} \]
\[ \begin{equation}
w_{k5}(\mu)={\partial f_k (\mu)\over \partial c}=t_k, \quad k=0, 1, \dots, N-1.
\end{equation} \tag{29} \]

Итерационная процедура уточнения среднеквадратичных оценок нелинейной регрессии (3) при использовании двух экспоненциальных составляющих описывается формулой (18). Процесс уточнения среднеквадратичных оценок $\hat{\mu}^{(i)}$ заканчивается при выполнении условия (19). Результатом расчетов, выполненных на этом шаге алгоритма рассматриваемого численного метода, является построенная аппроксимация кривой ползучести с двумя экспоненциальными составляющими:
\[ \begin{equation*}
\hat{y}_k=\hat{a}_1\bigl[1-\exp(-\hat{\alpha}_1t_k)\bigr]+
\hat{a}_2\bigl[1-\exp(-\hat{\alpha}_2t_k)\bigr]+\hat{c}t_k, \quad k=0, 1, 2,\dots, N-1,
\end{equation*} \]
в которой оценки параметров соответствуют элементам вектора $\hat{\mu}^{(i+1)}$:
\[ \begin{equation*}
\hat{a}_1=\hat{\mu}_1^{(i+1)},\quad \hat{\alpha}_1=\hat{\mu}_2^{(i+1)},\quad
\hat{a}_2=\hat{\mu}_3^{(i+1)},\quad \hat{\alpha}_2=\hat{\mu}_4^{(i+1)}, \quad
\hat{c}=\hat{\mu}^{(i+1)}_5.
\end{equation*} \]

Для вычисления дисперсий величин $a_1$, $\alpha_1$, $a_2$, $\alpha_2$ и $c$ находится матрица $V[\mu] = (W^\top W )^{-1}s^2_{\text{ост}}$, диагональные элементы которой $V_{11}(\mu)=D(a_1)$, $V_{22}(\mu)=D(\alpha_1)$, $V_{33}(\mu)=D(a_2)$, $V_{44}(\mu)=D(\alpha_2)$, $V_{55}(\mu)=D(c)$ и задают дисперсию величин $a_1$, $\alpha_1$, $a_2$, $\alpha_2$ и $c$.

На заключительном шаге третьего этапа решается задача проверки адекватности построенной модели результатам наблюдений и, как следствие, целесообразности использования в модели третьей экспоненциальной составляющей. Для этого следует применить критерии адекватности модели, представленные соотношениями (21) или (22).

Этап 4. Если условия (23) не выполняются, то следует перейти к четвертому этапу численного метода — выделению в соотношении для кривой ползучести (3) третьей экспоненциальной составляющей: $a_3 \bigl[1-\exp (-\alpha_3 t_k)\bigr]$.

Алгоритм вычисления среднеквадратичных оценок параметров третьей экспоненциальной составляющей аналогичен процедурам, описанным выше на втором и третьем этапах численного метода. В основе этого алгоритма лежит выборка результатов расчета:
\[ \begin{equation*}
y_k^{(3)}=y_k-\hat{y}_{1k}-y_k^{(1)}-y_k^{(2)},\quad k=0, 1, 2, \dots, {N-1}.
\end{equation*} \]
При уточнении среднеквадратичных оценок параметров нелинейной регрессии (3) с тремя экспоненциальными составляющими на основе итерационной процедуры (18) матрица Якоби $W(\mu)$ имеет размер $[N {\times} 7]$ и содержит уже семь столбцов $w_{kj}$, $ 1\leqslant j \leqslant 7$, первые четыре из которых описываются формулами (25), (26), (27), (28), последний седьмой — формулой (29), а пятый и шестой — формулами
\[ \begin{equation*}
\begin{array}{c}
w_{k5}(\mu)=\dfrac{\partial f_k (\mu)}{\partial a_3}=1-\exp(-\alpha_3t_k), \quad k=0, 1, \dots, N-1; \\
w_{k6}(\mu)=\dfrac{\partial f_k (\mu)}{\partial \alpha_3}=a_3t_k\exp(-\alpha_3t_k), \quad k=0, 1, \dots, N-1.
\end{array}
\end{equation*} \]

Далее алгоритм аналогичен случаю построения первых двух экспонент.

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

2. Построение математической модели теории неполной обратимости деформации ползучести

Как отмечено выше, исходной информацией для идентификации параметров и функций в (1) является серия кривых стационарной ползучести при постоянных напряжениях $\sigma(t)=\sigma^j$, $1 \leqslant j \leqslant m$, при $t \in [{0}, T]$ с последующей разгрузкой: $\sigma(t)=0$, $t>T$.

Этап 1. На первом этапе для каждой кривой ползучести при $\sigma(t)=\sigma^j$, $t \in [0,T]$, описываемой соотношением (2), по изложенной в п. 1 методике определяются оценки величин $\hat{a}^j_i$, $\hat{\alpha}^j_i$, $\hat{c}^j_i$ и их дисперсии, $1\leqslant i \leqslant s$, $1\leqslant j\leqslant m$. Поскольку стохастически нелинейные определяющие соотношения трудно использовать при решении соответствующих краевых задач, величины $\alpha_i$ в (1) будем считать детерминированными и полагать их значения равными усредненным значениям $\alpha^j_i$ при всех уровнях напряжений:
\[ \begin{equation*}
\alpha_i=\dfrac{1}{m}\sum \limits_{j=1}^m \alpha_i^j.
\end{equation*} \]
Аналогично, по этой же причине и величину $\beta$ будем также считать детерминированной (процедура ее вычисления приводится ниже). Тогда уравнения (1) будут стохастически линейными.

Этап 2. На втором этапе находятся аппроксимационные зависимости для величин $\hat{a}_i^j$, $\hat{c}^j$, $1 \leqslant j \leqslant m$, вида
\[ \begin{equation}
\hat{a}^*_i=\hat{b}_i\Bigl(\dfrac{\sigma}{\sigma_*}\Bigr)^{n_i}, \qquad
\hat{c}^*=\hat{\gamma}\Bigl(\dfrac{\sigma}{\sigma_*}\Bigr)^{n}
\end{equation} \tag{30} \]
по методу наименьших квадратов, т.е. в (1) полагаем, что
\[ \begin{equation*}
\varphi_i\Bigl(\dfrac{\sigma}{\sigma_*}\Bigr)=\Bigl(\dfrac{\sigma}{\sigma_*}\Bigr)^{n_i}, \qquad
f\Bigl(\dfrac{\sigma}{\sigma_*}\Bigr)=\hat{\gamma} \Bigl(\dfrac{\sigma}{\sigma_*}\Bigr)^{n},
\end{equation*} \]
где $n_i$ и $n$ — детерминированные параметры.

Этап 3. На третьем этапе в качестве базовой принимается кривая ползучести при $\sigma(t)=\sigma^j=\sigma_*$, $1 \leqslant j \leqslant m$. Тогда соотношение типа (2) для этой кривой ползучести имеет вид
\[ \begin{equation}
p(t)=\sum \limits_{i=1}^s \hat{b}_i\bigl[1-\exp(-\alpha_it)\bigr]+\hat{\gamma} t,
\end{equation} \tag{31} \]
где $\hat{b}_i$ и $\hat{\gamma}$ определяются из аппроксимации (30), а величина $\hat{\alpha}_i$ находится усреднением величин $\alpha_i^j$, $1\leqslant j\leqslant m$, по всем реализациям. Поскольку величины $\hat{b}_i$ и $\hat{\gamma}$ в (31) не совпадают с аналогичными значениями $a_i$ и $c$ в (2), оценки дисперсий величин $b_i$ и $\gamma$ определим, используя дисперсии $a_i$ и $c$, полученные ранее для кривой ползучести при $\sigma=\sigma_*$ на первом этапе исследований, полагая
\[ \begin{equation}
D[b_i]=\Bigl(\dfrac{\hat{b}_i}{\hat{a}_i}\Bigr)^2 D[a_i], \quad
D[\gamma]=\Bigl(\dfrac{\hat{\gamma}}{\hat{c}}\Bigr)^2 D[c].
\end{equation} \tag{32} \]

Полученным оценкам параметров $b_i$ и $\gamma$ придается универсальный характер, т.е. они используются при любых значениях напряжений $\sigma=\sigma(t)$ в определяющих соотношениях (1).

Поскольку оценки математических ожиданий случайных величин $b_i$ и $\gamma$ и их дисперсий известны, можно построить доверительные интервалы для деформации ползучести $p(t)$ при любых законах $\sigma=\sigma(t)$ на основании (1), в частности для серии кривых ползучести при $\sigma(t)=\sigma^j$, $1\leqslant j \leqslant m$, используя стохастически линейную зависимость вида
\[ \begin{equation*}
p(t)=\sum \limits_{i=1}^s b_i \Bigl(\dfrac{\sigma}{\sigma_*}\Bigr)^{n_i}
\bigl[1-\exp(-\alpha_it)\bigr]+\gamma \Bigl(\dfrac{\sigma}{\sigma_*}\Bigr)^nt,
\end{equation*} \]
из которой нетрудно определить дисперсию величины $p(t)$ в любой момент времени.

Этап 4. На четвертом этапе определяется величина $\beta$ в соотношениях (1), которая задает величину обратимой компоненты деформации ползучести, накопленной на первой стадии ползучести после полной разгрузки образца. В дальнейшем для сохранения стохастической линейности (1) вводим гипотезу, что величина $\beta$ также является детерминированной. Если на кривых ползучести при $\sigma=\sigma^j$ наблюдается ярко выраженная вторая стадия ползучести, т.е. величина $u^j(t)+v^j(t)$ достигла асимптотического значения
\[ \begin{equation*}
\lim_{t \to \infty} [u^j(t)+v^j(t) ]=\sum \limits_{i=1}^s b_i \Bigl(\dfrac{\sigma^j}{\sigma_*}\Bigr)^{n_i},
\end{equation*} \]
то величина $\beta^j$ для этой кривой ползучести находится по формуле
\[ \begin{equation*}
\beta^j={p^j (T )-p^j (T+t^* )\over \sum \limits_{i=1}^s b_i \bigl( {\sigma^j}/{\sigma_*}\bigr)^{n_i}},
\end{equation*} \]
где $t^*$ — время окончания наблюдения за деформацией ползучести после разгрузки ($\sigma^j(t)=0$, $t \in [T, t^*]$).

В случае если вторая стадия отсутствует, т.е. в наличии только первая стадия, величина $\beta^j$ определяется в соответствии с зависимостью
\[ \begin{equation*}
\beta^j={p^j (T )-p^j (T+t_1 )\over p^j(t_1)},
\end{equation*} \]
где $t_1 \in (0, T)$.

Далее в качестве величины $\beta$ принимается осредненное значение
\[ \begin{equation*}
\beta=\dfrac{1}{m}\sum \limits_{j=1}^m \beta^j.
\end{equation*} \]

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

3. Результаты расчетов и их анализ

Применение разработанной методики построения теории неполной обратимости деформации ползучести (1) реализовано на основе большого массива экспериментальных данных из работ [19, 31]. В качестве примера приведем результаты расчетов для сплавов ВЖ98 при температуре 900°C, ЭИ437А при температуре 700°C, стали ЭИ736 при 500°C, из [31] (здесь отсутствуют данные о деформации ползучести при полной разгрузке образцов, поэтому полагалось, что $\beta=0$) и сплава ЭП693 при температуре 700°C [19], где имеются кривые обратной ползучести ($\beta \neq 0$).

Таблица 1. Расчетные оценки параметров аппроксимаций кривых ползучести для сплавов ВЖ98, ЭИ437А, ЭП693 и стали ЭИ698 при фиксированных напряжениях
[Estimated parameters for the approximation of creep curves for alloys VZh98, EI437A, EP693, and steel EI698 at fixed stresses]

$\sigma$, MPa	$\alpha_1$, $h^{-1}$	$\hat{a}_1\cdot 10^4$	$\hat{c}\cdot 10^6$	$D[a_1]\cdot 10^9$	$D[c]\cdot 10^{13}$	$\sigma_*$, MPa
EI736 steel, 500°C
170	0.2	4.67	1.9	0.324	0.74	200
200	0.5	7.74	4.36	0.72	1.9
250	0.38	9.81	6.87	1.14	2.72
EI437A alloy, 700°C
240	0.16	2.38	3.1	0.25	0.55	300
300	0.37	3.44	6.71	0.124	0.39
350	0.28	5.71	11.36	0.162	0.44
VZh98 alloy, 900°C
10	0.27	4.48	10.2	0.818	1.97	15
15	0.24	11.07	19.5	4.23	9.76
20	0.25	13.5	30	4.35	10.3
25	0.3	17.8	49.5	8.9	21.4
EP693 alloy, 700°C
245.3	0.25	2.82	8.3	0.163	0.37	294.3
294.3	0.28	3.26	16.1	1.5	3.19
343.4	0.56	4.15	38.9	1.83	12.2

Таблица 2. Расчетные оценки параметров стохастической модели (1) [Estimated parameters of the stochastic model (1)]

Material	$\hat{\alpha}_1$	$\hat{b}_1 \cdot 10^4$	$\hat{\gamma} \cdot 10^6$	$D[b_1]\cdot 10^9$	$D[\gamma]\cdot 10^{13}$	$n_1$	$n$	$\beta$
VZh98 alloy, 900°C	0.27	9.0	20.0	2.80	10.0	1.46	1.71	0
EI736 steel, 500°C	0.36	7.0	4.0	0.59	1.60	1.86	3.2	0
EI437A alloy, 700°C	0.27	3.8	7.0	0.15	0.43	2.24	3.45	0
EP693 alloy, 700°C	0.36	3.4	18	1.63	4.0	1.14	4.6	0.64

Применение разработанной методики определения оценок параметров аппроксимации (2) показало, что для всех материалов достаточно одного экспоненциального слагаемого ($s=1$). В табл. 1 приведены оценки основных параметров (2) для каждой кривой ползучести при фиксированном напряжении. Далее в соответствии с методикой строились степенные аппроксимации (30) и определялись $\hat{b}_i$, $n_1$, $\hat{\gamma}$, $n$, а для определения дисперсий $D[b_1]$ и $D[\gamma]$ пользовались формулой (32), где $\hat{a}_1$ и $\hat{c}$ соответствуют базовой кривой ползучести при $\sigma=\sigma_*$. В табл. 2 приведены окончательные оценки всех параметров определяющих уравнений состояния (1) для всех исследуемых материалов.

Анализ данных табл. 2 свидетельствует, что величины $\alpha^j_i$, $D[b^j_i]$, $D[\gamma^j]$ $(1\leqslant j \leqslant m)$ монотонно не зависят от величины напряжения, что может служить некоторым обоснованием принятой гипотезы об их фиксированной величине в соотношениях (1).

Рис. 1. Результаты первичной аппроксимации кривых ползучести ЭИ437А при температуре 700°C с данными из табл. 1 (a) и по модели (1) с данными из табл. 2 (b): 1 — $\sigma=240$ МПа, 2 — $\sigma=300$ МПа, 3 — $\sigma=350$ МПа; сплошная линия — математическое ожидание; штриховые линии — границы доверительных интервалов; маркеры (точки) — экспериментальные данные для ползучести
[Figure. 1. Results of the primary approximation of creep curves for EI437A at a temperature of 700°C with data from table 1 (a), and according to the model (1) with data from table 2 (b): 1 — $\sigma=240$ MPa, 2 — $\sigma=300$ MPa, 3 — $\sigma=350$ MPa; solid line — mathematical expectation; dashed lines — confidence interval boundaries; markers (dots) — experimental data for creep]

На рис. 1, a в качестве иллюстрации приведены результаты аппроксимации после первичной статистической обработки кривых ползучести сплава ЭИ437А при температуре 700°C для трех уровней напряжений на основании функциональной зависимости, задаваемой соотношением (2). Исходя из некоррелированности случайных величин $a_i$ и $c$, дисперсия для деформации $p(t)$ для каждой кривой в любой момент времени с данными из табл. 1 рассчитывалась по формуле
\[ \begin{equation*}
D[p(t)]=\sum \limits_{i=1}^s D[a_i]\bigl(1-\exp(-\alpha_i t)\bigr)^2+D[c]t^2,
\end{equation*} \]
а затем c использованием классической методики [27] строились 99 % доверительные интервалы для математического ожидания деформации $p(t)$. Значения математических ожиданий величин $\hat{a}_i$, $\hat{c}$ и $\hat{\alpha}_i$ и дисперсий $D[a_i]$, $D[c]$ приведены в табл. 1. На рис. 1, b и рис. 2 приведены расчетные значения для деформации ползучести по реологической модели (1) с параметрами, приведенными в табл. 2, для четырех материалов: ЭИ437А (700°C), ЭИ736 (500°C), ВЖ98 (900°C) и ЭП693 (700°C). На всех рисунках сплошная линия — математическое ожидание, штриховые линии — границы 99 % доверительного интервала, маркеры (точки) — экспериментальные значения деформации ползучести. Для построения доверительных интервалов по модели (1) расчет дисперсии для деформации ползучести осуществлялся по формуле (см. формулу (32))
\[ \begin{equation*}
D[p(t)]=\sum \limits_{i=1}^s D[b_i]
\Bigl[\Bigl(\dfrac{\sigma}{\sigma_*}\Bigr)^{n_i}
\bigl(1-\exp(-\alpha_i t)\bigr)\Bigr]^2+
D[\gamma] \Bigl[\Bigl(\dfrac{\sigma}{\sigma_*}\Bigr)^{n}t\Bigr]^2.
\end{equation*} \]

Для всех приведенных вариантов расчета экспериментальные данные не выходят за пределы соответствующих доверительных интервалов для деформации ползучести, что свидетельствует о достоверности оценок параметров модели (1) и ее адекватности.

Рис. 2. Расчетные данные деформации ползучести по модели (1) для стали ЭИ736 (500°C) (a), сплава ВЖ98 (900°C) (b), сплава ЭП693 (700°C) (c): сплошная линия — математическое ожидание; штриховые линии — границы доверительных интервалов; маркеры (точки) — экспериментальные данные; цифры — напряжения в МПа
[Figure. 2. Calculated data of creep deformation according to the model (1) for EI736 steel (500°C) (a), VZh98 alloy (900°C) (b), and EP693 alloy (700°C) (c): solid line — mathematical expectation; dashed lines — confidence interval boundaries; markers (dots) — experimental data; numbers — stress in MPa]

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

Заключение

1. Разработан новый численный метод параметрической и структурной идентификации нелинейной теории неполной обратимости деформации ползучести, позволяющий свести задачу к нелинейному регрессионному анализу определения оценок случайных параметров реологической модели на основе временных рядов последовательности результатов наблюдений деформации ползучести при различных постоянных напряжениях.
2. Показано, что при использовании даже осредненных экспериментальных кривых ползучести построенная реологическая модель на их основе является стохастической вследствие процедуры параметрической идентификации нелинейной регрессионной модели.
3. Выполнена параметрическая и структурная идентификация теории неполной обратимости деформации ползучести для ряда сталей и сплавов, получены численные значения оценок параметров для этих материалов.
4. Выполнена проверка адекватности построенной математической теории неполной обратимости деформации ползучести экспериментальным данным для стали ЭИ736 (500°C) и сплавов: ЭИ437А (700°C), ВЖ98 (900°C) и ЭП693 (700°C). Наблюдается соответствие расчетных и опытных данных. Экспериментальные данные принадлежат соответствующим доверительным интервалам для деформации ползучести, что свидетельствует о достоверности полученных оценок параметров модели.
5. Разработано не имеющее аналогов программное обеспечение для предложенного численного метода идентификации параметров реологической модели.

Конкурирующие интересы. Конфликта интересов в отношении авторства и публикации этой статьи нет.
Авторский вклад и ответственность. Все авторы в равной степени принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.
Финансирование. Работа поддержана Министерством науки и высшего образования Российской Федерации в рамках государственного задания Самарского государственного технического университета (тема № FSSE-2023-0003).

About the authors

Vladimir P. Radchenko

Samara State Technical University

Email: radchenko.vp@samgtu.ru
ORCID iD: 0000-0003-4168-9660
https://www.mathnet.ru/person38375

Dr. Phys. & Math. Sci., Professor; Head of Dept.; Dept. of Applied Mathematics & Computer Science

Russian Federation, 443100, Samara, Molodogvardeyskaya st., 244

Vladimir E. Zoteev

Samara State Technical University

Email: zoteev.ve@samgtu.ru
ORCID iD: 0000-0001-7114-4894
https://www.mathnet.ru/person38585

Dr. Tech. Sci.; Professor; Dept. of Applied Mathematics & Computer Science

Russian Federation, 443100, Samara, Molodogvardeyskaya st., 244

Elena A. Afanasieva

Samara State Technical University

Author for correspondence.
Email: afanasieva.ea@samgtu.ru
ORCID iD: 0000-0001-7815-2723
https://www.mathnet.ru/person188683

Postgraduate Student; Dept. of Applied Mathematics & Computer Science

Russian Federation, 443100, Samara, Molodogvardeyskaya st., 244

References

Supplementary files