ThesisFull.tex

%%% Fedor Garbuzov
%%% fedor.garbuzov@gmail.com
%%% 2019

%!TeX TS-program = xelatex
%!TeX encoding = UTF-8
% !TeX spellcheck = russian-aot
\documentclass[12pt, a4paper]{report}

% Basic packages XeLaTex
\usepackage{xecyr}
\usepackage[english,russian]{babel} 
\usepackage{fontspec} 
\defaultfontfeatures{Ligatures={TeX},Renderer=Basic}
\setmainfont[Ligatures={TeX,Historic}]{Times New Roman} 
\setsansfont{Comic Sans MS}
\setmonofont{Courier New}

% Math
\usepackage{amsmath,euscript,mathrsfs}
\usepackage{amssymb,amsfonts,latexsym,mathtools}

\usepackage{caption,tabularx}

% Graphics
\usepackage{graphicx,setspace,subcaption}
\usepackage[usenames,dvipsnames,svgnames]{xcolor}

% Formatting
\usepackage[top=20mm,bottom=20mm,left=30mm,right=10mm]{geometry}
\setlength{\parindent}{1.25cm}
\onehalfspacing
\usepackage{enumitem}

% Hyperlinks
\usepackage[unicode]{hyperref}
\hypersetup{
	colorlinks=true,
	linkcolor={black!50!black},
	citecolor={blue!50!black},
	urlcolor={blue!50!black},
}

%\mathtoolsset{showonlyrefs=true}

% Dots in section titles
\renewcommand{\thechapter}{\arabic{chapter}.}
\renewcommand{\thesection}{\thechapter\arabic{section}.}
\renewcommand{\thesubsection}{\thesection\arabic{subsection}.}
\renewcommand{\theequation}{\arabic{chapter}.\arabic{equation}}
\renewcommand{\thefigure}{\arabic{chapter}.\arabic{figure}}

% Rename some sections
\addto\captionsrussian{\def\contentsname{Содержание}}
\addto\captionsrussian{\def\bibname{Список источников}}

%\renewcommand\contentsname{Содержание}
%\renewcommand\refname{Список источников}
%\renewcommand{\printtoctitle}[1]{\Huge\sffamily #1}

% Smaller chapter title
\usepackage[raggedright]{titlesec}
%\titleformat{\chapter}[display]
%{\bfseries\Large}%
%{\chaptertitlename\hspace{1ex}\thechapter.}%
%{2 pt}%
%{\bfseries\Large}%
%\titleformat*{\section}{\thesection \Large}
%\titleformat{\section}{}{}{\normalfont\bfseries}{\thesection}
%\titleformat{\section}{\bfseries\large}{\thesection.}{.5em}{}
%\titleformat{\subsection}{\bfseries}{\thesubsection.}{.5em}{}

%\titlespacing{\chapter}{0pt}{0pt}{40pt}

% Dots in table of contents
\usepackage{tocloft}
\renewcommand{\cftchapleader}{\cftdotfill{\cftdotsep}} % for chapters

% Graphics path for figures
\graphicspath{{Figures/}}

% Absolute value declaration
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\norm}{\|}{\|}
\DeclareMathOperator{\erf}{erf}
\DeclareMathOperator{\trace}{tr}
%\DeclareMathOperator{\tg}{tg}
%\DeclareMathOperator{\ctg}{ctg}
\DeclareMathOperator{\sech}{sech}

% Page breaking in multi-line formulae
\allowdisplaybreaks[1]

% Delimiters
\newcommand{\lb}{\left (}
\newcommand{\rb}{\right )}
\newcommand{\lset}{\left \{}
\newcommand{\rset}{\right \}}
\newcommand{\lsq}{\left [}
\newcommand{\rsq}{\right ]}

% Tensors
\newcommand{\vect}[1]{\underline{#1}}
\newcommand{\tens}[1]{\underline{\underline{#1}}}

% Additional commands
\newcommand{\divg}{\text{div}}

% Big 'O' notation
\renewcommand{\O}[1]{O \lb #1 \rb}

% Equation specific commands
\newcommand{\eqtext}[1]{\quad \text{#1} \quad}
\newcommand{\RA}{\quad \Rightarrow \quad}

% Derivatives (normal and partial)
\newcommand{\dd}[1]{\; \mathrm{d} #1}
\newcommand{\diff}[2]{\frac{\mathrm{d} #1}{\mathrm{d} #2}}
\newcommand{\diffn}[3]{\dfrac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
\newcommand{\pdd}[1]{\; \partial #1}
\newcommand{\pdiff}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdiffn}[3]{\dfrac{\partial^{#1} #2}{\partial #3^{#1}}}


\begin{document}

\begin{titlepage}
	\newgeometry{top=20mm,bottom=15mm,left=25mm,right=15mm}
	\begin{center}
		Министерство науки и высшего образования Российской Федерации\\
		Санкт-Петербургский политехнический университет Петра Великого\\
		Институт прикладной математики и механики
	\end{center}
	
	\vspace{5mm}
	
	\begin{flushleft}
		\rule{10cm}{0pt} {Работа допущена к защите}\\
		\rule{10cm}{0pt} Заведующий кафедрой\\
		\rule{10cm}{0pt} <<Прикладная математика>>\\
		\vspace{4mm}
		\rule{10cm}{0pt} \rule{3.6cm}{0.5pt} М.\,Е.~Фролов\\
		\vspace{2mm}
		\rule{10cm}{0pt} "\rule{0.8cm}{0.5pt}" \rule{3.7cm}{0.5pt} 2019 г.
	\end{flushleft}
	
	\vspace{22mm}
	
	\begin{center}
		{\bf ВЫПУСКНАЯ КВАЛИФИКАЦИОННАЯ РАБОТА МАГИСТРА}\\
		\vspace{5mm} 
		{\bf ПРОДОЛЬНЫЕ ВОЛНЫ ДЕФОРМАЦИИ В НЕЛИНЕЙНО УПРУГИХ ВОЛНОВОДАХ}\\
		\vspace{5mm}
		по направлению 01.04.02 «Прикладная математика и информатика»\\
		по образовательной программе 01.04.02\_01 <<Математическое моделирование>>
	\end{center}

	\vspace{20mm}
	\begin{flushleft}
	\begin{tabularx}{\linewidth}{Xr}
		Выполнил & \\
		студент гр. 23641/1 & Ф.\,E.~Гарбузов\\ 
		\vspace{3mm}
		Руководитель & \\ 
		проф., д.т.н. & Б.\,С.~Григорьев\\ 
		\vspace{3mm}
		Научный консультант  & \\ 
		к.ф.-м.н. & Я.\,М.~Бельтюков
	\end{tabularx} 
	\end{flushleft}
	
	\vspace{30mm}
	
	\begin{center}
		Санкт-Петербург\\2019
	\end{center}
	
\end{titlepage}



\newpage
\thispagestyle{empty}
\begin{center}
	\textbf{РЕФЕРАТ}
\end{center}

На 45 с., 13 рис., 1 табл.

\vspace{3mm}
\hspace{-\parindent}НЕЛИНЕЙНАЯ ТЕОРИЯ УПРУГОСТИ, НЕЛИНЕЙНЫЕ ВОЛНЫ, ОБЪЁМНЫЕ ВОЛНЫ ДЕФОРМАЦИИ, СОЛИТОНЫ, МОДЕЛЬ ТИПА БУССИНЕСКА, ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ, ПСЕВДОСПЕКТРАЛЬНЫЙ МЕТОД
\vspace{3mm}

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

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

\vspace{2cm}
\begin{center}
	\textbf{ABSTRACT}
\end{center}

45 pages, 13 pictures, 1 table
\hyphenation{BOUS-SINESQ}

\vspace{3mm}
\hspace{-\parindent}NONLINEAR ELASTICITY, NONLINEAR WAVES, BULK STRAIN WAVES, SOLITONS, \\BOUSSINESQ-TYPE MODEL, NUMERICAL SIMULATION, PSEUDOSPECTRAL METHOD
\vspace{3mm}

Two new Boussinesq-type models describing long nonlinear longitudinal strain waves in elastic rods of circular cross section derived. These models extended to include axially symmetric loading on the lateral surface and longitudinal pre-stretch. The dispersive properties and solitary wave solutions of these models compared with the ones of existing models.

A pseudospectral method applied to perform numerical simulation of the full three-dimensional equations of motion in nonlinearly elastic rod. A series of numerical experiments is conducted in order to excite and compare solitons of Boussinesq-type model and the full model.


\newpage
\setcounter{page}{4}
\tableofcontents



\chapter*{Введение}
\addcontentsline{toc}{chapter}{Введение}

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

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

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

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

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

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


\chapter{Обзор литературы и постановка задачи}

\section{Нелинейные волновые уравнения и методы их решения}
Изучение нелинейных волновых процессов является важной научной задачей, берущей начало ещё в XIX веке в связи с исследованием волн, возникающих на поверхности воды. С тех пор нелинейные волны были обнаружены во многих других физических системах, имеющих самую различную природу, а исследования в этой области ведутся и по сей день. Для описания нелинейных волновых явлений была сформирована единая теория, систематическому изложению которой посвящено множество книг, например,~\cite{Dodd,Ablowitz2011}. Приведём здесь некоторые важные положения, необходимые для дальнейшего изложения.

Исследование нелинейных волн началось в связи с открытием шотландского инженера Дж. Скотта Рассела, наблюдавшего в 1834 году на поверхности канала уединённую волну, бежавшую несколько миль почти без затухания. Впоследствии Рассел неоднократно воспроизводил это явление в экспериментах и определил соотношение между глубиной канала, скоростью и амплитудой этой волны, однако первая математическая модель, описывающая уединённые волны, была получена в работе Ж.~Буссинеска в 1872 году.
Уравнение Буссинеска включает в себя помимо линейных волновых членов также нелинейное и дисперсионное слагаемые:
\begin{equation}\label{1_boussinesq}
u_{tt} - u_{xx} - 3 (u^2)_{xx} - u_{xxxx} = 0.
\end{equation}
Здесь $x$~--~безразмерная пространственная переменная, $t$~--~безразмерное время, $u$~--~нормированная высота поднятия воды над уровнем невозмущённой поверхности, а нижний индекс означает частную производную по соответствующей переменной.
Уравнение \eqref{1_boussinesq} имеет однопараметрическое решение в виде двух уединённых волн, бегущих в противоположных направлениях:
\begin{equation}\label{1_soliton_bouss}
u(x,t) = A \cosh^{-2}\lsq\sqrt{\frac{A}{2}}\lb x \pm t\sqrt{1 + 2A} \rb\rsq,
\end{equation}
где амплитуда $A$ является свободным параметром.
Для той же задачи Д.~Кортевегом и Г.~де~Фризом в 1895 году было выведено другое уравнение:
\begin{equation}\label{1_kdv}
u_{t} + 6 u u_x + u_{xxx} = 0,
\end{equation}
которое тоже имеет решение в виде уединённой волны.

Существование решений нелинейных уравнений в виде устойчивых уединённых волн обусловлено так называемым балансом нелинейности и дисперсии. Нелинейное слагаемое, нестрого говоря, стремится сделать фронт волны более крутым и в конечном счёте опрокинуть его, а дисперсионное слагаемое, наоборот, стремится сделать волну более пологой. 
Стоит отметить, что в реальных физических системах как правило присутствует диссипация энергии. Простейшую нелинейную систему с дисперсией и диссипацией описывает уравнение Бюргерса-Кортевега-де~Фриза:
\begin{equation}\label{1_burgers-kdv}
u_{t} + 6u u_x + u_{xxx} = \nu u_{xx}.
\end{equation}
Существует множество других классических уравнений теории нелинейных волн, например, нелинейное уравнение Шрёдингера, однако в настоящей работе они не затрагиваются. Отметим, что приведённые выше уравнения имеют множество других похожих на себя уравнений. Так, уравнением \emph{типа Буссинеска} будем называть всякое уравнение вида \eqref{1_boussinesq}, возможно, с другим нелинейным или дисперсионным слагаемым, например:
$$
u_{tt} - u_{xx} - 3 (u^3)_{xx} + u_{xxtt} = 0.
$$

Важнейшие свойства уединённых волн, являющихся решением уравнения Кортевега-де Фриза \eqref{1_kdv}, были открыты в 1965 году Н.~Забужским и М.~Краскалом~\cite{Zabusky}. Оказалось, что уединённые волны сталкиваются <<упруго>>, то есть после взаимодействия полностью восстанавливают свою форму.
На рисунке \ref{fig:kdv_collision} изображены результаты численного эксперимента, где видно, как волна с большей амплитудой обгоняет волну с меньшей амплитудой, причём в результате столкновения уединённые волны не изменили свою форму, а лишь претерпели фазовый сдвиг.
Для того, чтобы подчеркнуть <<упругий>> характер взаимодействия, свойственный частицам, Забужский и Краскал назвали такие уединённые волны \emph{солитонами}. 
\begin{figure}[h]
	\centering
	\includegraphics[width=0.6\linewidth]{kdv_collision.png}
	\caption{Попутное столкновение двух солитонов (из книги Р. Додда~\cite{Dodd}).}
	\label{fig:kdv_collision}
\end{figure}

Работа Забужского и Краскала дала толчок к аналитическим исследованиям уравнения Кортевега-де Фриза (КдФ), которые привели к возникновению в 1967 году метода обратной задачи рассеяния (МОЗР), позволяющего получить точное аналитическое решение уравнения по заданному начальному условию~\cite{GardnerIST}. Позже этот метод был обобщён и на некоторые другие нелинейные волновые уравнения, однако лишь очень небольшое количество нелинейных систем являются полностью интегрируемыми (решаемыми с помощью МОЗР)~\cite{Zakharov}.

Для численного эксперимента в работе 1965 года Забужский и Краскал использовали простейшую явную конечно-разностную схему второго порядка по временной и пространственной переменной. Впоследствии был предложен ряд других более сложных неявных конечно-разностных схем и псевдоспектральных методов, обзору которых посвящена работа~\cite{Taha}. Численное моделирование особенно важно для неинтегрируемых систем, к числу которых относятся, например, все уравнения типа Буссинеска, кроме \eqref{1_boussinesq}. Для таких уравнений были предложены конечно-разностные схемы, в которых на каждом шаге для обработки нелинейного слагаемого проводится несколько внутренних итераций~\cite{Christov, Kolkovska}.
Помимо метода конечных разностей к нелинейным уравнениям применяют и другие численные методы, например, конечных объёмов и конечных элементов~\cite{Dutykh_Boussinesq, Karczewska}. Отметим, что для численного решения уравнений, обладающих гладкими решениями, в областях простой геометрии одним из лучших методов является псевдоспектральный метод, широко применяющийся в нелинейных задачах~\cite{Gottlieb_Orszag, Canuto2007}.




\section{Нелинейная динамика твёрдого тела}
В этом разделе приведены важные для дальнейшего изложения сведения по нелинейной динамике твёрдого тела, систематическому описанию которой посвящено множество книг, например,~\cite{LurieNL}. 

Динамика упругой сплошной среды, занимающей объём $\Omega$, описывается уравнениями движения, которые в векторном виде в случае однородного тела представляются следующим образом:
\begin{equation}\label{1_dynamics}
\rho\ddot{\vect{U}}(\vect{x}, t) = \divg \tens{P} + \vect{F}, \quad \vect{x}\in\Omega.
\end{equation}
Здесь $\rho$ -- плотность материала, $\vect{U}(\vect{x},t)$ -- вектор перемещений, $\vect{x}$ -- координаты точки среды в отсчётной конфигурации, $\tens{P}$ -- первый тензор напряжений Пиолы-Кирхгофа, $\vect{F}$ -- плотность массовых сил, точка обозначает частную производную по времени, а дивергенция берётся по координатам в отсчётной конфигурации. Тензор напряжений $\tens{P}$ выражается через тензор деформации $\tens{E}$ следующим образом:
\begin{equation} \label{1_piola}
\tens{P} = (\tens{I} + \nabla\vect{U}) \cdot \pdiff{W}{\tens{E}},
\end{equation}
где $W$ -- плотность энергии деформации, а тензор деформации связан с градиентом перемещения:
\begin{equation}\label{1_strain}
2\tens{E} = (\nabla\vect{U})^T + \nabla\vect{U} + (\nabla\vect{U})^T\cdot\nabla\vect{U} \,.
\end{equation}
Заметим, что в линейной теории деформация предполагается бесконечно малой и нелинейное слагаемое в \eqref{1_strain} отбрасывается. Для завершения постановки задачи уравнения \eqref{1_dynamics} -- \eqref{1_strain} необходимо дополнить соотношением, связывающем энергию и деформацию, а также граничными условиями:
\begin{equation}\label{1_bc}
\vect{U} = \vect{U}_b, \;\; \vect{x}\in S_U; \qquad \tens{P}\cdot\vect{n} = \vect{P}_b, \;\; \vect{x}\in S_P; \qquad S_U\cup S_P = \partial\Omega.
\end{equation}
%Случай, когда $\vect{P}_b = 0$ будем называть условием свободной поверхности.

Энергия деформации $W$ однородного и изотропного тела может быть разложена в ряд по инвариантам $I_i$ тензора деформации:
\begin{equation}\label{1_murnaghan}
W = \frac{\lambda + 2\mu}{2}I_1(\tens{E})^2 - 2\mu I_2(\tens{E}) + \frac{l+2m}{3}I_1(\tens{E})^3 - 2m I_1(\tens{E}) I_2(\tens{E}) + n I_3(\tens{E}) + \dots,
\end{equation}
при этом коэффициенты в этом разложении характеризуют упругость материала и называются модулями упругости ($\lambda$ и $\mu$ --- коэффициенты Ламе, а $l$, $m$, $n$ --- модули Мурнагана). Заметим, что первые два слагаемых в приведённом разложении являются слагаемыми второго порядка относительно компонент тензора $\tens{E}$, а следующие три --- третьего порядка.
В разложении \eqref{1_murnaghan} для линейно упругого материала удерживаются только слагаемые второго порядка, а для слабо нелинейного материала Мурнагана~\cite{Murnaghan} учитываются ещё и слагаемые третьего порядка. Существуют другие нелинейно упругие материалы, например, материал Муни-Ривлина или Огдена, однако они предназначены в первую очередь для описания резино- и пеноподобных материалов, подверженных большим деформациям~\cite{Bergstrom}. Отметим, что нелинейно упругие материалы иногда называют гиперупругими.

Помимо классической постановки задачи в виде дифференциальных уравнений в частных производных \eqref{1_dynamics}, \eqref{1_bc}, существует вариационная постановка на основе принципа Гамильтона, гласящего, что истинная траектория системы $\vect{U}$ является стационарной точкой функционала действия $\mathcal{S}$:
\begin{equation}\label{1_action}
\delta\mathcal{S} = \delta\int_{t_1}^{t_2} dt \lsq \int_{\Omega} \lb \frac12\rho\dot{\vect{U}}^2 - W + \vect{F}\cdot\vect{U} \rb dx + \int_{S_P} \vect{P}_b\cdot\vect{U} ds \rsq = 0.
\end{equation}
В \eqref{1_action} варьирование происходит по перемещениям $\vect{U}$. Отметим, что существует обобщённый принцип Гамильтона, где в функционал действия включаются соотношения \eqref{1_piola} и \eqref{1_strain}, а варьирование осуществляется не только по перемещениям, но и по деформациям $\tens{E}$ и напряжениям $\tens{P}$~\cite{Yu}.
%Отметим, что некоторые авторы использовали <<неполный>> принцип наименьшего действия, при котором в функционал действия включаются только уравнения движения (объёмный интеграл по $\Omega$), а граничные условия ставятся в виде уравнений \eqref{1_bc}.



\section{Нелинейные волны деформации в твёрдых упругих волноводах}
Изучение нелинейных волн деформации в твёрдых телах, в том числе солитонов деформации, является важной темой современного изучения волн~\cite{S_book, P_book}.
Разработка теории началась в 1970\babelhyphen{nobreak}х годах с исследования волн в упругом стержне круглого сечения, поскольку такая геометрия волновода является наиболее простой.

Исторически первым исследованием стала работа Г.~Нариболи и А.~Седова, которым удалось вывести уравнение Бюргерса-Кортевега-де~Фриза для длинных продольных волн в бесконечном вязкоупругом осесимметричном стержне со свободной от напряжений поверхностью~\cite{NS}. Для этого уравнения нелинейной теории упругости~\eqref{1_dynamics} и граничные условия~\eqref{1_bc}, записанные в цилиндрической системе координат $(x,r,\varphi)$, были упрощены с помощью:
\begin{itemize}[noitemsep,topsep=1pt]
	\item предположения о малости радиуса стержня $a \ll 1$,
	\item разложения перемещений в степенной ряд по радиусу стержня: 
	\begin{align}\label{1_nariboli}
	U(x,r,t) &= U_0(x,t) + a^2 U_2(x,r,t) + \mathcal{O}(a^4),\\
	V(x,r,t) &= -a \nu r \pdiff{U_0}{x} + a^3 V_3(x,r,t) + \mathcal{O}(a^5),
	\end{align}
	где $U$ --- продольное перемещение вдоль оси стержня, совпадающей с осью $x$, $V$ --- радиальное перемещение, а $\nu$ --- коэффициент Пуассона,
	\item предположения о малых деформациях $U, \, V \sim \varepsilon \ll 1$.
\end{itemize}
Позже Л.~Островский и А.~Сутин получили модель типа Буссинеска, %:
%\begin{equation}\label{1_ostrovsky_bou}
%u_{tt} - u_{xx} - 3 (u^2)_{xx} - u_{ttxx} = 0,
%\end{equation}
используя принцип Гамильтона и нижеследующие гипотезы, позволившие упростить функционал действия задачи~\cite{OS}: 
\begin{equation}\label{1_ostrovsky_hyp}
U(x,r,t) = U(x,t), \qquad V(x,r,t) = -\nu r \pdiff U x.
\end{equation}
Первая из этих гипотез называется гипотезой плоских сечений и означает, что поперечные сечения стержня остаются плоскими после деформации, а вторая гипотеза аналогична гипотезе Кирхгофа-Лява в теории тонких пластин и оболочек~\cite{Love}.
А.\,М.~Самсонов, используя подход Островского и Сутина, предложил модель типа Буссинеска с двумя дисперсионными слагаемыми и обобщил её на случай с меняющимися вдоль оси стержня радиусом и модулями упругости~\cite{S_book}. Коэффициенты модели Самсонова с двумя типами дисперсионных членов были позже уточнены в работах А.\,М.~Самсонова и А.\,В.~Порубова~\cite{SP}.
А.\,В.~Порубовым и М.~Веларде предложена дисперсионно-диссипативная модель для длинных волн в упругом стержне, помещённом в вязкоупругую среду~\cite{PV}.
Модель типа Буссинеска с тремя типами дисперсионных членов обсуждалась В.\,И.~Ерофеевым, однако коэффициент при нелинейном слагаемом в его модели отличается от соответствующего коэффициента у Островского и Самсонова~\cite{E_book}. Все выводы моделей типа Буссинеска в упомянутых исследованиях основывались на представлении Мурнагана для энергии упругой деформации и последующем упрощении полного функционала действия задачи с использованием некоторых гипотез. 

Несколько другой подход к задаче применён в работе X. Дая и X. Фана, которым удалось упростить полные уравнения движения с граничными условиями в виде свободной от напряжений поверхности стержня, сведя их к системе из двух связанных уравнений~\cite{DF}. Для этого была введена система масштабов для переменных и функций так, что масштаб перемещений $h$ и радиус стержня $a$предполагались малыми по сравнению с характерной длиной волны~$l$: $\varepsilon = h/l \ll 1$, $\delta = a^2/l^2 \ll 1$.
Полные уравнения были упрощены при помощи разложения перемещений в степенной ряд по радиальной координате и отбрасывания членов порядка $O(\varepsilon^2, \varepsilon\delta, \delta^2)$.
В другой работе X.~Дай и З.~Цай применили аналогичный асимптотический вывод для описания волн в предварительно растянутом гиперупругом стержне, сделанном из материала Муни-Ривлина~\cite{DC}.

Во всех приведённых выше работах твёрдое тело считалось непрерывным. Однако помимо непрерывной модели существует решёточная (дискретная) модель, согласно которой твёрдое тело представляется в виде системы часиц некоторой массы, соединённых пружинами. В рамках такой модели К.\,Р. Хуснутдинова и др., предполагая пружины нелинейно упругими, получили систему разностно-дифференциальных уравнений, которая в континуальном пределе сводится к уравнению типа Буссинеска~\cite{KSZ}. Заметим, что уравнение модели было получено из полных уравнений движения с помощью асимптотических методов без использования упрощающих гипотез. Интересно, что в этом исследовании была выведена модель типа Буссинеска с тремя дисперсионными слагаемыми, а также система связанных уравнений типа Буссинеска для волн в слоистом волноводе с неидеальным контактом.
В недавних исследованиях модели типа Буссинеска использовались для изучения распространения длинных продольных уединённых волн деформации в волноводе с расслоением \cite{KS, KT1, KT2}, а некоторые соответствующие экспериментальные наблюдения были опубликованы в \cite{JAP2010, JAP2012}. 
%Отметим исследования по распространению нелинейных волн на дефектах в решетках \cite{ML} и струнах \cite{AM}. Некоторые другие модели типа Буссинеска были получены для описания распространения поперечных волн большой амплитуды~\cite{DS}.

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

Вывод моделей в настоящей работе выполнен с помощью пакета символьных вычислений Mathematica, а результаты работы частично опубликованы автором в сотрудничестве с К.\,Р.~Хуснутдиновой и И.\,В.~Семёновой~\cite{Garbuzov}.




\chapter{Слабо нелинейные продольные волны деформации в тонких волноводах}

\section{Формулировка задачи}
Рассмотрим стержень круглого поперечного сечения радиуса $R$. Введём цилиндрическую систему координат $(x, r, \varphi)$, где $x$ --- осевая координата, $r$ --- продольная, $\varphi$ --- угловая, как показано на рисунке \ref{fig:rod}. Положим стержень бесконечным вдоль оси $x$. Используя Лагранжев подход, введём вектор перемещения точек тела: $\vect{U} = (U, V, W)$, где $U$ --- осевое (продольное) перемещение, $V$ --- радиальное (поперечное) перемещение, а $W$ --- вращение.
\begin{figure}[h]
	\centering
	\includegraphics[width=0.33\textwidth]{1_RodSchematic}
	\caption{Стержень с круглым поперечным сечением.}
	\label{fig:rod}
\end{figure}

Следуя предыдущим исследованиям, которые обсуждались в главе 1, будем рассматривать стержень, сделанный из материала Мурнагана, энергия упругой деформации которого представляется в виде:
\begin{equation}\label{2_murnaghan_pot}
W = \frac{\lambda + 2\mu}{2}I_1(\tens{E})^2 - 2\mu I_2(\tens{E}) + \frac{l+2m}{3}I_1(\tens{E})^3 - 2m I_1(\tens{E}) I_2(\tens{E}) + n I_3(\tens{E}),
\end{equation}
где $I_1(\tens{E}) = \trace \tens{E},\  I_2(\tens{E}) = \lsq(\trace \tens{E})^2 - \trace \tens{E}^2\rsq/\,2,\ I_3(\tens{E}) = \trace \tens{E}$ являются инвариантами тензора деформации Грина
$\tens{E} = \lb (\nabla\vect{U})^T + \nabla\vect{U} + (\nabla\vect{U})^T\cdot\nabla\vect{U} \rb / \, 2$; 
$\lambda$,~$\mu$ --- коэффициенты Ламе; $l$,~$m$,~$n$ --- модули Мурнагана. Здесь и далее в тексте работы все частные производные берутся по координатам в отсчётной конфигурации. Отметим, что модель Мурнагана является общей моделью слабо нелинейных упругих деформаций.

Рассмотрим задачу, в которой отсутствует кручение стержня, а продольное и поперечное перемещения $U$ и $V$ не зависят от угла $\varphi$:
\begin{equation}\label{2_assumptions}
U = U(x,r,t), \quad V = V(x,r,t), \quad W = 0.
\end{equation}
Уравнения движения, в условиях \eqref{2_assumptions} и отсутствия массовых сил, принимают вид
\begin{align}
\label{2_eq1_0}
&\rho  \frac{\partial^2 U(x,r,t)}{\partial t^2}-\frac{\partial P_{x x}}{\partial x}-\frac{\partial P_{xr}}{\partial r}-\frac{P_{xr}}{r} = 0,\\
\label{2_eq2_0}
&\rho  \frac{\partial^2 V(x,r,t)}{\partial t^2} - \frac{\partial P_{rx}}{\partial x}-\frac{\partial P_{rr}}{\partial r}-\frac{P_{rr} - P_{\varphi\varphi}}{r} = 0,
\end{align}
а третье уравнение представляется в виде тождества $0\equiv 0$. Здесь $P_{\alpha \beta}$ обозначает компоненту первого тензора Пиолы-Кирхгофа, а $\rho$ --- плотность материала.

Зададим на поверхности стержня осесимметричное напряжение $\vect{P_b} = (T(x,t), P(x,t), 0)$. Тогда граничные условия имеют вид:
\begin{align}
P_{rr} &= P(x, t) \quad \mbox{при} \quad r = R \label{2_bc_rr},\\
P_{xr} &= T(x, t) \quad \mbox{при} \quad r = R \label{2_bc_rx}.
\end{align}
Поскольку компонента $ P_{\varphi r} \equiv 0 $, третье граничное условие $P_{\varphi r} = 0$ при $r = R$ выполняется автоматически.



\section{Вывод уравнения типа Буссинеска с внешним воздействием}
\label{sec:derivation}
\subsection{Вывод с помощью степенных разложений по радиусу}
Подход к выводу уравнения модели в этом разделе опирается на метод, описанный в~\cite{DF}. Упростим этот метод с помощью разложений, использованных для вывода линейной модели~\cite{bostrm2000}. Таким образом, будем искать решение в виде степенного ряда по радиальной координате:
\begin{align}
\label{u_series}
U(x,r,t) &= U_0(x,t) + r^2 U_2(x,t) + r^4 U_4(x,t) + \dots \, ,\\
\label{v_series}
V(x,r,t) &= r V_1(x,t) + r^3 V_3(x,t) + r^5 V_5(x,t) + \dots \, .
\end{align}
Здесь продольное перемещение разложено в ряд по чётным степеням радиуса, в то время как поперечное перемещение по нечётным, что обусловлено рассмотрением именно продольных колебаний. 
Отметим, что в отличие от~\cite{DF}, мы сведём задачу к одному уравнению типа Буссинеска, учтём напряжение, приложенное к поверхности стержня, а также рассмотрим предварительно растянутый стержень.

Введём масштабные множители, выделяющие среди прочих задачу о распространении длинных по сравнению с радиусом стержня волн малой амплитуды. Тогда безразмерные переменные и функции определяются следующими выражениями:
\begin{equation} \label{scales1}
\tilde t = \frac{t}{L/c}, \quad \tilde x = \frac{x}{L}, \quad \tilde r = \frac{r}{L}, \quad \tilde U = \frac{U}{\varepsilon L}, \quad \tilde V = \frac{V}{\varepsilon L}, \quad \tilde P = \frac{P}{E \varepsilon}, \quad \tilde T = \frac{T}{E \varepsilon\delta},
\end{equation}
из которых следует, что
$\displaystyle \tilde U_n =  \frac{L^n U_n}{\varepsilon L},  \ \tilde V_n =  \frac{L^n V_n}{\varepsilon L}$ для $n \geqslant 0$.
Здесь $L$ является характерной длиной волны, $c$ --- скорость линейной волны, $E$ --- модуль Юнга, $\varepsilon \ll 1$ --- малый параметр амплитуды, $\displaystyle \delta = \frac{R}{L} \ll 1$ --- второй малый параметр, а тильда обозначает безразмерную величину. В дальнейшем мы будем использовать выражения для модуля Юнга и коэффициента Пуассона через коэффициенты Ламе:
\begin{equation}\label{young_mod}
E = \frac{\mu(3\lambda + 2\mu)}{\lambda + \mu}, \quad \nu = \frac{\lambda}{2 (\lambda + \mu)}.
\end{equation}
С учётом \eqref{scales1} разложения \eqref{u_series}, \eqref{v_series} представляются в виде:
\begin{align}
\label{u_series_scaled}
\widetilde U(\tilde x, \tilde r, \tilde t) &= \widetilde{U}_0(\tilde x, \tilde t) + \tilde{r}^2\widetilde{U}_2(\tilde x, \tilde t) + \tilde{r}^4\widetilde{U}_4(\tilde x, \tilde t) + O(\tilde r^6),\\
\label{v_series_scaled}
\widetilde V(\tilde x, \tilde r, \tilde t) &= \tilde{r} \widetilde{V}_1(\tilde x, \tilde t) + \tilde{r}^3\widetilde{V}_3(\tilde x, \tilde t) + \tilde{r}^5\widetilde{V_5}(\tilde x, \tilde t) + O(\tilde r^7).
\end{align}
Радиальная координата $\tilde r$ точек стержня принимает значения от $0$ до $\delta$ и, следовательно, является малой величиной.
В дальнейшем мы опустим тильду над безразмерными величинами.

Подставляя \eqref{u_series_scaled} и \eqref{v_series_scaled} в уравнения движения \eqref{2_eq1_0}, \eqref{2_eq2_0}, получаем
\begin{equation} 
\label{eq1_1}
\begin{split}
	&\rho c^2 U_{0tt} - (\lambda + 2\mu) U_{0xx} - 2(\lambda + \mu) V_{1x} - 4\mu U_2 + \Phi_1(U_0, V_1, U_2) \varepsilon\\
	&\qquad+ \left[\rho c^2 U_{2tt} - (\lambda + 2\mu)U_{2xx} - 4(\lambda + \mu)V_{3x} - 16\mu U_4\right] r^2 + O(\varepsilon^2, \varepsilon r^2, r^4) = 0,
\end{split}
\end{equation}
\begin{equation} \label{eq2_1}
\begin{split}
	&r \big( \rho c^2 V_{1tt} - \mu V_{1xx} - 2(\lambda + \mu)U_{2x} - 8(\lambda + 2\mu)V_3 - \Phi_2(U_0, V_1, U_2, V_3)\varepsilon \\
	&\qquad- \left[\rho c^2 V_{3tt} - \mu V_{3xx} - 4(\lambda + \mu)U_{4x} - 24(\lambda + 2\mu)V_5 \right] r^2  + O(\varepsilon^2, \varepsilon r^2, r^4)\big) = 0,
\end{split}
\end{equation}
где индексы $x$ и $t$ обозначают частные производные по соответствующим переменным, а выражения для нелинейных функций $\Phi_1$ и $\Phi_2$ приведены в Приложении~1.
Функции $ U_2 $, $ V_3 $, $ U_4 $ могут быть получены из \eqref{eq1_1} и \eqref{eq2_1}, приравнивая к нулю коэффициенты при различных степенях~$r$:
\begin{align}
\label{U2}
U_2 &= \frac{1}{4\mu} \left[ \rho c^2 U_{0tt} - (\lambda + 2\mu) U_{0xx} - 2(\lambda + \mu) V_{1x} \right] + \varepsilon f_2(x,t) + O(\varepsilon^2),\\
\label{V3}
V_3 &= \frac{1}{8(\lambda + 2\mu)} \left[ \rho c^2 V_{1tt} - 2(\lambda + \mu) U_{2x} - \mu V_{1xx} \right] + \varepsilon f_3(x,t) + O(\varepsilon^2),\\
\label{U4}
U_4 &= \frac{1}{16\mu}\left[\rho c^2 U_{2tt} - (\lambda + 2\mu) U_{2xx} - 4(\lambda + \mu) V_{3x}\right] + O(\varepsilon).
%V_5 &=& \frac{1}{24(\lambda + 2\mu)} \left(\rho c^2 V_{3tt} - 4(\lambda+\mu)U_{4x} - \mu V_{3xx}\right) + O(\varepsilon).
\end{align}
Выражения для функций $f_2$ и $f_3$ приведены в Приложении 1.

Затем, подставляя функции $ U_2 $, $ V_3 $, $ U_4 $ в граничные условия \eqref{2_bc_rr}, \eqref{2_bc_rx}, которые в безразмерном виде должны выполняться при $r = \delta$, получаем уравнения:
\begin{equation} \label{2_bc_rr_subst}
\begin{split}
2 (\lambda + \mu) V_1 + \lambda U_{0x} &+\varepsilon \Psi_1(U_0, V_1) + \frac{\delta^2}{8} \bigg[ (\lambda + 3\mu) U_{0xxx}- \frac{\rho c^2(\lambda + 3\mu)}{\lambda + 2\mu} U_{0xtt} \\
&+ \frac{2\rho c^2(2\lambda + 3\mu)}{\lambda + 2\mu} V_{1tt} + 2\lambda V_{1xx}\bigg] 
+ O(\varepsilon^2, \varepsilon\delta^2, \delta^4) =  \frac{\mu(3\lambda + 2\mu)}{\lambda + \mu} P,
\end{split}
\end{equation}
\begin{equation} \label{2_bc_rx_subst}
\begin{split}
\rho  c^2 U_{0tt} -2 \lambda  V_{1x}-(\lambda +2 \mu ) U_{0xx} - \varepsilon \Psi_2(U_0, V_1)
+ \frac{\delta^2}{8}\bigg[(3\lambda + 4\mu)U_{0xxxx} + \frac{\rho^2 c^4}{\mu}U_{0tttt}\\
- \frac{\rho c^2\left(\lambda^2 + 7\lambda\mu + 8\mu^2\right)}{\mu(\lambda + 2\mu)} U_{0xxtt} + 2(3\lambda + 2\mu)V_{1xxx} - \frac{2 \rho  c^2 \left(\lambda ^2+4 \lambda  \mu +2 \mu ^2\right)}{\mu(\lambda + 2\mu)} V_{1xtt} \bigg]\\
+\,O(\varepsilon^2, \varepsilon\delta^2, \delta^4)
= \frac{2\mu(3\lambda + 2\mu)}{\lambda + \mu} T,
\end{split}
\end{equation}
где нелинейные члены выражаются следующим образом: 
\begin{align}
	\nonumber
	\Psi_1 &= (4l + 2m + 3\lambda + 3\mu) V_1^2 + (4l - 2m + n + \lambda) U_{0x} V_1 + \frac{1}{2} (2l + \lambda) U_{0x}^2, \\
	\nonumber
	\Psi_2 &= \left((4l - 2m + n + \lambda) V_1^2 + 2(2l + \lambda) U_{0x} V_1 + \frac12(2l + 4m + 3\lambda + 6\mu) U_{0x}^2 \right)_x.
\end{align}
Отметим, что при $\varepsilon = 0$ уравнения \eqref{2_bc_rr_subst} и \eqref{2_bc_rx_subst} сводятся к уравнениям, полученным в линейной задаче~\cite{bostrm2000}. Эта система связанных уравнений представляет собой довольно сложную модель, однако она может быть сведена к одному уравнению типа Буссинеска.

Существует два естественных способа вывода модели типа Буссинеска. В первом способе исключение функции $V_1$ из уравнений \eqref{2_bc_rr_subst} и \eqref{2_bc_rx_subst} осуществляется с помощью асимптотического выражения, следующего из уравнения \eqref{2_bc_rr_subst}:
\begin{equation} \label{v1_asympt}
V_1(x, t) = \frac{\mu(3\lambda + 2\mu) P}{2(\lambda + \mu)^2} - \frac{\lambda U_{0x}}{2(\lambda + \mu)} + \varepsilon f(x,t) + \delta^2 g(x,t) + O(\varepsilon^2, \varepsilon\delta^2, \delta^4),
\end{equation}
где неизвестные функции $f$ и $g$ ищутся из условия равенства нулю коэффициентов при $\varepsilon$ и $\delta^2$ в \eqref{2_bc_rr_subst}. Их вид представлен в Приложении 1. Затем, подстановка $V_1$ в \eqref{2_bc_rx_subst} приводит к следующему уравнению типа Буссинеска относительно $U_0$:
\begin{equation} \label{eq_u0_asympt}
\begin{split}
\rho c^2 U_{0tt} &- \frac{\mu(3\lambda + 2\mu)}{\lambda + \mu}\left(U_{0xx} + \frac{\lambda}{\lambda + \mu} P_x + 2T\right) - \varepsilon \left(\gamma_1 U_{0x}^2 + \gamma_2 U_x P + \gamma_3 P^2 \right)_x \\
&+ \delta^2 \bigg[\frac{\rho ^2 c^4 U_{0tttt}}{8\mu} + \frac{\mu (3\lambda + 2\mu)^2 U_{0xxxx}}{8(\lambda + \mu)^2} - \frac{\rho c^2 \left(7\lambda^2 + 10\lambda\mu + 4\mu^2\right) U_{0xxtt}}{8(\lambda + \mu)^2} + F \bigg]\\
&\hspace{100mm}+ O(\varepsilon^2, \varepsilon\delta^2, \delta^4) = 0.
\end{split}
\end{equation}
Здесь нелинейные коэффициенты $\gamma_i$ и функция $F$ представляются в виде:
\begin{align} 
	\nonumber
	\gamma_1 &= \frac{3n(\lambda + \mu)\lambda^2 + 2\mu \left[9\lambda^3 + 24\mu\lambda^2 + 21\mu^2\lambda + m(3\lambda + 2\mu)^2 + 2\mu^2 (l + 3\mu)\right]}{4 (\lambda + \mu)^3},\\
	\nonumber
	\gamma_2 &= \frac{\left[3\lambda^3 + 5\lambda^2\mu + 2\lambda\mu^2 + 4l\mu^2 + 2\lambda m(3\lambda + 2\mu) - 2\lambda n (\lambda + \mu)\right] \mu(3\lambda + 2\mu)}{2(\lambda + \mu)^4},\\
	\nonumber
	\gamma_3 &= \frac{\left[n (\lambda +\mu )-2 \left(\lambda ^2+\lambda  \mu -2 l \mu \right)-2 m (2 \lambda +\mu )\right] \mu^2 (3\lambda + 2\mu)^2}{4(\lambda + \mu)^5},\\
	\nonumber
	F &= \frac{3\lambda + 2\mu}{8\mu(\lambda + \mu)^3}\left[ \mu (4\lambda^2 + 5\lambda\mu + 2\mu^2) P_{xxx} - \rho c^2(\lambda^2 + \lambda\mu + \mu^2) P_{xtt}\right].
\end{align}

Другой метод основан на исключении $V_1$ из \eqref{2_bc_rr_subst} и \eqref{2_bc_rx_subst} таким образом, каким это сделано в~\cite{bostrm2000} для линейной задачи. В линейном случае этот подход не использует асимптотическое выражение \eqref{v1_asympt} и приводит к уравнению того же типа, что и \eqref{eq_u0_asympt}, но с другими коэффициентами при дисперсионных слагаемых. Запишем уравнения \eqref{2_bc_rr_subst} и \eqref{2_bc_rx_subst} в виде:
\begin{align}%\label{key}
L_{1} V_1 + \varepsilon\Psi_1(U_0, V_1) &= a_1 P + M_1 U_0 + O(\varepsilon^2, \varepsilon\delta^2, \delta^4),\\
L_{2} V_1 + \varepsilon\Psi_2(U_0, V_1) &= a_2 T + M_2 U_0 + O(\varepsilon^2, \varepsilon\delta^2, \delta^4),
\end{align}
где $a_1$ и $a_2$ --- константы; $L_{1}$, $L_{2}$ и $M_{1}$, $M_{2}$ --- линейные дифференциальные операторы, действующие на $V_1$ и $U_0$ соответственно в уравнениях \eqref{2_bc_rr_subst} и \eqref{2_bc_rx_subst}. Теперь, применяя $L_{2}$ к первому уравнению, $L_{1}$ ко второму и вычитая одно уравнение из другого, получаем:
\begin{equation}
\varepsilon [L_{2}\Psi_1(U_0, V_1) - L_{1}\Psi_2(U_0, V_1)] =  L_{2}\lb a_1 P + M_1 U_0\rb - L_{1}\lb a_2 T + M_2 U_0\rb + O(\varepsilon^2, \varepsilon\delta^2, \delta^4).
\end{equation}
Здесь $V_1$ исключена из линейной части уравнений точно, а не асимптотически. Чтобы исключить её из нелинейной части, воспользуемся выражением \eqref{v1_asympt} и получим следующее уравнение:
\begin{gather} \label{eq_u0_bostr}
\begin{split}
%\nonumber
\rho c^2 U_{0tt} &- \frac{\mu(3\lambda + 2\mu)}{\lambda + \mu} \left(U_{0xx} + \frac{\lambda}{\lambda + \mu} P_x + 2T\right)
- \varepsilon \left(\gamma_1 U_{0x}^2 + \gamma_2 U_x P + \gamma_3 P^2 \right)_x \\
&\hspace{15mm}+ \delta^2 \bigg[\frac{\rho ^2 c^4 (\lambda^2 + 5\lambda\mu + 5\mu^2) U_{0tttt}}{8\mu(\lambda+2\mu)(\lambda+\mu)} - \frac{\rho c^2 \left(6\lambda^2 + 21\lambda \mu + 14\mu^2\right) U_{0xxtt}}{8(\lambda + 2\mu)(\lambda + \mu)} \\
&\hspace{55mm} + \frac{\mu(3\lambda + 2\mu) U_{0xxxx}}{4(\lambda + \mu)} + G\bigg] + O(\varepsilon^2, \varepsilon\delta^2, \delta^4) = 0,
\end{split}\\
\nonumber
G  = \frac{\mu(3\lambda + 2\mu)}{8(\lambda + \mu)^2} \left[(3\lambda + 2\mu) P_{xxx} - \frac{\rho c^2(\lambda^2 + 4\lambda\mu + 2\mu^2)}{\mu(\lambda + 2\mu)} P_{xtt} - \frac{2\rho c^2 (2\lambda + 3\mu)}{\lambda + 2\mu} T_{tt} - 2\lambda T_{xx}\right]. 
\end{gather}
Отметим, что в линейном приближении при $\varepsilon = 0$ уравнение \eqref{eq_u0_bostr} совпадает с уравнениями, выведенными для линейной задачи в~\cite{bostrm2000}.
Из \eqref{eq_u0_asympt} и \eqref{eq_u0_bostr}, задавая $\varepsilon = 0$, $\delta = 0$ и $P = T = 0$, получаем скорость линейной продольной волны в бесконечно тонком стержне:
\begin{equation}
\label{lin_velocity}
c = \ \sqrt{\frac{\mu(3\lambda + 2\mu)}{\rho(\lambda+\mu)}} = \sqrt{\frac{E}{\rho}}.
\end{equation}

Теперь перепишем оба выведенных уравнения типа Буссинеска \eqref{eq_u0_asympt} и \eqref{eq_u0_bostr} в унифицированной форме и выразим коэффициенты Ламе через модуль Юнга $ E $ и коэффициент Пуассона $ \nu $:
\begin{equation}\label{eq_u0_fin}
\begin{split}
U_{0tt} &- U_{0xx} - 2\left(\nu P_{x} + T\right) - \frac{\varepsilon}{2 E} \left(\beta_1U_{0x}^2 + 2 \beta_2 U_{0x} P + \beta_3 P^2 \right)_x \\
&+ \delta^2 \left(\alpha_1^{(i)} U_{0tttt} + \alpha_2^{(i)} U_{0xxtt} + \alpha_3^{(i)} U_{0xxxx} + F^{(i)}\right) + O(\varepsilon^2, \varepsilon\delta^2, \delta^4) = 0, \quad i = 1,2,
\end{split}
\end{equation}
где $i=1$ соответствует уравнению \eqref{eq_u0_asympt}, а $i=2$ уравнению \eqref{eq_u0_bostr}, а коэффициенты принимают следующий вид:
\begin{align}
\label{alpha_1}
&\alpha_1^{(1)} = \alpha_3^{(1)} = \frac{1 + \nu}{4}, \quad \alpha_2^{(1)} = -\frac{1 + \nu + \nu^2}{2}, \\
\label{alpha_2}
&\alpha_1^{(2)} = \frac{5 - 5\nu - 6\nu^2 + 4\nu^3}{8(1-\nu)},\quad \alpha_2^{(2)} = -\frac{7 - 7\nu - 2\nu^2}{8(1-\nu)}, \quad \alpha_3^{(2)} = \frac14,\\
\label{beta_1}
&\beta_1 = 3E + 2l(1 - 2\nu)^3 + 4m(1 + \nu)^2 (1 - 2\nu) + 6n\nu^2,\\
\label{beta_2}
&\beta_2 = 2 (1 + \nu) \left[2 l (1 - 2 \nu)^3 + \nu \left(E + 4m \left(1 - \nu - 2\nu^2\right) - 2n (1 - 2\nu)\right) \right],\\
\label{beta_3}
&\beta_3 = 2(1 + \nu)(1 - 2 \nu) \Big[ (1 + \nu)(1 - 2\nu) [4l \left(1 - 2\nu\right) - 2m(1 + 2\nu) + n] - 2\nu E \Big]\\
\label{F1}
&F^{(1)} = \frac{1}{4} \left[(1 + \nu + 2\nu^2) P_{xxx} - (1 - \nu + 2\nu^2 + 4\nu^3) P_{xtt}\right],\\
\label{F2}
&F^{(2)} = \frac{1}{4} \bigg[ (1 + \nu) P_{xxx} - \frac{1 + \nu - 2\nu^2 - 2\nu^3}{1 - \nu} P_{xtt} - \frac{3 - 5\nu - 4\nu^2 + 4\nu^3}{2(1 - \nu)} T_{tt} - 2\nu T_{xx}\bigg].
\end{align}
Дифференцируя \eqref{eq_u0_fin} по $x$, получаем два уравнения для продольной <<деформации>> $u = U_{0x}$:
\begin{equation}\label{eq_u0_def_fin}
\begin{split}
u_{tt} - u_{xx} &- 2\left(\nu P_{xx} + T_x\right) - \frac{\varepsilon}{2 E} \left(\beta_1 u^2 + 2 \beta_2 u P + \beta_3 P^2\right)_{xx}\\
& + \delta^2 \left(\alpha_1^{(i)} u_{tttt} + \alpha_2^{(i)} u_{xxtt} + \alpha_3^{(i)} u_{xxxx} + F^{(i)}_x \right) + O(\varepsilon^2, \varepsilon\delta^2, \delta^4) = 0, \quad i = 1,2.
\end{split}
\end{equation}

Три различных асимптотических модели следует из уравнений \eqref{eq_u0_def_fin} в зависимости от соотношения между двумя малыми параметрами $\varepsilon$ и $\delta$. Во-первых, если нелинейность сильно слабее дисперсии, т.е. $ \varepsilon\ll\delta^2\ll1 $, мы можем асимптотически свести \eqref{eq_u0_def_fin} к линейным уравнениям:
\begin{equation}\label{eq_u0_def_3}
\begin{split}
u_{tt} - u_{xx} - 2\left(\nu P_{xx} + T_x\right) + \delta^2 \left(\alpha_1^{(i)} u_{tttt} + \alpha_2^{(i)} u_{xxtt} + \alpha_3^{(i)} u_{xxxx} + F^{(i)}_x \right) + O(\delta^4) = 0, \\
\quad i = 1,2,
\end{split}
\end{equation}
из которых следует, что эволюция волн будет происходить главным образом под влиянием дисперсии.
Во-вторых, если нелинейность намного сильнее дисперсии, т.е. $ \delta^2\ll\varepsilon\ll1 $, мы получаем уравнение
\begin{equation}\label{eq_u0_def_2}
u_{tt} - u_{xx} - 2\left(\nu P_{xx} + T_x\right) - \frac{\varepsilon}{2 E} \left(\beta_1 u^2 + 2 \beta_2 u P + \beta_3 P^2\right)_{xx} + O(\varepsilon^2) = 0,
\end{equation}
означающее, что эволюция волн определяется нелинейностью.
Наконец, если нелинейное и дисперсионные слагаемые уравновешивают друг друга, т.е. $ \varepsilon \sim \delta^2 $, мы получаем <<модель максимального баланса>> (''maximal balance model'' согласно терминологии, используемой в~\cite{Ablowitz2011}):
\begin{equation}\label{eq_u0_def_1}
\begin{split}
u_{tt} - u_{xx} - 2&\left(\nu P_{xx} + T_x\right) - \varepsilon \bigg[\frac{1}{2E} \left(\beta_1 u^2 + 2 \beta_2 u P + \beta_3 P^2\right)_{xx}\\
& - \frac{\delta^2}{\varepsilon} \left(\alpha_1^{(i)} u_{tttt} + \alpha_2^{(i)} u_{xxtt} + \alpha_3^{(i)} u_{xxxx} + F^{(i)}_x \right)\bigg] + O(\varepsilon^2) = 0, \quad i = 1,2.
\end{split}
\end{equation}
Последняя асимптотическая модель \eqref{eq_u0_def_1} (обе её версии $i = 1,2$) является уравнением типа Буссинеска. Хорошо известно, что такие уравнения имеют решения в виде солитонов сжатия~\cite{S_book}.

Отбросим члены порядка $O(\varepsilon^2)$ в уравнениях \eqref{eq_u0_def_1} и запишем их в размерном виде, не меняя обозначения для размерных переменных:
\begin{equation}\label{eq_dim}
\begin{split}
u_{tt} - c^2 u_{xx} - \frac{2}{\rho}\bigg(\nu&P_{xx} + \frac1R T_x\bigg) - \left(\frac{\beta_1}{2\rho} u^2 + \frac{\beta_2}{\rho E} u P + \frac{\beta_3}{2\rho E^2} P^2\right)_{xx}\\
& + R^2 \bigg(\frac{\alpha_1^{(i)}}{c^2} u_{tttt} + \alpha_2^{(i)} u_{xxtt} + c^2\alpha_3^{(i)} u_{xxxx} + G^{(i)} \bigg) = 0, \quad i = 1,2,
\end{split}
\end{equation}
где $c^2 = E/\rho$, а размерные функции $G^{(i)}$ представляются в виде:
\begin{align}
G^{(1)} &= \frac{1 + \nu + 2\nu^2}{4\rho} P_{xxxx} - \frac{1 - \nu + 2\nu^2 + 4\nu^3}{4E} P_{xxtt}, \label{G1}\\
G^{(2)} &= \frac{1 + \nu}{4\rho} P_{xxxx} - \frac{1 + \nu - 2\nu^2 - 2\nu^3}{4E(1 - \nu)} P_{xxtt} - \frac{3 - 5\nu - 4\nu^2 + 4\nu^3}{8ER(1 - \nu)} T_{xtt} - \frac{\nu}{2\rho R} T_{xxx}. \label{G2}
\end{align}

Уравнения (\ref{eq_u0_def_1}) были выведены для случая сильных поверхностных напряжений, когда соответствующие слагаемые находятся в ведущем порядке по малому параметру $\varepsilon$. Если напряжения сравнительно невелики:
$P = \varepsilon \hat P$, $T = \varepsilon \hat T$,
тогда в \eqref{eq_u0_def_1} напряжения <<сдвигаются>> в следующий порядок по $\varepsilon$, что в размерном виде приводит следующим уравнениям:
\begin{equation}\label{eq_dim_weak_trac}
\begin{split}
u_{tt} - c^2 u_{xx} - \frac{2}{\rho}\left(\nu P_{xx} + \frac1R T_x\right)  &- \left(\frac{\beta_1}{2\rho} u^2 \right)_{xx} \\
+ R^2& \bigg(\frac{\alpha_1^{(i)}}{c^2} u_{tttt} + \alpha_2^{(i)} u_{xxtt} + c^2\alpha_3^{(i)} u_{xxxx}  \bigg) = 0, \quad i = 1,2.
\end{split}
\end{equation}
Отметим, что в случае условия свободной поверхности, т.е. $P = T = 0$, уравнения \eqref{eq_dim} и~\eqref{eq_dim_weak_trac} сводятся к следующему:
\begin{equation}\label{eq_dim_free_surf}
u_{tt} - c^2 u_{xx} = \frac{\beta_1}{2\rho}\left(u^2\right)_{xx} -  R^2 \bigg(\frac{\alpha_1^{(i)}}{c^2} u_{tttt} + \alpha_2^{(i)} u_{xxtt} + c^2\alpha_3^{(i)} u_{xxxx}\bigg), \quad i = 1,2.
\end{equation}

Сравним оба уравнения \eqref{eq_dim_free_surf} с <<уравнением с двумя дисперсиями>>, полученным Самсоновым и Порубовым~\cite{SP}:
\begin{equation}\label{eq_dim_SP}
u_{tt} - c^2 u_{xx} =  \frac{\beta_1}{2 \rho} (u^2)_{xx} - \frac{\nu (1-\nu) R^2}{2} u_{xxtt} + \frac{\nu c^2 R^2}{2} u_{xxxx},
\end{equation}
и <<регуляризованным>> уравнением, выведенным Островским и Сутиным~\cite{OS}:
\begin{equation}\label{eq_dim_OS}
u_{tt} - c^2 u_{xx} =  \frac{\beta_1}{2 \rho} (u^2)_{xx} + \frac{\nu^2 R^2}{2} u_{xxtt}.
\end{equation}
Все четыре модели имеют одинаковое нелинейное слагаемое, однако дисперсионные слагаемые отличаются. Уравнения \eqref{eq_dim_SP} и \eqref{eq_dim_OS} могут быть записаны в форме уравнений \eqref{eq_dim_free_surf} с помощью следующих дисперсионных коэффициентов:
\begin{align} \nonumber
&\alpha_1^{(3)} = 0,& &\alpha_2^{(3)} = \frac{(1-\nu)\nu}{2},&  &\alpha_3^{(3)} = -\frac \nu 2,\\
\nonumber
&\alpha_1^{(4)} = 0,& &\alpha_2^{(4)} = -\frac{\nu^2}{2},&  &\alpha_3^{(4)} = 0.
\end{align}

Отметим, что все четыре приведённые выше модели не являются асимптотически точными уравнениями, т.е. в безразмерной форме они содержат как члены $ O(1) $, так и $ O(\varepsilon) $. Следовательно, все эти уравнения могут быть <<регуляризованы>> (сведены) к одному уравнению, в котором есть только одно дисперсионное слагаемое, используя соотношение в главном порядке $ u_{tt} = c^2 u_{xx} + \mbox{<\,малые члены\,>} $. Коэффициент при этом дисперсионном слагаемом определяется суммой дисперсионных коэффициентов $\alpha_j$ и одинаков для всех четырех уравнений:
\begin{equation} \label{alpha_sum}
\alpha_1^{(i)} + \alpha_2^{(i)} + \alpha_3^{(i)} = -\frac{\nu^2}{2}, \quad i = \overline{1,4},
\end{equation}
что означает, что эти уравнения асимптотически эквивалентны.

Поскольку модель с одним дисперсионным членом проще, чем модель с тремя дисперсионными членами, представляется целесообразным получить регуляризованную модель с внешним воздействием. Из безразмерного уравнения \eqref{eq_u0_def_1} следует асимптотическое соотношение 
\begin{equation}\label{key}
u_{tt} = u_{xx} + 2(\nu P_{xx} + T_{x}) + O(\varepsilon),
\end{equation}
с помощью которого можно выразить $u_{tttt}$ и $u_{xxxx}$ через $u_{xxtt}$. Получаемая таким образом из (\ref{eq_u0_def_1}, $i=1$) модель в размерной форме имеет вид:
\begin{equation}\label{2_eq_fin_reg}
\begin{split}
&u_{tt} - c^2 u_{xx} - \frac{2}{\rho}\bigg(\nu P_{xx} + \frac1R T_x\bigg) - \left(\frac{\beta_1}{2\rho} u^2 + \frac{\beta_2}{\rho E} u P + \frac{\beta_3}{2\rho E^2} P^2\right)_{xx} - \frac{\nu^2 R^2}{2} u_{xxtt}\\
& + \frac{R^2}{4} \lb\frac{1-\nu}{\rho}P_{xxxx} - \frac{1-3\nu+4\nu^3}{E}P_{xxtt}\rb + \frac{(1+\nu)R}{2}\lb \frac{1}{E}T_{xtt} - \frac{1}{\rho}T_{xxx} \rb = 0, \quad i = 1,2.
\end{split}
\end{equation}
Уравнение \eqref{2_eq_fin_reg} является обобщением уравнения \eqref{eq_dim_OS} на случай ненулевых напряжений на поверхности.

Отметим, что в некоторых исследованиях, в частности, Самсонова и Порубова~\cite{SP, S_book, P_book}, для вывода модели типа Буссинеска использовались асимптотические разложения перемещений по малому параметру, а не степенные разложения по поперечной координате (радиусу). В следующем пункте мы кратко покажем, что, используя такой подход, можно получить уравнение (\ref{eq_dim},~$i=1$).



\subsection{Вывод с помощью асимптотического разложения}
В этом пункте предложен вывод уравнения \eqref{eq_u0_def_fin} из полной постановки задачи \eqref{2_eq1_0} и \eqref{2_eq2_0}, используя менее жёсткие предположения о форме асимптотических разложений перемещений и, следовательно, обосновывая разложения \eqref{u_series} в виде степенного ряда по радиальной переменной, использованные в предыдущем разделе.

Введём безразмерные переменные так, как это было сделано ранее в \eqref{scales1}, изменив лишь масштаб радиуса:
\begin{equation}\label{2_scales_modif}
\tilde r = \frac{r}{\delta L}
\end{equation}
чтобы безразмерный радиус $\tilde r$ был порядка $1$, а не $\delta$. Как и ранее, опустим тильду над безразмерными величинами в дальнейшем изложении. Будем искать безразмерные перемещения в виде асимптотических разложений по малому параметру $\delta$:
\begin{align}
\label{expand_delta1}
U(x, r, t) &= U_0(x, t) + U_2(x, r, t) \delta^2 + U_4(x, r, t) \delta^4 + O(\delta^6),\\
\label{expand_delta2}
V(x, r, t) &= V_1(x, r, t)\delta + V_3(x, r, t) \delta^3 + V_5(x, r, t) \delta^5 + O(\delta^7).
\end{align}
Здесь для краткости уже использовано предположение о том, что $U_0$ не зависит от $r$, которое следует из соотношения в главном порядке по $\delta$ в уравнении \eqref{2_eq1_0}. Кроме того, здесь мы сразу опускаем нечётные степени по $\delta$ в $U_0$ и чётные в $V_0$, что может быть доказано в ходе подробного, но громоздкого вывода, который здесь приводить не будем.

Подставляя разложения \eqref{expand_delta1} и \eqref{expand_delta2} в уравнения движения \eqref{2_eq1_0} \eqref{2_eq2_0}, получаем
\begin{eqnarray}\label{eq1_2}
\begin{split}
&\rho c^2 U_{0tt} - (\lambda + 2\mu) U_{0xx} - (\lambda + \mu) \left(V_{1xr} + \frac{V_{1x}}{r}\right) - \mu\left(U_{2rr} + \frac{U_{2r}}{r}\right) + \varepsilon\widetilde{\Phi}_1(U_0, V_1, U_2) \\
&+ \delta^2 \left[\rho c^2 U_{2tt} - (\lambda + 2\mu) U_{2xx} - (\lambda + \mu) \left(V_{3xr} + \frac{V_{3x}}{r}\right) - \mu\left(U_{4rr} + \frac{U_{4r}}{r}\right)\right] + O(\varepsilon^2, \varepsilon\delta^2, \delta^4) = 0,
\end{split} 
\end{eqnarray}
\begin{equation}\label{eq2_2}
\begin{split}
(\lambda + 2\mu) \left(\frac{V_1}{r^2} - \frac{V_{1r}}{r} - V_{1rr}\right) &+ \varepsilon \widetilde{\Phi}_2(U_0, V_1) + \delta^2 \bigg[\rho c^2 V_{1tt} - \mu V_{1xx} - (\lambda + \mu)U_{2xr} \\
& + (\lambda + 2\mu)\left(\frac{V_3}{r^2} - \frac{V_{3r}}{r} - V_{3rr}\right) \bigg] + O(\varepsilon^2, \varepsilon\delta^2, \delta^4) = 0,
\end{split}
\end{equation}
где функции $\widetilde{\Phi}_1$ и $\widetilde{\Phi}_2$ включают все нелинейные члены (для краткости не приводятся здесь). Нижний индекс $x$, $r$ и $t$, как и раньше, обозначают частную производную по соответствующей переменной.

Граничные условия при $r = 1$ принимают вид:
\begin{equation} \label{2_bc_rr_2}
\begin{split}
\lambda U_{0x} + (\lambda + 2\mu)V_{1r} + \lambda V_1 + \varepsilon\widetilde{\Psi}_1(U_0, V_1, U_2, V_3) &+ \delta^2\big[\lambda  U_{2x} + (\lambda + 2\mu)V_{3r} + \lambda V_3\big]\\
&+ O(\varepsilon^2, \varepsilon\delta^2, \delta^4) =  \frac{\mu(3\lambda + 2\mu)}{\lambda + \mu} P(x,t),
\end{split}
\end{equation}
\begin{equation} \label{2_bc_rx_2}
\mu \left(U_{2r} + V_{1x}\right) + \varepsilon\widetilde{\Psi}_2(U_0, V_1, U_2, V_3) + \delta^2 \mu \left(U_{4r}+V_{3x}\right) + O(\varepsilon^2, \varepsilon\delta^2, \delta^4) = \frac{\mu(3\lambda + 2\mu)}{\lambda + \mu} T(x,t).
\end{equation}
Собирая коэффициенты при одинаковых степенях $\delta$ в уравнениях \eqref{eq1_2} и \eqref{eq2_2}, получаем систему нелинейных обыкновенных дифференциальных уравнений по переменной $r$, где все нелинейные слагаемые порядка $n$ умножены на $\varepsilon^{n-1}$. Мы решаем эти уравнения с граничными условиями, следующими из \eqref{2_bc_rr_2} и \eqref{2_bc_rx_2}, используя асимптотические разложения функций по $\varepsilon$. Покажем эту процедуру на примере функции $V_1$, которую представим в следующем виде:
\begin{equation}\label{v1_eps_expantion}
V_1(x,r,t) = f(x,r,t) + \varepsilon g(x,r,t) + \O{\varepsilon^2},
\end{equation}
где $f$ и $g$ неизвестные функции. Подставляя разложение \eqref{v1_eps_expantion} в уравнение \eqref{eq2_2}, получаем ОДУ относительно $f$ в главном порядке по $\varepsilon$:
\begin{equation}\label{v1_f_ode}
f_{rr} + \frac{f_r}{r} - \frac{f}{r^2} = 0,
\end{equation}
общее решение которого имеет вид:
\begin{equation}\label{v1_f_gen_sln}
f(x,r,t) = C_1(x,t)r + \frac{C_2(x,t)}{r} .
\end{equation}
Неизвестные функции $C_1$ и $C_2$ находятся из граничного условия, следующего из \eqref{2_bc_rr_2}, а также из условия симметрии, согласно которому поперечное перемещение равно нулю в центре стержня:
\begin{align}
\label{v1_f_bc1}
(\lambda + 2\mu)f_r + \lambda f = -\lambda U_{0x} + \frac{\mu(3\lambda + 2\mu)}{\lambda + \mu} P(x,t) \quad\text{при}\quad r = 1,\\
\label{v1_f_bc2}
f = 0 \quad\text{при}\quad r = 0.
\end{align}
Из \eqref{v1_f_bc2} следует, что $C_2 \equiv 0$, а $C_1$ находится из \eqref{v1_f_bc1}, что даёт итоговое выражение для~$f$:
\begin{equation}\label{v1_f_sln}
f(x,r,t) = \frac{r}{2(\lambda + \mu)} \left(\frac{ \mu(3\lambda + 2\mu)}{\lambda + \mu} P - \lambda U_{0x}\right).
\end{equation}
Используя \eqref{v1_f_sln}, получаем уравнение относительно функции~$g$ в следующем порядке по~$\varepsilon$ с граничными условиями в виде:
\begin{align}
\label{v1_g_ode}
\quad g_{rr} + \frac{g_r}{r} - \frac{g}{r^2} &= 0,\\
\label{v1_g_bc1}
(\lambda + 2\mu)g_r + \lambda g &= a_1 U_{0x}^2 + a_2 U_{0x} P + a_3 P^2 \hspace{-2cm}&\text{при}\quad r = 1,\\
\label{v1_g_bc2}
g &= 0 \hspace{-2cm}&\text{при}\quad r = 0.
\end{align}
Решение задачи \eqref{v1_g_ode} -- \eqref{v1_g_bc2} имеет форму:
\begin{equation}\label{v1_g_sln}
g(x,t) = \frac{r (a_1 U_{0x}^2 + a_2 U_{0x} P + a_3 P^2)}{2(\lambda + \mu)}.
\end{equation}

Функции $U_2$, $V_3$ и $U_4$ исключаются аналогичным образом, при использовании ещё одного условия $U_r = 0$ при $r = 0$. В конечном итоге получаем:
\begin{align}
\begin{split}
V_1(x, r, t) &= \frac{r}{2(\lambda + \mu)} \left(\frac{ \mu(3\lambda + 2\mu)}{\lambda + \mu} P - \lambda U_{0x}  + \varepsilon (a_1 U_{0x}^2 + a_2 U_{0x} P + a_3 P^2) \right) + O(\varepsilon^2),
\end{split}\\
\begin{split}
U_2(x, r, t) &= \frac{r^2}{4\mu} \left(\rho c^2 U_{0tt} - 2\mu U_{0xx} - \frac{\mu(3\lambda + 2\mu)}{\lambda + \mu} P_x\right) \\
&\quad + \varepsilon r^2 \big[U_{0x} \left(a_4 U_{0tt} + a_5 U_{0xx} + a_6 P_x\right) + P(a_7 U_{0tt} + a_8 U_{0xx} + a_9 P_x)\big] + O(\varepsilon^2),
\end{split}\\
\begin{split}
V_3(x, r, t) &= r \left(b_1(r) P_{tt} + b_2(r) U_{0xtt} + b_3(r) P_{xx} + b_4(r) U_{0xxx}\right) + O(\varepsilon),
\end{split}\\
\begin{split}
U_4(x, r, t) &= r^4 a_{10} U_{0tttt} + r^2 \left(b_5(r) U_{0xxxx} + b_6(r) U_{0xxtt} + b_7(r) P_{xtt} + b_8(r) P_{xxx}\right) + O(\varepsilon).
\end{split}
\end{align}
Здесь функции $b_i(r) = b_{i}^{(2)}r^2 + b_{i}^{(0)}$, коэффициенты $a_i$ и $b_{i}^{(j)}$ зависят от упругих модулей $\lambda, \mu, l, m, n$ и плотности $\rho$.
Итоговое уравнение, следующее из \eqref{2_bc_rx_2} после подстановки $V_1$, $U_2$, $V_3$ и $U_4$, совпадает с ранее выведенным уравнением~\eqref{eq_u0_asympt}, из которого следует уравнение~(\ref{eq_dim},~$i=1$).


\section{Вывод уравнения типа Буссинеска с внешним воздействием в растянутом стержне}
Некоторые исследователи рассматривали задачу о распространении длинных продольных волн в предварительно растянутом стержне~\cite{DF}. Вывод уравнений, предложенный в предыдущем параграфе, позволяет легко учесть это растяжение.

Рассмотрим распространение продольных волн в равномерно растянутом вдоль своей оси стержне. Продольное перемещение в растянутом состоянии представляется в виде:
\begin{equation}\label{pre_stretch_u}
U^*(x) = \kappa x,
\end{equation}
где $\kappa$ --- постоянная.
Обезразмерим предварительное растяжение с помощью того же масштаба, что $U$ в \eqref{scales1}:
\begin{equation} \label{scales_pre_stretch}
\tilde U^* = \frac{U^*}{\varepsilon L} = \tilde\kappa \tilde x, \quad \mbox{где} \quad  \tilde \kappa = \frac{\kappa}{\varepsilon}.
\end{equation}
Более того, будем предполагать, что в предварительно растянутом состоянии на стержень не действуют внешние напряжения, приложенные к его поверхности. Решая уравнения движения \eqref{2_eq1_0} и \eqref{2_eq2_0} с граничными условиями \eqref{2_bc_rr} и \eqref{2_bc_rx} при $P = T = 0$, записанными в безразмерном виде с помощью \eqref{scales1} и \eqref{scales_pre_stretch}, получаем поперечное перемещение $\tilde V^*$ в растянутом стержне:
\begin{equation}\label{pre_stretch_v}
%\begin{split}
\tilde V^*(\tilde r) = -\,\frac{\lambda \tilde\kappa \tilde r}{2(\lambda+\mu)} \bigg(1 + \varepsilon \frac{\tilde\kappa \left(2\mu^2 (\lambda + 2l) + \lambda^2 (3\lambda + 6m - 2n) + \lambda\mu(5\lambda + 4m - 2n)\right)}{4\lambda(\lambda + \mu)^2} 
+ O(\varepsilon^2)\bigg).
%\end{split}
\end{equation}
Введём новые безразмерные разложения перемещений (тильды опущены):
\begin{align}
\label{pre_stretch_u_series_scaled}
U(x,r,t) &= U^*(x) + U_0 + r^2 U_2 + r^4 U_4 + O(r^6),\\
\label{pre_stretch_v_series_scaled}
V(x,r,t) &= V^*(r) + r V_1 + r^3 V_3 + r^5 V_5 + O(r^7).
\end{align}
Следуя выводу уравнений \eqref{eq_u0_def_fin} и используя разложения \eqref{pre_stretch_u_series_scaled}, \eqref{pre_stretch_v_series_scaled} взамен \eqref{u_series_scaled}, \eqref{v_series_scaled}, получаем уравнения:
\begin{equation}\label{pre_stretch_eq_u0_def_fin}
\begin{split}
u_{tt} - \left(1 + \varepsilon\kappa\frac{\beta_1}{E}\right)u_{xx} - 2 \left[\left(\nu + \varepsilon\kappa \frac{\beta_2}{2E} \right) P_{xx} + T_x\right] - \varepsilon \left(\frac{\beta_1}{2 E} u^2 + \frac{\beta_2}{E} u P + \frac{\beta_3}{2E} P^2\right)_{xx}\\
+ \delta^2 \left(\alpha_1^{(i)} u_{tttt} + \alpha_2^{(i)} u_{xxtt} + \alpha_3^{(i)} u_{xxxx} + F^{(i)}_x \right) + O(\varepsilon^2, \varepsilon\delta^2, \delta^4) = 0, \quad i = 1,2,
\end{split}
\end{equation}
где использованы обозначения, введённые в предыдущих разделах. Отметим, что здесь $u = U_{0x}$ является возмущением относительно растянутого состояния, в то время как в уравнениях \eqref{eq_u0_def_fin} оно обозначает возмущение относительно недеформированного состояния.

Предполагая, что нелинейные и дисперсионные слагаемые одного порядка ($\varepsilon \sim \delta^2$) и отбрасывая малые члены в \eqref{pre_stretch_eq_u0_def_fin}, получаем уравнение, которое в размерном виде представляется следующим образом:
\begin{equation}\label{pre_stretch_eq_dim}
\begin{split}
u_{tt} - \left(c^2 + \kappa\frac{\beta_1}{\rho}\right) u_{xx} - \frac{2}{\rho}\left[\left(\nu + \kappa \frac{\beta_2}{2E} \right) P_{xx} + \frac1R T_x\right] - \left(\frac{\beta_1}{2\rho} u^2 + \frac{\beta_2}{\rho E} u P + \frac{\beta_3}{2\rho E^2} P^2\right)_{xx}\\
+ R^2 \left(\frac{\alpha_1^{(i)}}{c^2} u_{tttt} + \alpha_2^{(i)} u_{xxtt} + c^2\alpha_3^{(i)} u_{xxxx} + G^{(i)} \right) = 0, \quad i = 1,2.
\end{split}
\end{equation}
Здесь коэффициенты $\alpha_j^{(i)}$, $\beta_j$ и функции $G^{(i)}$ Задаются формулами \eqref{alpha_1} -- \eqref{beta_3} и (\ref{G1}) -- (\ref{G2}) соответственно.
Предварительное растяжение изменило скорость длинных линейных волн, квадрат которой равен коэффициенту при $u_{xx}$. 
Это явление называется акустоэластическим эффектом и изучалось в~\cite{HughesKelly, ADO}. Отметим, что акустоэластический эффект используется для экспериментального определения упругих модулей Мурнагана.

Насколько известно автору, обе модели, описываемые уравнениями (\ref{pre_stretch_eq_dim}), а также их упрощённые версии (\ref{eq_dim}), (\ref{eq_dim_free_surf}) и \eqref{2_eq_fin_reg}, получены впервые. В следующем параграфе мы проанализируем свойства полученных уравнений и сравним их с уравнениями выведенными ранее.




\section{Дисперсионные свойства и солитонные решения}

На рисунке \ref{fig:disp} представлены дисперсионные кривые четырёх упрощённых (с нулевыми напряжениями на поверхности и без предварительного растяжения) линеаризованных уравнения типа Буссинеска, приведённые в предыдущих разделах, а также нижние три ветви точного дисперсионного соотношения Похгаммера-Кри для линейной задачи.
Дисперсионные соотношения этих моделей имеют вид:
\begin{align}
&\frac{2p}{R}\left(q^2 + k^2\right) J_1(pR) J_1(qR) - \left(q^2 - k^2\right)^2 J_0(pR) J_1(qR) -4k^2 p q J_1(pR) J_0(qR) = 0,\\
& \alpha_1^{(i)} \overline{\omega}^4 - \left(1 - \alpha_2^{(i)} \overline{k}^2\right) \overline{\omega}^2 + \overline{k}^2 \left(1 + \alpha_3^{(i)} \overline{k}^2 \right) = 0, \quad i = 1,2,\\
& \left(1 - \frac{(1-\nu)\nu}{2} \overline{k}^2\right) \overline{\omega}^2 - \overline{k}^2 \left(1 - \frac{\nu  \overline{k}^2}{2}\right) = 0,\\
& \left(1 + \frac{\nu^2}{2} \overline{k}^2\right) \overline{\omega}^2 - \overline{k}^2 = 0,
\end{align}
для решения Похгаммера-Кри и уравнений (\ref{eq_dim_free_surf}, $i=1,2$), \eqref{eq_dim_SP} и \eqref{eq_dim_OS} соответственно. Здесь $\overline{k} = k R$, $\overline{\omega} = \omega R / c$, где $k$ и $\omega$ --- волновое число и волновая частота соответственно, $J_i$ --- функция Бесселя первого рода, а параметры $p$ и $q$ выражаются следующим образом:
\begin{equation*}
p^2 = \frac{\rho \overline{\omega}^2}{\lambda + 2\mu} - \overline{k}^2, \quad q^2 = \frac{\rho \overline{\omega}^2}{\mu} - \overline{k}^2.
\end{equation*}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.85\linewidth]{2_DispBlack}
	\caption{Дисперсионные кривые для стержня с $\nu = 0.34$.}
	\label{fig:disp}
\end{figure} 

Все модели достаточно хорошо описывают нижнюю ветвь дисперсионной кривой в длинноволновой области, однако наиболее точной является модель (\ref{eq_dim_free_surf}, $i=2$). Уравнение Самсонова -- Порубова \eqref{eq_dim_SP} обладает коротковолновой неустойчивостью, в то время как остальные три модели не имеют такого эффекта. Отметим, что коротковолновая неустойчивость затрудняет численный счёт, поскольку высокочастотные гармоники в таком случае могут неограниченно возрастать. Полученные в настоящей работе уравнения \eqref{eq_dim_free_surf} в отличие от других уравнений улавливают вторую ветвь дисперсионной кривой, правда описывают её очень неточно: помимо большого отличия по значению, эти кривые имеют всюду положительный наклон, тогда как точная кривая имеет отрицательный наклон в области длинных волн, что соответствует отрицательной групповой скорости.

Все четыре уравнения (\ref{eq_dim_free_surf}, $i=1,2$), \eqref{eq_dim_SP} и \eqref{eq_dim_OS} имеют семейство солитонных решений, вывод которых описан в Приложении 2:
\begin{equation}\label{soliton}
u_i(x,t) = A\ {\rm sech}^2\ \left[B_{i} \left(x\pm t \sqrt{c^2+\frac{A \beta_1}{3 \rho}}\right) \right], \quad i = \overline {1,4}.
\end{equation}
Здесь амплитуда $A$ является свободным параметром, причём при $A<0$ такая волна называется солитоном сжатия, а при $A>0$ --- солитоном разрежения. Для заданной амплитуды~$A$, соответствующие солитонные решения имеют одинаковую скорость, но разные параметры длины~$B_i$:
\begin{align}
\label{Bi}
B_i &= \sqrt{\frac{3A\beta_1 E}{-4\left[(A\beta_1 + 3E)^2\alpha_1^{(i)} + 3E(A\beta_1 + 3E)\alpha_2^{(i)} + 9E^2\alpha_3^{(i)}\right] R^2}} \, , \quad i=1,2,\\
\label{B3}
B_3 &= \sqrt{\frac{A\beta_1}{\left[6\nu E + 2 A \beta_1 (\nu - 1)\right] \nu R^2}}\, , \\
\label{B4}
B_4 &= \sqrt{\frac{A\beta_1}{(6E + 2A\beta_1)\nu^2 R^2}},
\end{align}
для уравнений \eqref{eq_dim_free_surf}, \eqref{eq_dim_SP} и \eqref{eq_dim_OS} соответственно.

На рисунке \ref{fig:soliton} в левой части изображены четыре солитона сжатия, задаваемых формулами \eqref{soliton} -- \eqref{B4} и имеющих амплитудный параметр $A = -0.05$.  <<Регуляризованный>> солитон \eqref{B4} и солитон (\ref{Bi}, $i=1$) практически полностью совпадают и являются самыми длинными из всех.
Однако солитоны такой амплитуды вызывают напряжения близкие к пределу упругости для полистирола. В экспериментах с полистироловым стержнем, описанных в~\cite{Garbuzov}, амплитуда очень мала: $A \sim 10^{-3} - 10^{-4}$, и, следовательно, в во всех четырёх формулах параметр длины примерно равен
\begin{equation}\label{B}
B = \sqrt{\frac{A\beta_1}{6\nu^2 E R^2}},
\end{equation}
а соответствующее солитонное решение изображено в правой части рисунка~\ref{fig:soliton} для $A = -0.001$.
Отметим, что для полистирола, упругие характеристики которого представлены в таблице \ref{tab:ps}, коэффициент $\beta_1$, задаваемый формулой~\eqref{beta_1} отрицателен и, следовательно, параметр $B$ вещественен. Более того, из формулы \eqref{B} следует, что при малых значениях $A$ тип солитона (сжатия или растяжения) определяется именно знаком коэффициента~$\beta_1$.
\begin{figure}[h]
	\centering
	\includegraphics[width=0.54\linewidth]{3a_FourSolitonsBlack}
	\includegraphics[width=0.42\linewidth]{3b_SingleSoliton}
	\caption{Графики функций $-u_i(x,t)$, задаваемых формулой \eqref{soliton}, в стержне радиуса $R = 5$~мм, сделанном из полистирола при $A = -0.05$ (слева) и $A= -0.001$ (справа). Упругие модули полистирола приведены в таблице~\ref{tab:ps}.}
	\label{fig:soliton}
	\vspace{-3mm}
\end{figure}
\begin{table}[h]
	\captionsetup{justification=raggedleft,singlelinecheck=false}
	\caption{Упругие модули полистирола~\cite{HughesKelly}.}
	\vspace{-7mm}
	\begin{center}
		\begin{tabular}{|c|c|c|c|c|c|}
			\hline
			\rule[-1ex]{0pt}{3ex} Модуль Юнга & Коэффициент & \multicolumn{3}{|c|} {Модули Мурнагана, Н/м\textsuperscript{2} } & Плотность \\
			\cline{3-5}
			\rule[-1ex]{0pt}{3ex} $E$, Н/м\textsuperscript{2} & Пуассона, $\nu$ & $l$ & $m$ & $n$ & $\rho$, кг/м\textsuperscript{3}  \\
			\hline
			\rule[-1ex]{0pt}{3ex} $3.7\cdot10^9$ & $0.34$ & $-18.9\cdot10^{9}$ & $-13.3\cdot10^{9}$ & $-10\cdot10^{9}$ & 1060 \\
			\hline
		\end{tabular}
	\end{center}
	\label{tab:ps}
\end{table}

На рисунке~\ref{fig:rod_deformed} показано, как выглядит деформированный стержень, когда по нему бежит солитон сжатия. Увеличенные перемещения позволяют увидеть эффект Пуассона (утолщение тела при продольном сжатии), а также нелинейную зависимость продольного перемещения от~$r$ (вертикальные линии немного изгибаются в месте сжатия). На всякий случай напомним, что на рисунке~\ref{fig:soliton} изображён график продольной деформации на оси стержня, а не форма границы стержня на рисунке~\ref{fig:rod_deformed}. 
\begin{figure}[hh]
	\centering
	\includegraphics[width=0.6\linewidth]{4a_DeformedRod1}
	\includegraphics[width=0.6\linewidth]{4b_DeformedRod2}
	\caption{Стержень из полистирола при прохождении солитона деформации~(\ref{soliton},~${i=1}$) с амплитудным параметром $A=-0.03$. Продольные перемещения увеличены в 30 раз, поперечные в 60 раз.}
	\label{fig:rod_deformed}
\end{figure}

Репараметризуем солитонное решение \eqref{soliton} через скорость $v$ вместо амплитуды $A$:
\begin{align}
\label{soliton_v}
u_i(x, t) &= \frac{3 \rho \left(v^2 - c^2\right)}{\beta_1} {\rm sech}^{2}\left[\tilde{B_i}(x\pm v t)\right], \qquad v = \sqrt{c^2+\frac{A\beta_1}{3 \rho}} \, ,\\
\label{Bimod}
\tilde{B}_i &= \sqrt{\frac{c^2(v^2- c^2)}{-4\left(\alpha_1^{(i)} v^4 + \alpha_2^{(i)} c^2 v^2 + \alpha_3^{(i)}c^4\right)R^2}} \, , \quad i =1,2,\\
\label{B3mod}
\tilde{B}_3 &= \sqrt{\frac{v^2- c^2}{2\nu R^2 [c^2 - (1 - \nu) v^2 ]}} \, ,\\
\label{B4mod}
\tilde{B}_4 &= \sqrt{\frac{v^2- c^2}{2\nu^2 v^2 R^2}} \, .
\end{align}
Солитонное решение существует, только если параметр длины $\tilde{B}$ вещественен, что приводит, предполагая $\nu < 0.5$, к следующим ограничениям на скорость солитона:
\begin{itemize}
	\item $\tilde{B}_i^2 > 0 \implies$ $\displaystyle v^2 < \frac{-\alpha_2^{(i)} - \sqrt{\alpha_2^{(i)2} - 4\alpha_1^{(i)}\alpha_3^{(i)}}}{2\alpha_1^{(i)}} c^2$ или $\displaystyle c^2 < v^2 < \frac{-\alpha_2^{(i)} + \sqrt{\alpha_2^{(i)2} - 4\alpha_1^{(i)}\alpha_3^{(i)}}}{2\alpha_1^{(i)}} c^2$, $i=1,2$,
	\item $\tilde{B}_3^2 > 0 \implies$ $\displaystyle c^2 < v^2 < \frac{c^2}{1-\nu}$, 
	\item $\tilde{B}_4^2 > 0 \implies$ $c^2 < v^2$,
\end{itemize}
причём $\displaystyle 0 < \frac{-\alpha_2^{(i)} - \sqrt{\alpha_2^{(i)2} - 4\alpha_1^{(i)}\alpha_3^{(i)}}}{2\alpha_1^{(i)}} \leqslant 1$ и $\displaystyle \frac{-\alpha_2^{(i)} + \sqrt{\alpha_2^{(i)2} - 4\alpha_1^{(i)}\alpha_3^{(i)}}}{2\alpha_1^{(i)}} \geqslant 1$ при $\nu\in[0, 0.5]$, $i=1,2$. 
Таким образом, в то время как первые три модели дают ограниченный диапазон скоростей, регуляризованная модель не накладывает ограничения сверху на скорость солитона. Также первые две модели, в отличие от двух других, допускают существование солитонов противоположного знака. 

Помимо математического ограничения на параметры солитона (вещественность параметра длины $B$) есть физическое ограничение: солитон не должен вызывать пластических деформаций, поскольку модель строилась только для упругих деформаций. Так, для полистирола и широкого круга других материалов именно физическое ограничение сильнее всего сужает диапазон допустимых скоростей, поскольку, как видно из формулы \eqref{soliton_v}, чем больше разность ${v^2-c^2}$, тем больше амплитуда и меньше длина солитона.


\chapter{Численное решение уравнений нелинейной теории упругости}
В главе 2 из полных нелинейных уравнений движения 
%\eqref{2_eq1_0}, \eqref{2_eq2_0} и граничных условий \eqref{2_bc_rr}, \eqref{2_bc_rx} 
выведены упрощённые модели типа Буссинеска, описывающие длинные продольные волны деформации, которые теперь интересно сравнить с численным решением полных уравнений.

В настоящей работе основной интерес представляют непрерывные гладкие решения нелинейных уравнений. Наилучшим средством для их нахождения является псевдоспектральный метод, с помощью которого мы будем решать как пространственно одномерные уравнения типа Буссинеска, так и полные трёхмерные уравнения движения.
Отметим, что пседвоспектральный метод чаще всего используется только для пространственной дискретизации, в то время как дискретизация по времени осуществляется с помощью стандартных методов решения обыковенных дифференциальных уравнений, например, Рунге-Кутты 4-го порядка.
%Опишем метод подробнее.

Семейство спектральных методов, одним из которых является псевдоспектральный метод или метод коллокации, основано на поиске решения задачи в некотором подпространстве, имеющем конечный базис, в качестве которого чаще всего выбирается базис Фурье (набор синусов и косинусов) или семейство ортогональных многочленов, например, Чебышёва или Лежандра. 
Отличительной особенностью псевдоспектрального метода является выполнение уравнений не на всей области, а лишь в конечном наборе точек, называемых точками коллокации. Получаемые в результате значения решения в точках коллокации затем интерполируются на всю область задачи. 
Выбор точек коллокации определяется используемым базисом: для базиса Фурье точки равномерно распределены по области, а для ортогональных многочленов выбираются точки соответствующей интерполяционной квадратуры, задаваемые нулями определённых многочленов. Конкретный вид многочлена, у которого ищутся нули, определяется типом квадратуры: Гаусса, Гаусса-Радо или Гаусса-Лобатто. Более подробному описанию этого метода применительно к решаемым уравнениям посвящён следующий параграф.



\section{Численная схема}
\subsection{Одномерное уравнение типа Буссинеска}
Сформулируем начально-краевую задачу для одномерного регуляризованного уравнения Буссинеска с внешним воздействием:
\begin{align}
\label{3_bq_reg}
\begin{split}
&u_{tt} - c^2 u_{xx} - g_1 P_{xx} - g_2 T_x - \left(\tilde{\beta}_1 u^2 + \tilde{\beta}_2 u P + \tilde{\beta}_3 P^2\right)_{xx} + \tilde\alpha u_{xxtt}\\
&\qquad\qquad\qquad\qquad\qquad + \gamma_1 P_{xxxx} + \gamma_2 P_{xxtt} +\gamma_3 T_{xtt} + \gamma_4 T_{xxx} = 0, \quad x\in(0, L),
\end{split}\\
\label{3_bq_iv}
&u(x, 0) = \phi(x), \qquad u_t(x, 0) = \psi(x),%\\
%\label{3_bq_bc}
%&u(0, t) = u(L, t).
\end{align}
где на границах области поставлены условия симметрии, поскольку модель Буссинеска выводилась для бесконечного стержня. Уравнение \eqref{3_bq_reg} совпадает по форме с \eqref{2_eq_fin_reg}, где некоторые коэффициенты переобозначены для лаконичности записи.

Будем решать регуляризованное уравнение типа Буссинеска с помощью псевдоспектрального метода Фурье. Для этого сначала отобразим область задачи $(0, L)$ в $(0, 2\pi)$ с помощью замены $\tilde x = Sx$, где $S = 2\pi/L$, и запишем новое уравнение, опустив тильду над $x$:
\begin{equation}
\label{3_bq_reg_scaled}
\begin{split}
&u_{tt} - S^2 c^2 u_{xx} - S^2 g_1 P_{xx} - S g_2 T_x - S^2\left(\tilde{\beta}_1 u^2 + \tilde{\beta}_2 u P + \tilde{\beta}_3 P^2\right)_{xx} + S^2\tilde\alpha u_{xxtt}\\
&\qquad\qquad\qquad\qquad + S^4\gamma_1 P_{xxxx} + S^2\gamma_2 P_{xxtt} + S\gamma_3 T_{xtt} + S^3\gamma_4 T_{xxx} = 0, \quad x\in(0, 2\pi),
\end{split}
\end{equation}
Будем искать приближённое решение $u_N$ в виде комбинации из $N$ гармоник:
\begin{equation}\label{3_spect_approx}
u_N(x, t) = \sum_{k=-N/2}^{N/2-1} \widehat u(k, t) e^{ikx}, \quad x\in(0, 2\pi).
\end{equation}
Ключевая идея псевдоспектрального метода заключается в поиске таких коэффициентов $\widehat u(k, t)$, что точное решение $u$ совпадает с приближённым решением $u_N$ в точках коллокации $x_j$:
\begin{equation}\label{3_pseudo_sp_approx}
u_N(x_j, t) = u(x_j, t), \quad x_j = \frac{2\pi j}{N}, \quad j=0,\dots, N-1.
\end{equation}
Из условия \eqref{3_pseudo_sp_approx} естественным образом вытекает применение дискретного преобразования Фурье (ДПФ) для нахождения $\widehat u(k, t)$, а для обратной операции --- обратного ДПФ:
\begin{align}\label{3_dft}
\widehat u(k, t) &= \mathcal{F} u = \frac1N\sum_{j=0}^{N-1} u(x_j, t) e^{-ikx_j}, \quad k=-\frac N 2, \dots, \frac N 2 - 1,\\
\label{3_idft}
u(x_j, t) &= \mathcal{F}^{-1} \widehat u = \sum_{k=-N/2}^{N/2-1} \widehat u(k, t) e^{ikx_j}, \quad j=0, \dots, N-1.
\end{align}
Отметим, что использование ДПФ выгодно с вычислительной точки зрения, поскольку для его реализации есть алгоритм быстрого преобразования Фурье.

Спектральное представление \eqref{3_spect_approx} позволяет быстро находить производную по пространственной переменной:
\begin{equation}\label{3_spect_deriv}
u_N'(x, t) = \sum_{k=-N/2}^{N/2-1} ik \widehat u(k, t) e^{ikx}.
\end{equation}
Применим спектральное представление \eqref{3_spect_approx} к уравнению \eqref{3_bq_reg_scaled} и начальным условиям \eqref{3_bq_iv}. Приравнивая коэффициент при каждой гармонике к нулю, получаем систему обыкновенных дифференциальных уравнений с начальными условиями в виде:
\begin{gather}\label{3_bq_reg_dft}
\begin{split}
\lb 1 - \tilde{\alpha} S^2 k^2\rb\ddot{\widehat u} =&  -S^2 k^2 \lb c^2 \widehat u + \tilde \beta_1 \widehat{u^2} + \tilde\beta_2\widehat{uP} + \tilde{\beta}_3\widehat{P^2} \rb + S^2 g_1 \widehat{P_{xx}} - S g_2 \widehat{T_x} \\
&- \lb S^4\gamma_1 \widehat{P_{xxxx}} + S^2\gamma_2 \widehat{P_{xxtt}} + S\gamma_3 \widehat{T_{xtt}} + S^3\gamma_4 \widehat{T_{xxx}}\rb,
\end{split}\\
\label{3_bq_dft_iv}
\widehat{u}(k, 0) = \widehat{\phi}(k), \qquad \dot{\widehat{u}}(k, 0) = \widehat{\psi}(k), \qquad k=-\frac N 2, \dots, \frac N 2 - 1,
\end{gather}
где точка над функцией обозначает производную по времени.
Для нахождения Фурье образа нелинейного слагаемого $u^2$, а также $uP$, к функции $\widehat u$ применяется обратное преобразование, затем полученная функция $u(x_j, t)$ возводится в квадрат в узлах сетки или домножается на $P(x_j, t)$ и переводится обратно в пространство Фурье:
$$
\widehat{u^2} = \mathcal{F}\lbrack\lb\mathcal{F}^{-1}\widehat u\rb^2\rbrack, \qquad \widehat{uP} = \mathcal{F}\lbrack P\cdot \mathcal{F}^{-1}\widehat u\rbrack.
$$
Такой метод при использовании быстрого преобразования Фурье позволяет эффективно вычислять $\widehat{u^2}$ и $\widehat{uP}$.
Решать систему \eqref{3_bq_reg_dft}, \eqref{3_bq_dft_iv} будем с помощью метода Рунге-Кутты 4-го порядка. Применяя обратное ДПФ к решению этой системы, получаем приближение исходной задачи \eqref{3_bq_reg}, \eqref{3_bq_iv}.

С помощью псевдоспектрального метода уравнения в частных производных часто сводятся к системе ОДУ относительно значений функции в точках коллокации $u_N(x_j, t)$, а не коэффициентов $\widehat u(k, t)$. В таком случае нелинейные слагаемые вычисляются простым поточечным произведением значений функций в точках коллокации, а вычисление пространственных производных осуществляется через умножение вектора значений $u_N(x_j, t)$ на матрицу производной $D_N$:
\begin{equation}\label{3_phys_deriv}
\vect{u}_N'(t) = D_N \vect{u}_N(t), \qquad \vect{u}_N(t) = \{u_N(x_j, t)\}_{j=0}^{N-1}.
\end{equation}
Однако для регуляризованного уравнения Буссинеска, где присутствует смешанная производная $u_{xxtt}$, применить такой метод вычисления производной затруднительно. Матрицы $D_N$ для различных базисов приведены в~\cite{Canuto2007}.
%Для базиса Фурье она записывается следующим образом:
%\begin{equation}\label{3_diff_mat_fourier}
%(D_N)_{jl} = \begin{cases}
%\frac12(-1)^{j+l} \ctg\lbrack\frac{(j-l)\pi}{N}\rbrack, &j\neq l,\\
%0, &j=l,
%\end{cases}
%\end{equation}
%а для базиса Лежанжра в точках Гаусса-Лобатто и Гаусса-Радо так:
%\begin{equation}\label{3_diff_mat_legendre}
%(D_N)_{jl} = \begin{cases}
%\dfrac{L_N(x_j)}{L_N(x_l)} \dfrac{1}{x_j-x_l}, &j\neq l,\\
%-\dfrac{(N+1)N}{4}, &j=l=0,\\
%\dfrac{(N+1)N}{4}, &j=l=N,\\
%0, &\mbox{в остальных случаях},
%\end{cases}
%\end{equation}

Аналогичным образом можно дискретизировать уравнения (\ref{eq_dim}) с тремя дисперсионными слагаемыми, однако возникающая система обыкновенных дифференциальных уравнений 4-го порядка оказывается жёсткой и не поддаётся решению явными методами. В настоящей работе много времени было уделено реализации псевдоспектрального метода для полных уравнений, описанию которого посвящён следующий пункт, поэтому задачу о численном моделировании уравнений \eqref{eq_dim} было решено оставить на будущее.

\subsection{Полные трёхмерные уравнения}
%Разрабатывая численный метод для решения полных уравнений, мы задались целью сделать его пригодным для как можно более широкого круга численных экспериментов. Во-первых, метод должен подходить для моделирования эволюции некоторого начального импульса в течение длительного времени, для чего необходимо иметь возможность задачу на периодическом интервале. Во-вторых, метод должен позволять моделировать волны, возникающие в результате внешнего воздействия, как на боковую поверхность, так и на торец стержня (решение начально-краевой задачи на некотором конечном, но достаточно длинном интервале). %В третьих, мы хотели бы учитывать в будущем не только осесимметричную деформацию, т.е. случай, когда перемещение зависит от угловой координаты $\varphi$.

%Для того, чтобы выполнить поставленные цели, мы реализовали метод с дискретизацией по всем трём координатам. Для возможности моделирования уравнений в длинных стержнях с высокой точностью без рассмотрения гармоник или полиномов большой степени, мы реализовали многодоменный метод, позволяющий разбить область по осевой координате на несколько подобластей (доменов). Например, если нужно иметь 1000 точек коллокации по оси $x$, то, чтобы избежать учёт полиномов 1000-й степени, можно разбить область на 50 подобластей и на каждой подобласти использовать полиномы до 20-й степени. {\color{red} Поправить вводную часть после описания экспериментов}

Сформулируем начально-краевую задачу для полных уравнений, описывающих динамику однородного стержня круглого сечения длиной $L$ и радиусом $R$:
\begin{align}
\label{3_eqns_full}
&\ddot{\vect{U}}(x, r, \varphi, t) = %\frac1\rho
\rho^{-1} \,\divg\tens{P}, \qquad 0<x<L, \quad 0<r<R,\quad 0<\varphi<2\pi,\\
\label{3_bc_full}
&\tens{P} \cdot \vect{n} = \vect{P_b}, \quad r = R,\\
\label{3_iv_full}
&\vect{U}(x,r,\varphi,0) = \vect{U}_0(x,r,\varphi), \quad \vect{\dot U}(x,r,\varphi,0) = \vect{\dot U}_0(x,r,\varphi),\\
\label{3_stress}
\begin{split}
&\tens{P} = \lb \tens{I} + \nabla\vect{U} \rb \cdot \big[ \lambda \lb\trace{\tens{E}}\rb \tens{I} + 2\mu \tens{E} \\
& \hspace{30mm}+ l \lb\trace{\tens{E}}\rb^2 \tens{I} - m\lb \lb\trace{\tens{E}}\rb^2\tens{I} - 2\lb\trace{\tens{E}}\rb\tens{E} - \lb\trace{\tens{E}^2}\rb\tens{I}\rb + n\lb\tens{E}^*\rb^T\big],
\end{split}\\
\label{3_strain}
&\tens{E} = \frac12 \lb \nabla\vect{U}^T + \nabla\vect{U} + \nabla\vect{U}^T \cdot \nabla\vect{U} \rb,
\end{align}
где $\tens{E}^*$ --- союзная к $\tens{E}$ матрица, $\tens{I}$ --- единичный тензор, а по $x$ ставятся периодические граничные условия. Напомним, что уравнения \eqref{3_eqns_full} совпадают с \eqref{2_eq1_0}, \eqref{2_eq2_0} при отсутствии кручения и независимости перемещений от $\varphi$.

Будем искать решение в некотором конечномерном пространстве, задаваемым системой базисных функций $\Phi_{n_1}(x)$, $\Psi_{n_2}(r)$, $\Theta_{n_3}(\varphi)$. 
Представим приближённое решение $\vect{U}_{N}$ в спектральном виде и потребуем, чтобы приближённое решение совпадало с точным в точках коллокации:
\begin{align}
\label{3_spect_approx_full}
&\vect{U}_{N}(x,r,\varphi,t) = \sum_{n_1,n_2,n_3} \vect{\widehat{U}}(n_1,n_2,n_3,t) \Phi_{n_1}(x) \Psi_{n_2}(r) \Theta_{n_3}(\varphi),\\
\label{3_pseudo_sp_approx_full}
\begin{split}
&\vect{U}_{N}(x_{j_1}, r_{j_2}, \varphi_{j_3}, t) = \vect{U}(x_{j_1}, r_{j_2}, \varphi_{j_3}, t), \quad j_1 = 0,\dots, N_1-1, \quad j_2 = 0,\dots, N_2-1,\\
&\hspace{72.5mm} j_3 = 0,\dots, N_3-1.
\end{split}
\end{align}
В уравнении \eqref{3_eqns_full} отсутствуют смешанные производные по временной и пространственным переменным, поэтому подстановка \eqref{3_spect_approx_full} в \eqref{3_eqns_full} приводит к системе ОДУ второго порядка относительно~$\vect{U}$:
\begin{equation}\label{3_ode_inner}
\vect{\ddot{U}}_{N,k} = F_{k}(\vect{U}_{N,k_0}, \dots, \vect{U}_{N,k_N}), \quad k = k_0, \dots, k_N.
\end{equation}
Здесь $k$ --- мультииндекс: $k = (j_1,j_2,j_3)$, $k_0=(0,0,0)$, $k_N = (N_1-1,N_2-1,N_3-1)$, а также использовано обозначение: $\vect{\widehat U}_{N,j_1,j_2,j_3} = \vect{\widehat U}_{N}(x_{j_1}, r_{j_2}, \varphi_{j_3}, t)$. В граничных узлах необходимо сделать добавку, чтобы выполнялись граничные условия \eqref{3_bc_full}:
\begin{equation}\label{3_ode_boundary}
\vect{\ddot{U}}_{N,k} = F_{k}(\vect{U}_{N,k_0}, \dots, \vect{U}_{N,k_N}) + \vect{P}_{b,k} - \vect{P}(\vect{U}_{N,k_0}, \dots, \vect{U}_{N,k_N}), \quad r_k = R.
\end{equation}

Функция $F_k$ содержит в себе множество пространственных производных, которые можно вычислить по формуле~\eqref{3_phys_deriv} с помощью умножения на матрицу дифференцирования $D_N$. Оценим сложность этого метода на примере вычисления спектральной производной:
\begin{align} \label{3_spect_deriv_example}
&\lb\pdiff{U_N}{x}\rb_{j_1, j_2, j_3} = \lb D_{N_1}\rb_{j_1 l} U_{N,l,j_2,j_3},
\begin{split} 
j_1 = 0,\dots, N_1-1, \quad &j_2 = 0,\dots, N_2-1,\\
&j_3 = 0,\dots, N_3-1,
\end{split}
\end{align}
где в правой части выполняется суммирование по повторяющемуся индексу $l$. Такое вычисление производной требует $O(N_1^2N_2N_3)$ операций. Если в задаче используется базис Фурье или Чебышёва, то дифференцирование можно провести в пространстве коэффициентов $\widehat{U}_N$, для перехода к которым используется быстрое преобразование Фурье, что потребует $O\lb N_1\log(N_1) N_2 N_3\rb$ операций. Если используется другой базис, например, Лежандра, то такой возможности ускорить вычисление производной нет, однако практическая реализация алгоритма показала, что метод с использованием быстрого преобразования Фурье эффективнее умножения на матрицу лишь начиная со значения $N_1 \sim 1000$. Это обусловлено тем, что константа алгоритма умножения матриц значительно меньше константы быстрого преобразования Фурье.

Для того, чтобы  иметь возможность относительно быстро получать решения уравнений на сетке с большим количеством точек, мы применили \emph{многодоменный} псевдоспектральный метод, он же метод \emph{спектрального элемента}.
В рамках этого метода область задачи $\Omega$ разбивается на $M$ подобластей (доменов) $\Omega_m$:
\begin{align*}
&\Omega = \{(x,r,\varphi)\in\mathbb{R}^3 \, | \, 0<x<L, \; 0<r<R, \; 0<\varphi<2\pi\} = \bigcup\limits_{i=1}^{M} \Omega_m,\\
&\Omega_m = \{(x,r,\varphi)\in\mathbb{R}^3 \, | \, 0\leqslant L_{m}<x<L_{m+1}\leqslant L, \; 0<R_m<r<R_{m+1}\leqslant R, \\
&\hspace{40mm} 0\leqslant \varphi_m<\varphi<\varphi_{m+1}\leqslant 2\pi\}.
\end{align*}
На каждом домене осуществляется дискретизация согласно псевдоспектральному методу, а затем домены сшиваются так, как это происходит в методе конечных элементов. Теперь для вычисления производной \eqref{3_spect_deriv_example} требуется $O(M N_{1,M}^2 N_{2,M} N_{3,M})$, где $N_{1,M}$, $N_{2,M}$, $N_{3,M}$ --- размерность базиса на одном домене. Выбирая базисы не слишком большой размерности на каждом элементе, мы сможем добиться более высокой скорости вычисления пространственных производных.

Пример сетки с двумя доменами вдоль оси $x$ представлен на рисунке \ref{fig:grid}, где для дискретизации по $x$ и по $r$ используется базис Лежандра, причём по $x$ в качестве точек коллокации выбираются точки Гаусса-Лобатто, включающие границы, а по $r$ --- точки Гаусса-Радо, чтобы исключить центральную точку стержня ($r=0$). По $\varphi$ используется базис Фурье с равномерно распределёнными по интервалу $(0, 2\pi)$ точками.
Узлы такой сетки распределены по объёму неравномерно со сгущением около оси стержня и разрежением у его поверхности (рисунок \ref{fig:grid} слева). В работе~\cite{Boyd2011} предложен ряд способов как равномернее распределить узлы и, следовательно, повысить точность аппроксимации при том же количестве узлов, однако эти методы значительно сложнее и среди них нет наилучшего.

\begin{figure}[h]
	\centering
	\begin{subfigure}{0.5\textwidth}
		\centering
		\includegraphics[width=.8\textwidth]{Figures/Grid3D}
	\end{subfigure}
	\begin{subfigure}{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth]{Figures/Grid2D}
	\end{subfigure}
	\caption{Пример трёхмерной (слева) и двумерной (справа) сеток из двух доменов, обозначенных разными цветами.}
	\label{fig:grid}
\end{figure}

Для численного моделирования мы будем использовать многодоменный метод на сетке, являющейся прямым произведением одномерных сеток по каждой из осей, как это описано выше. Несмотря на отмеченный недостаток, такая сетка является значительно более простой в реализации. Более того, в настоящей работе рассматриваются решения уравнений лишь в осесимметричном случае, поэтому мы не проводим дискретизацию по $\varphi$, и, следовательно, проблемы неравномерного распределения узлов не возникает. Пример двумерной сетки по $x$ и $r$ приведён на рисунке \ref{fig:grid} справа.

Продемонстрируем ускорение, которое позволяет добиться многодоменный метод. Для этого будем решать полные уравнения в стержне радиусом $R=5$ и длиной $L=100$ с начальным условием в виде некоторой гладкой волны на сетках из 10 точек по $r$ и разным количеством точек по $x$ (от 10 до 70). Затем полученные решения сравним с решением, построенным на сетке $40\times 800$ точек по $r$ и $x$ соответственно. Результат представлен на рисунке \ref{fig:error}. При одном и том же значении относительной погрешности (примерно $1.5\cdot 10^{-4}$) метод с пятью доменами работает более, чем в три раза быстрее, чем однодоменный.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.48\linewidth]{Figures/ErrorBlack}
	\caption{Относительное отклонение решений по $L_2$-норме на грубых сетках от решения на мелкой сетке $40\times 800$ точек.}
	\label{fig:error}
\end{figure}



\section{Численное моделирование}

\subsection{Образование солитона из длинной волны}
Будем решать начально-краевую задачу для полных уравнений со свободными от напряжений границами ($\vect{P}_b=0$), где начальное условие задано в виде распространяющейся вправо со скоростью $c$ волны:
\begin{align}
\label{3_ic_u}
& U_0 (x, r) = A_0 W \erf\lb \frac{x - L/2}{W}\rb,&  &\dot{U}_0 (x, r) = -c \pdiff{U_0}{x}\\
\label{3_ic_v}
& V_0(x, r) = -\nu r \pdiff{U_0}{x},& &\dot{V}_0 (x, r) = -c \pdiff{V_0}{x}
\end{align}
На рисунке~\ref{fig:evolution1} слева изображена эволюция этой волны и зарождение солитона при $A_0=-0.01$, $W=20$ в стержне из полистирола, характеристики которого представлены в таблице~\ref{tab:ps}. Исходная волна испускает небольшие осцилляции и отрывается от них, что согласуется с моделью Буссинеска, согласно которой солитон распространяется со скоростью недоступной линейным волнам. На рисунке~\ref{fig:evolution1} справа представлено сравнение оторвавшейся в численном моделировании волны (<<численного>> солитона) с солитонном \eqref{soliton_v} из модели Буссинеска.
Если скорость солитона Буссинеска, являющуюся свободным параметром в \eqref{soliton_v}, подобрать так, чтобы она совпадала со скоростью численного солитона, то последний имеет почти ту же длину, что и солитон Буссинеска, но на 6\% большую амплитуду.
Отметим, что формула \eqref{soliton_v} задаёт четыре солитона различной длины, однако при амплитуде около $0.01$ они практически идентичны.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.48\linewidth]{Figures/SolEvolSmall}
	\includegraphics[width=0.48\linewidth]{Figures/SolCompareSmall}
	\caption{Слева: решение полных уравнений с начальными условиями \eqref{3_ic_u}, \eqref{3_ic_v} при $A_0=-0.01$, $W=20$. На графике изображены профили продольной деформации $-U_x(x - ct, 0, t)$ в центре стержня в различные моменты времени в движущейся со скоростью $c=\sqrt{E/\rho}$ системе координат. Масштаб деформации показан чёрным прямоугольником. 
	%На границах расчётной области задан поглощающий слой, чтобы волны, выходящие за левую границу области не входили справа.
	\\
	Справа: сравнение профиля солитона, полученного в численном моделировании, и солитона Буссинеска \eqref{soliton_v}, имеющего ту же скорость.}
	\label{fig:evolution1}
\end{figure}

На рисунке \ref{fig:evolution2} представлены результаты моделирования, аналогичного тому, что показано на рисунке \ref{fig:evolution1}, только с параметром $A_0=-0.05$. Исходная волна рождает два солитона, причём больший солитон движется быстрее меньшего. Сравнение большего солитона с регуляризованным солитоном Буссинеска показано на рисунке \ref{fig:evolution2} справа. При одинаковых скоростях численный солитон значительно длиннее и имеет примерно на 45\% большую амплитуду, чем солитон Буссинеска с такой же скоростью.

Как отмечалось в главе 2, деформации порядка $0.01$ сопоставимы с пределом упругости. Если решать задачу с параметром $A_0 = -0.001$, что соответствует порядку деформаций, возникающих в экспериментах, описанных в~\cite{Garbuzov}, то получающийся в результате солитон будет иметь всего на 0.3\% большую амплитуду, чем солитон Буссинеска с той же скоростью.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.48\linewidth]{Figures/SolEvolSmall2}
	\includegraphics[width=0.48\linewidth]{Figures/SolCompareSmall2}
	\caption{Слева: решение полных уравнений при $A_0=-0.05$, $W=20$. Справа: сравнение профиля большего численного солитона и регуляризованного солитона Буссинеска \eqref{soliton_v}, \eqref{B4mod}.}
	\label{fig:evolution2}
\end{figure}

Теперь решим аналогичную задачу для регуляризованного уравнения типа Буссинеска, где начальное условие задано формулой:
\begin{align}
\label{3_bq_ic}
& u_0 (x) = 2A_0 \exp\lb \frac{(x - L/2)^2}{W^2}\rb,&  &\dot{u}_0 (x) = -c \pdiff{u_0}{x}.
\end{align}
Отметим, что уравнение Буссинеска записано относительно функции $u = U_{x}$, а выражение поперечного перемещения $V$ через продольное $U$ уже <<вшито>> в модель.

На рисунке~\ref{fig:evol_compare} слева представлено сравнение эволюции волны согласно полным уравнениям и согласно модели Буссинеска.
В первом случае ($A_0=-0.01$) обе модели показывают качественно одинаковый результат, однако модель Буссинеска дала солитон примерно на 8\% большей амплитуды и, как следствие, бегущий быстрее солитона в полной модели. Во втором случае ($A_0=-0.05$) в модели Буссинеска родилось три солитона, а в полной модели два. Меньший солитон в обеих моделях идентичен, а больший солитон в модели Буссинеска на 26\% выше, чем солитон в полной модели. В случае малой амплитуды ($A_0=-0.001$) видимого различия между полной моделью и моделью Буссинеска не наблюдается, а разница амплитуды появляющихся солитонов составляет примерно 0.5\%.
\begin{figure}[h]
	\centering
	\includegraphics[width=0.49\linewidth]{Figures/SolEvolCompareSmallBlack}
	\includegraphics[width=0.49\linewidth]{Figures/SolEvolCompareSmallBlack2}
	\caption{Сравнение решения регуляризованного уравнения Буссинеска с начальным  условием~\eqref{3_bq_ic} и полных уравнений с начальными  условиями~\eqref{3_ic_u}, \eqref{3_ic_v} с параметрами $A_0=-0.01$ (слева) и $A_0=-0.05$ (справа).}
	\label{fig:evol_compare}
\end{figure}

Отличительной особенностью нелинейных волн является зависимость скорости от амплитуды волны. 
На основе численных решений полных уравнений с начальными волнами разной амплитуды мы построили зависимость скорости от амплитуды образовавшегося солитона и сравнили её с зависимостью в модели Буссинеска, задаваемой формулой \eqref{soliton_v} (см. рисунок \ref{fig:vel_ampl}). Зависимость скорости солитона Буссинеска для небольших амплитуд почти линейна, в то время как в полной модели рост скорости солитона существенно замедляется с увеличением амплитуды. На рисунке \ref{fig:vel_ampl} отчётливо виден асимптотический характер модели Буссинеска: чем меньше амплитуда волны, использовавшаяся как малый параметр $\varepsilon$ при выводе в параграфе \ref{sec:derivation}, тем больше соответствие полной модели.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.55\linewidth]{Figures/VelAmplBlack}
	\caption{Сравнение зависимостей скорости от амплитуды в рамках полных уравнений и уравнения Буссинеска. Горизонтальная линия --- скорость линейных продольных волн $c$.}
	\label{fig:vel_ampl}
\end{figure}


%\subsection{Образование солитона из бегущего по поверхности напряжения}
%Выведенная модель типа Буссинеска отличается от полученных ранее учётом напряжения на поверхности. 
%\begin{figure}[h]
%	\centering
%	\includegraphics[width=0.49\linewidth]{Figures/SolGenForceCompare}
%	\includegraphics[width=0.49\linewidth]{Figures/SolGenForce}
%	\caption{}
%	\label{fig:gen_compar1}
%\end{figure}


\subsection{Образование солитона из удара по торцу стержня}
В предыдущем пункте моделирование волн проводилось в стержне со свободной от напряжений поверхностью, однако представляет интерес изучение волн, вызываемых внешним воздействием, например, ударом по торцу стержня. Эти результаты могут быть полезны для того, чтобы возбудить уединённые волны в эксперименте. 

Модель типа Буссинеска строилась для бесконечного стержня, следовательно, вообще говоря, она применима для описания волн лишь вдали от торцов стержня. Моделирование волн вблизи торцов возможно с помощью полных уравнений с нулевыми начальными условиями (отсутствие возмущения) и следующими условиями на левой и правой границах стержня:
\begin{align}
\label{3_impact_stress}
&P_{xx} = 0.2 \cosh^{-2}\lb\frac t 2\rb\mbox{ГПа}, \qquad P_{rx} = 0&  &\mbox{при}\quad x = 0,\\
\label{3_impact_stress2}
&P_{xx} = P_{rx} = 0&  &\mbox{при}\quad x = L
\end{align}
%\begin{figure}[h]
%	\centering
%	\vspace{-10mm}
%	\includegraphics[width=.9\linewidth]{Figures/Pxx}
%	\caption{Схема удара. Зависимость нормального напряжения от времени на торце стержня.}
%	\label{fig:p_xx}
%\end{figure}

Чтобы иметь возможность сравнить результаты полной модели с моделью Буссинеска, было выбрано начальное условие в виде узкой волны, похожей по амплитуде и длине на ту, что образовывается в первый момент времени в полной модели. Для этого мы использовали начальное условие в виде:
\begin{align}
\label{3_bq_ic_2}
& u_0 (x) = -0.072 \exp\lb \frac{x^2}{5^2}\rb,&  &\dot{u}_0 (x) = -c \pdiff{u_0}{x}.
\end{align}

Сравнение решения полных уравнений с решением уравнения типа Буссинеска представлено рисунке \ref{fig:impact}. Быстрое воздействие \eqref{3_impact_stress} формирует короткую волну, половину которой мы видим в момент времени $t=0$. Затем эта волна, уменьшившись в амплитуде и породив множество мелких осцилляций, формирует солитон, отрывающийся от них.
\begin{figure}[h]
	\centering
	\includegraphics[width=\linewidth]{Figures/Impact_mod}
	\caption{Моделирование удара по торцу. Слева: решение полных уравнений с граничными условиями \eqref{3_impact_stress}, \eqref{3_impact_stress2}. Справа: решение резуляризованного уравнения Буссинеска с начальным условием \eqref{3_bq_ic_2}.}
	\label{fig:impact}
\end{figure}

Обе модели дают похожий результат лишь с той разницей, что полная модель учитывает большее количество различных мелких вибраций, порождаемые возникшей волной. Однако по скорости и форме родившиеся в обеих моделях солитоны достаточно близки: уравнение Буссинеска даёт солитон на 15\% большей амплитуды, при том что величина возмущений большая --- порядка $0.05$.


\subsection{Образование солитона из удара по боковой поверхности}
Представляет интерес возможность возбуждения солитона с помощью воздействия не только на его торец, но и на боковою поверхность. Для этого зададим непродолжительное локализованное на небольшом участке в центре стержня нормальное напряжение:
\begin{equation}\label{3_surf_imp}
P_{rr} = 0.3 \cosh^{-2}\lb\frac{t}{10}\rb \cosh^{-2}\lb\frac{x}{10}\rb\mbox{ГПа}, \quad P_{xr} = 0, \quad \mbox{при} \quad r=R.
\end{equation}
Напомним, что в уравнении Буссинеска \eqref{3_bq_reg} нормальное граничное напряжение $P_{rr}$ обозначено как функция $P(x,t)$, а касательное напряжение $P_{xr}$ как $T(x,t)$.

Результаты моделирования согласно полной модели и согласно модели Буссинеска представлены на рисунке \ref{fig:surf_impact}. Моделирование проводилось на отрезке $-2000<x<2000$, однако в силу симметричности результата относительно нуля мы показываем только правую часть этого промежутка. 
\begin{figure}[h]
	\centering
	\includegraphics[width=\linewidth]{Figures/SurfImpact}
	\caption{Моделирование удара по поверхности \eqref{3_surf_imp} в полных уравнениях и в регуляризованной модели Буссинеска.}
	\label{fig:surf_impact}
\end{figure}

Непродолжительное воздействие на небольшую область стержня около $x=0$ формирует две волны сжатия, бегущих в разные стороны. Эта волна затемпотом разделяется на мелкие осцилляции и солитон, отрывающийся от них. Полная модель предсказывает появление солитона амплитудой на 9\% больше, чем в модели Буссинеска.

Ещё раз подчеркнём, что деформации порядка $10^{-2}$ достаточно велики. Для меньших деформаций разница между моделями ещё меньше, из чего можно сделать вывод о применимости уравнения типа Буссинеска для моделирования генерации солитонов из коротких волн, возникающих вследствие удара.




\chapter*{Заключение}
\addcontentsline{toc}{chapter}{Заключение}

В настоящей работе получены следующие результаты:
\begin{enumerate}
	\item Для описания продольных длинных волн в стержнях круглого сечения выведены две новые асимптотические модели типа Буссинеска, отличающиеся друг от друга коэффициентами при дисперсионных слагаемых. Эти модели обобщены на случай ненулевой осесимметричной нагрузки на боковой поверхности, а также на случай предварительно растянутого стержня. Установлено, что новые модели:
	\begin{itemize}[noitemsep,topsep=1pt]
		\item имеют семейство решений в виде солитонов деформации;
		\item асимптотически эквивалентны моделям, выведенным ранее;
		\item обладают лучшими дисперсионными свойствами, чем модели, выведенные ранее.
	\end{itemize}

	\item Построен метод, позволяющий численно моделировать полные трёхмерные уравнения движения стержня в рамках нелинейной теории упругости. Проведена серия численных экспериментов, в ходе которых:
	\begin{itemize}[noitemsep,topsep=1pt]
		\item показано, что в полной модели, как и в модели Буссинеска, формируются солитоны;
		\item построена зависимость скорости от амплитуды солитона в полной модели и показано, что чем меньше амплитуда солитона, тем более схожи его свойства в полной модели и модели Буссинеска;
		\item установлено, что солитоны деформации могут быть возбуждены с помощью воздействия на поверхность стержня (удара), причём модель Буссинеска даёт схожие с полной моделью солитоны.
	\end{itemize}
\end{enumerate}

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

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

\newpage
\section*{Список публикаций}
\addcontentsline{toc}{chapter}{Список публикаций}
\begin{enumerate}[{label=[\arabic{*}]}]
	%\parindent0em 
	\item Garbuzov F.\,E., Khusnutdinova K.\,R., Semenova I.\,V., On Boussinesq-type models for long longitudinal waves in elastic rods, Wave Motion, 2019, 88, 129--143.
	\item Samsonov A.\,M., Semenova I.\,V., Garbuzov F.\,E.,  Nonlinear guided bulk waves in heteroge-neous elastic structural elements, Int. J. Nonlin. Mech., 2017, 94, 343--350.
	\item Semenova I.\,V., Belashov A.\,V., Garbuzov F.\,E., Samsonov A.\,M., Semenov A.\,A., Bulk strain solitons as a tool for determination of the third order elastic moduli of composite materials, Proceedings of SPIE 10329, 2017, 103291W.
	\item Samsonov A.\,M., Belashov A.\,V., Garbuzov F.\,E., Gula I.\,A., Semenov A.\,A., Semenova I.\,V., Solitary strain waves in a nanostructured rod, Proc. Intern. Conf. Days on Diffraction, 2017, 288--293.
	\item Гарбузов Ф.\,Е., Самсонов А.\,М., Семёнов А.\,А., Шварц А.\,Г., Определение упругих модулей	3-го порядка по параметрам объёмных солитонов деформации, ПЖТФ, 2016, 42\,(2), 121--123.
\end{enumerate}



\begin{thebibliography}{99}
	\addcontentsline{toc}{chapter}{Список источников}
	% Nonlinear waves: history review and general books 
	\bibitem{Dodd} Додд Р., Эйлбек Дж., Гиббон Дж., Моррис Х. Солитоны и нелинейные волновые уравнения. М.: «Мир», 1988.
	%\bibitem{Ablowitz} Абловиц М., Сигур Х. Солитоны и метод обратной задачи. М.: Мир, 1987.
	\bibitem{Ablowitz2011} Ablowitz M.\,J., Nonlinear dispersive waves: asymptotic analysis and solitons, Cambridge: Cambridge University Press, 2011.
	%\bibitem{Russel} Scott Russell J. \textit{Report on waves}, Rept. 14th meetings of the British Assoc. for the Advancement of Science. London: John Murray, 1844. P. 311–390.
	%\bibitem{Boussinesq} Boussinesq J. The?orie de l’intumescence liquid apple?e onde solitarie ou de translation se propagent dans un canal rectangularie // Comptes Rendus. 1871. Vol. 72. P. 755–759.
	%\bibitem{Korteweg} Korteweg D.J., de Vries G. On the change of form of long waves advancing in a rectangular channel, and on a new type of long stationary waves // Phil. Mag. 1895. Vol. 39. P. 422–443.
	%\bibitem{BBM} T.B. Benjamin, J.L. Bona, and J.J. Mahony, Model equations for long waves in nonlinear dispersive systems, \textit{Philos. Trans. R. Soc. London, Ser. A} 272(1220) (1972) 47-78.
	\bibitem{Zabusky} Zabusky N.J., Kruskal M.D. Interaction of «solitons» in a collisionless plasma and the reccurence of initial states, Physical Review Letters, 1965, 15 (6), 240--243.
	\bibitem{GardnerIST} Gardner C.\,S., Greene J.\,M., Kruskal M.\,D., Miura R.\,M., Method for Solving the Korteweg-de Vries Equation, Physical Review Letters, 1967, 19 (19), 1095--1097.
	\bibitem{Zakharov} Захаров В.\,Е., Манаков С.\,В., Новиков С.\,П., Питаевский Л.\,П., Теория солитонов: метод обратной задачи, М.: <<Наука>>, 1980. 
	
	% Numerical methods review
	\bibitem{Taha} Taha T.\,R., Ablowitz M.\,I., Analytical and numerical aspects of certain nonlinear evolution equations. III. Numerical, Korteweg-de Vries equation, Journal of Computational Physics, 1984, 55 (2), 231--253.
	\bibitem{Christov} Christov C.\,I., Conservative Difference Scheme for Boussinesq Model of Surface Waves, Proceedings ICFD 5, 1996, 343--349.
	\bibitem{Kolkovska} Kolkovska N., Dimova M., A new conservative finite difference scheme for Boussinesq paradigm equation, Central European Journal of Mathematics, 2012, 10 (3), 1159--1171.
	%\bibitem{Dutykh_KdV} Dutykh D., Katsaounis T., Mitsotakis D., Finite volume methods for unidirectional dispersive wave models, электронная публикация, arXiv/HAL
	\bibitem{Dutykh_Boussinesq} Dutykh D., Katsaounis T., Mitsotakis D., Finite volume schemes for Boussinesq type equations, электронная публикация, https://arxiv.org/abs/1101.1728
	\bibitem{Karczewska} Karczewska A., Rozmej P., Szczecinski M., Boguniewicz B., A finite element method for extended KdV equations, International Journal of Applied Mathematics and Computer Science, 2016, 26 (3), 555--567.
	\bibitem{Gottlieb_Orszag} Gottlieb D., Orszag S.\,A., Numerical Analysis of Spectral Methods: Theory and Applications, Philadelphia: SIAM, 1977.
	%\bibitem{Trefethen} Trefethen L.\,N., Spectral Methods in MATLAB, SIAM, Philadelphia, 2000.
	\bibitem{Canuto2007} Canuto C. et al., Spectral Methods. Evolution to Complex Geomenties and Applications to Fluid Dynamics, Berlin: Springer-Verlag, 2007.
	
	% Nonlinear mechanics %
	\bibitem{LurieNL} Лурье А.\,И., Нелинейная теория упругости, М.: <<Наука>>, 1980.
	\bibitem{Murnaghan} Murnaghan F.\,D., Finite deformation of an elastic solid, New York: John Wiley and Sons, 1951.
	\bibitem{Bergstrom} Bergstr\"{o}m J., Mechanics of Solid Polymers, Amsterdam: William Andrew Publishing, 2015.
	\bibitem{Yu} Yu Y.-Y., Generalized Hamilton's Principle and Variational Equation of Motion in Nonlinear Elasticity Theory, With Application to Plate Theory, Journal of the Acoustical Society of America, 1964, 36, 111--120.
	
	% Nonlinear waves
	%\bibitem{Maugin} G.A. Maugin, \textit{Nonlinear waves in elastic crystals}, Oxford University Press, Oxford, 1999.
	%\bibitem{Dai} Dai~H.-H., Z. Cai, Phase transition in a slender cylinder composed of an incompressible elastic material. I. Asymptotic model equation, \textit{Proc. Roy. Soc. A} 462 (2006) 419-438.
	%\bibitem{M}
	%A. Mayer, 
	%Surface acoustic waves in nonlinear elastic media, \textit{Phys. Reports} 256 (1995) 257 - ?.
	%Nonlinear surface acoustic waves: Theory, \textit{Ultrasonics} 48 (2008) 478-481.
	%\bibitem{HL} P. Hess, A.M. Lomonosov, Solitary surface acoustic waves and bulk solitons in nanosecond and picosecond laser ultrasonics, \textit{Ultrasonics} 50 (2010) 167-171.
	%\bibitem{E1} J. Engelbrecht, A. Salupere and K. Tamm, Waves in microstructured solids and the Boussinesq paradigm, \textit{Wave Motion} 48 (2011) 717-726.
	%\bibitem{P} A. Pau, F. Lanza di Scalea, Nonlinear guided wave propagation in prestressed plates, \textit{J. Acoust. Soc. Am.} 137 (2015) 1529-1540.
	%\bibitem{E2} T. Peets, K. Tamm, J. Engelbrecht, On the role of nonlinearities in the Boussinesq-type wave equations, \textit{Wave Motion} 71 (2017) 113-119.
	\bibitem{S_book} Samsonov A.\,M., Strain solitons in solids and how to construct them, Boca Raton: Chapman \& Hall/CRC, 2001.
	\bibitem{P_book} Porubov A.\,V., Amplification of nonlinear strain waves in solids, Singapore: World Scientific, 2003.
	\bibitem{NS} Nariboli G.\,A., Sedov A., Burgers-Korteweg de Vries equation for viscoelastic rods and plates, J. Math. Anal. Appl., 1970, 32 (3), 661--677.
	\bibitem{OS} Ostrovsky~L.\,A., Sutin~A.\,M., Nonlinear elastic waves in rods, PMM, 1977, 41, 531--537.
	\bibitem{Love} Love A\,.E.\,H., A treatise on the mathematical theory of elasticity, London: Cambridge University Press, 1927.
	%\bibitem{S1} Samsonov~A.\,M., Structural optimization in nonlinear wave propagation problems. In: \textit{Structural Optimization under Dynamical Loading. Seminar and Workshop for Junior Scientists}, U. Lepik ed., Tartu University Press, 75-76 (1982).
	%\bibitem{S2} Samsonov~A.\,M., Soliton evolution in a rod with variable cross section, \textit{Sov. Physics - Doklady} 29 (1984) 586-587.
	\bibitem{SP} Порубов А\,.В., Самсонов А.\,М. Письма в ЖТФ, 1993, 19(12), 26--29.
	\bibitem{PV}
	Porubov~A.\,V., M.G. Velarde, Dispersive - dissipative solitons in nonllinear solids, Wave Motion, 2000, 31(3), 197--207.
	\bibitem{E_book} Ерофеев В.\,И., Кажаев В.\,В., Семерикова Н.\,П., Волны в	стержнях. Дисперсия. Диссипация. Нелинейность. М.: Физматлит, 2002. -- 208 с.
	\bibitem{DF} Dai~H.-H., Fan~X., Asymptotically approximate model equations for weakly nonlinear long waves in compressible elastic rods and their comparisons with other simplified model equations, Maths. Mechs. Solids, 2004, 9, 61--79.
	\bibitem{DC} Dai~H.-H., and Z. Cai, Uniform asymptotic analysis for transient waves in a pre-stressed compressible hyperelastic rod, Acta Mechanica, 2000, 139, 201--230.
	\bibitem{KSZ} Khusnutdinova~K.\,R., Samsonov~A.\,M., A.S. Zakharov, Nonlinear layered lattice model and generalized solitary waves in imperfectly bonded structures, Phys. Rev. E, 2009, 79(5), 056606.
	\bibitem{KS} Khusnutdinova~K.\,R., Samsonov~A.\,M., Fission of a longitudinal strain solitary wave in a delaminated bar, Phys. Rev. E, 2008, 77, 066603.
	\bibitem{KT1} Khusnutdinova~K.\,R., Tranter~M.\,R., Modelling of nonlinear wave scattering in a delaminated elastic bar, Proc. R. Soc. A, 2015, 471, 20150584.
	\bibitem{KT2} Khusnutdinova~K.\,R., Tranter~M.\,R., On radiating solitary waves in bi-layers with delamination and coupled Ostrovsky equations, Chaos, 2017, 27, 013112.
	\bibitem{JAP2010} Dreiden~G.\,V., Khusnutdinova~K.\,R., Samsonov~A.\,M., and Semenova~I.\,V., Splitting induced generation of soliton trains in layered waveguides, J. Appl. Phys., 2010, 107, 034909.
	\bibitem{JAP2012} Dreiden G.\,V., Khusnutdinova~K.\,R., Samsonov~A.\,M., and Semenova~I.\,V., Bulk strain solitary waves in bonded layered polymeric bars with delamination, J. Appl. Phys., 2012, 112, 063516.
	%\bibitem{ML} J.-F. Mercier, B. Lombard, A two-way model for nonlinear acoustic waves in a non-unfirm lattice of Helmholtz resonators, \textit{Wave Motion} 72 (2017) 260-275.
	%\bibitem{AM} R. Arredondo and J.P. McHugh, Mean displacement near an interface in a nonlinear string, \textit{SIAM J. Appl. Math.} 78 (2018) 1470-1488.
	%\bibitem{DS} M. Destrade and G. Saccomandi, Nonlinear transverse waves in deformed dispersive solids, \textit{Wave Motion} 45 (2008) 325-336.
	\bibitem{Garbuzov} Garbuzov F.\,E., Khusnutdinova K.\,R., Semenova I.\,V., On Boussinesq-type models for long longitudinal waves in elastic rods, Wave Motion, 2019, 88 129--143.
	
	%\bibitem{DSS} M. Destrade, G. Saccomandi, I. Sigura, Methodical fitting for mathematical models of rubber-like materials, \textit{Proc. R. Soc. A} 473 (2017) 20160811.
	%\bibitem{Ogden} R.W. Ogden, \textit{Non-linear elastic deformations}, Dover, Mineola, New York, 1984.
	\bibitem{bostrm2000} Bostr\"{o}m A., On wave equations for elastic rods, ZAMM, 2000, 80 (4), 245--251. 
	\bibitem{HughesKelly} Hughes D.\,S., Kelly J.\,L., Second order elastic deformation of solids, Phys. Rev., 1953, 92,  1145--1149.
	\bibitem{ADO} Abiza Z., Destrade M., Ogden R.\,W., Large acoustoelastic effect, Wave Motion, 2012, 49, 364--374.
	\bibitem{Boyd2011} Boyd J.\,P., Yu F., Comparing seven spectral methods for interpolation and for solving the Poisson equation in a disk: Zernike polynomials, Logan–Shepp ridge polynomials, Chebyshev–Fourier Series, cylindrical Robert functions, Bessel–Fourier expansions, square-to-disk conformal mapping and radial basis functions, Journal of Computational Physics, 2011, 230\,(4) 1408--1438.
%	\bibitem{PAT} A.V. Pichugin, H. Askes, A. Tyas, Asymptotic equivalence of homogenisation procedures and fine-tuning of continuum theories, \textit{J. Sound and Vibration} 313 92008) 858-874.
%	\bibitem{ADKM} I.V. Andrianov, V.D. Danishevsky, J.D. Kaplunov and B. Markert, Wide frequency higher-order dynamic model for transient waves in a lattice, In: I.V. Andrianov et al. ed., ``Problems of Nonlinear Mechanics and Physics of Materials", Springer, 2019.
%	\bibitem{TP1988} G. V. Dreiden, Yu. I. Ostrovsky, A. M. Samsonov, I. V. Semenova, E. V. Sokurinskaya, Formation and propagation of strain solitons in nonlinearly elastic solid, \textit{Techn. Phys.} 58 (1988) 2040-2017. (in Russian)
%	\bibitem{Dreiden-TP-2008} Dreiden~G.\,V., Samsonov~A.\,M., Semenova~I.\,V., Evolution of bulk strain solitons in long polymeric waveguides, \textit{Techn. Phys.} 53 (2008) 540-546.
%	\bibitem{TPL-Polycarb}  Dreiden~G.\,V., Samsonov~A.\,M., Semenova~I.\,V., Bulk elastic strain solitons in polycarbonate, \textit{Techn. Phys. Lett.} 37 (2011) 500-502.
%	\bibitem{JAP2008} Dreiden~G.\,V., Khusnutdinova~K.\,R., Samsonov~A.\,M., Semenova~I.\,V., Comparison of the effect of cyanoacrylate- and polyurethane-based adhesives on a longitudinal strain solitary wave in layered polymethylmethacrylate waveguides, \textit{J. Appl. Phys.} 104 (2008) 086106.
%	\bibitem{APL2014} Dreiden~G.\,V., Samsonov~A.\,M., Semenova~I.\,V., and A.G. Shvartz, Strain solitary waves in a thin-walled waveguide, \textit{Appl. Phys. Lett.} 105 (2014) 211906.
%	\bibitem{APL2018} A.V. Belashov, Y.M. Beltukov, N.V. Petrov, Samsonov~A.\,M., Semenova~I.\,V., Indirect assessment of bulk strain soliton velocity in opaque solids, \textit{Appl. Phys. Lett.} 112 (2018) 121903.
%	\bibitem{Strain2010} Dreiden~G.\,V., Khusnutdinova~K.\,R., Samsonov~A.\,M., Semenova~I.\,V., Longitudinal strain solitary wave in a two-layered polymeric bar, \textit{Strain} 46 (2010) 589-598.
%	\bibitem{SPIE2018} A.V. Belashov, Y.M. Beltukov, Semenova~I.\,V., Pump-probe digital holography for monitoring of long bulk nonlinear strain waves in solid waveguides, \textit{Proc. SPIE} 10678 (2018) 1067810.
\end{thebibliography}


\chapter*{\vspace{-25mm}Приложение 1. \\Длинные формулы}
\addcontentsline{toc}{chapter}{Приложение 1. Длинные формулы}
\vspace{-5mm}
Нелинейные функции в уравнении \eqref{eq2_1}:
\begin{equation} \nonumber
\begin{split}
&\begin{split}
\Phi_1 =& \, 2\left[(-4\lambda - 4\mu + n - 4m) V_1 - 2(\lambda + 2\mu + m) U_{0x}\right] U_2 \\
&- \left[ 2(2l + \lambda) V_1 + (3\lambda + 6\mu + 2l + 4m) U_{0x} \right] U_{0xx} \\
& - \left[ (2\lambda + 2\mu + 8l + n) V_1 + 2(\lambda + \mu + 2l + m) U_{0x} \right] V_{1x},
\end{split} \\
&\begin{split}
\Phi_2 =& \, \frac12 \left[2 (2\lambda + 2\mu + 8l + n) U_{2x} + (4\lambda + 4\mu + 4m - n) V_{1xx} + 32(2\lambda + 3\mu + 2l + 2m) V_3 \right] V_1 \\
& + 2(\lambda + \mu + 2l + m) U_{0x} U_{2x} + 2(\mu + m)\left[ U_{0xx} + 4 V_{1x} \right] U_2 + (\lambda + 2\mu + m)(U_{0x} V_{1x})_x\\
& + \frac14(12\lambda + 20\mu + 12m - n) V_{1x}^2 + 8(\lambda + 2l) V_3 U_{0x} + (4\lambda + 12\mu + 4m + n) U_2^2.
\end{split}
\end{split}
\end{equation}
Функции из формул \eqref{U2}, \eqref{V3}:
\begin{align}
\nonumber
\begin{split}
f_2(x, t) =& \, \frac{1}{8 \mu ^2} \bigg(U_{0xx} V_1 \left(8\mu^2 - 2\mu(-4\lambda + 4l - 4m + n) + \lambda(4\lambda + 4m - n)\right) \\
&+ 2V_1 V_{1x} (2 (\lambda + \mu) (2\lambda + \mu) - 8l\mu + 4m(\lambda + \mu) - n(\lambda + 2\mu)) \\
&+ \rho c^2 U_{0tt}V_1 (-4(\lambda + \mu) - 4m + n) + 2U_{0x}U_{0xx} \left(-2\mu^2 + \mu(\lambda - 2(l + m)) + \lambda (\lambda + m)\right) \\
&+ 4 U_{0x} V_{1x} \left((\lambda + \mu)^2 - 2l\mu + \lambda m\right) - 2\rho c^2 U_{0x}U_{0tt} (\lambda + 2\mu + m)\bigg)
\end{split}\\
\nonumber
\begin{split}
f_3(x, t) =& -\frac{1}{8 (\lambda +2 \mu )}\bigg(-\frac{2 V_1 (2 (\lambda +l+m)+3 \mu ) \left(-\rho  s^2 V_{1tt}+2 (\lambda +\mu ) U_{2,x}+\mu  V_{1xx}\right)}{\lambda +2 \mu }\\
&-\frac{(\lambda +2 l) U_{0x} \left(-\rho  s^2 V_{1tt}+2 (\lambda +\mu ) U_{2x}+\mu  V_{1xx}\right)}{\lambda +2 \mu }+8 l V_1 U_{2x} + (4l + 2m)U_{0x} U_{2x}\\
&+2 U_2 (\mu +m) \left(U_{0xx}+4 V_{1,x}\right)+m (U_{0x}V_{1x})_x + \lb 3\lambda + 5\mu + 3m -\frac{n}{4}\rb V_{1x}^2\\
&+\lb 2\lambda + 2\mu + 2m-\frac{n}{2}\rb V_1 V_{1xx}+(2\lambda + 2\mu + n) V_1 U_{2x} + \lambda(U_{0x}V_{1x})_x\\
&+ 2\mu(U_{0x} V_{1x})_x + (2\lambda + 2\mu)U_{0x}U_{2x} + (4 (\lambda +3 \mu )+4 m+n) U_2^2 \bigg)
\end{split}
\end{align}
Функции из асимптотического представления $V_1$~\eqref{v1_asympt}:
\begin{align}
\nonumber
\begin{split}
f(x,t) &= \frac{1}{16 (\lambda +\mu )^3 (\lambda +2 \mu )}\bigg(\mu  \rho  s^2 (3 \lambda +2 \mu ) (2 \lambda +3 \mu ) P_{tt}+\lambda  \mu  (\lambda +2 \mu ) (3 \lambda +2 \mu ) P_{xx}\\
&\hspace{25mm}+(\lambda +\mu ) \left(\mu  (\lambda +2 \mu ) (4 \lambda +3 \mu ) U_{0xxx}-\rho c^2 \left(3 \lambda ^2+7 \lambda  \mu +3 \mu ^2\right) U_{0xtt}\right)
\end{split}\\
\nonumber
\begin{split}
g(x,t) &= -\frac{1}{8(\lambda + \mu)}\bigg(-\frac{2U_{0x}(\lambda + 4l - 2m + n) \left(\lambda(\lambda + \mu) U_{0x} - \mu P (3\lambda + 2\mu)\right)}{(\lambda + \mu)^2} + \\
&\hspace{20mm}\frac{(3 (\lambda +\mu )+4 l+2 m) \left(\mu  P (3 \lambda +2 \mu )-\lambda  (\lambda +\mu ) U_{0,x}\right){}^2}{(\lambda +\mu )^4} + (4l + 2\lambda)  U_{0x}^2\bigg)
\end{split}
\end{align}

\chapter*{\vspace{-20mm}Приложение 2.\\Солитонные решения уравнения типа Буссинеска}
\addcontentsline{toc}{chapter}{Приложение 2. Солитонные решения уравнения типа Буссинеска}
Все уравнения типа Буссинеска, обсуждавшиеся в настоящей работе, могут быть записаны в виде:
\begin{equation}
u_{tt} - c^2 u_{xx} = \tilde{\beta} (u^2)_{xx} + \tilde\alpha_1 u_{tttt} + \tilde\alpha_2 u_{ttxx} + \tilde\alpha_3 u_{xxxx},
\end{equation}
где $c$, $\tilde\beta$ и $\tilde\alpha_i$ --- некоторые постоянные. Будем искать решения в виде волн, бегущих влево или вправо:$ u = u(\xi), \; \mbox{где} \; \xi = x \pm vt, $
что сводит исходное уравнение к обыкновенному дифференциальному уравнению:
\begin{equation}
(v^2 - c^2) u'' = \tilde\beta (u^2)'' + (\tilde\alpha_1 v^4 + \tilde\alpha_2 v^2 + \tilde\alpha_3) u^{IV}.
\label{ODE}
\end{equation}
Интегрируя это уравнения по $\xi$ дважды и требуя, чтобы на бесконечности не было возмущений $u, u', u'', u''' \to 0$ при $\xi \to \pm \infty$, получаем уравнение:
\begin{equation}
u'' = \frac{(v^2-c^2) u - \tilde\beta u^2}{\tilde\alpha_1 v^4 + \tilde\alpha_2 v^2 + \tilde\alpha_3},
\end{equation}
которое может рассматриваться как уравнение движения частицы единичной массы в поле потенциальной силы. Интеграл энергии имеет вид:
\begin{equation}
E = \frac 12 \left (u' \right )^2 - \frac{3 (v^2-c^2) u^2 - 2 \tilde\beta u^3}{6 (\tilde\alpha_1 v^4 + \tilde\alpha_2 v^2 + \tilde\alpha_3)},
\end{equation}
а солитонное решение соответствует нулевому уровню энергии $E=0$. Разделение переменных позволяет получить солитонное решение в виде:
\begin{equation} \label{4_soliton}
u = \frac{3 (v^2 - c^2)}{2 \tilde\beta} \sech^2 \left [\sqrt{\frac{v^2 - c^2}{4 (\tilde\alpha_1 v^4 + \tilde\alpha_2 v^2 + \tilde\alpha_3)}} (x \pm v t)\right ]
\end{equation}
для таких значений параметра $v$, что это решение вещественнозначное.
Решение \eqref{4_soliton} можно репараметризовать через амплитуду $A$:
\begin{equation}
u = A \sech^2 \left [\sqrt{\frac{3A\tilde{\beta}}{2 (\tilde\alpha_1 (3c^2+2A\beta)^2 + 3\tilde\alpha_2 (3c^2+2A\beta) + 9\tilde\alpha_3)}} \lb x \pm t\sqrt{c^2+\frac23A\beta} \rb\right ].
\end{equation}
	
\end{document}