Accounting for viscous and thermal effects in time in computational problems of acoustics

Введение

Как известно [1], термические и вязкие эффекты в воздухе малы, что позволяет описывать распространение акустических волн с помощью волнового уравнения. Однако в узких полостях и вблизи границ термические и вязкие потери могут быть значительны. Прямой способ учета термовязких эффектов заключается в точном или приближенном решении уравнений Навье-Стокса, как это делается, например, в [2–4].

Решение уравнений Навье-Стокса с вычислительной точки зрения представляет собой весьма трудоемкую задачу. Более эффективным является использование системы уравнений пограничного слоя вблизи границ и волнового уравнения в остальной области [5, 6]. Такой подход сопряжен с вычислительными сложностями – в пограничном слое необходимо генерировать более плотную сетку и осуществлять сращивание двух разных решений в граничных областях. Наиболее простой и эффективный метод заключается в учете термовязких эффектов с помощью граничного условия специального вида. В [1, 7] термовязкие граничные условия получены для стационарной задачи и имеют следующий вид:

$\begin{array}{c}{\delta }_V\frac{i-1}{2}{\Delta }_{\tau }\tilde{p}(\omega )-{\delta }_T\frac{{\omega }^2}{c_0^2}\frac{(i-1)(\gamma -1)}{2}\tilde{p}(\omega )+\\ +\frac{\partial \tilde{p}(\omega )}{\partial n}=\mathrm{0,}\end{array}$      (1)

где ω = 2πf – частота, c₀ – скорость звука, $\tilde{p}(\omega )$ – Фурье-образ давления,

${\delta }_V=\sqrt{\frac{2\nu }{\omega }}$, ${\delta }_T=\sqrt{\frac{2\kappa }{\omega {\rho }_0{c}_p}}$

представляют собой толщины соответственно вязкого и термического слоя, _ρ0 – равновесная плотность воздуха, γ = c_p/c_V – показатель адиабаты, ν – кинематическая вязкость, c_p – удельная теплоемкость при постоянном давлении, κ – коэффициент теплопроводности. а _Δτ – касательный лапласиан (оператор Лапласа–Бельтрами).

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

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

Постановка задачи

Пусть акустическое давление p(x, t) удовлетворяет в некоторой конечной области Ω волновому уравнению (см. рис. 1):

 $\begin{array}{cc}\frac{1}{{ñ}_0^2}\frac{{\partial }^{2}p(x,t)}{\partial t^2}-\Delta p(x,t)={\rho }_{0}f(x,t),& x\in \Omega \end{array}$,        (2)

где f(x, t) – функция источника, Δ – оператор Лапласа, x º (x, y, z) – координатный вектор.

Рис. 1. Примерный вид области Ω.

Разобьем границу области $\partial$ Ω на две части. Часть границы, на которой термовязкие эффекты пренебрежимо малы, обозначим как $\partial$ _Ωn. На ней выполняется граничное условие Неймана:

 $\begin{array}{cc}\frac{\partial p(x,t)}{\partial n}=\mathrm{0,}& x\in \partial {\Omega }_n\end{array}$,             (3)

где n – вектор внутренней нормали.

Часть границы, на которой термовязкие эффекты значительны, обозначим как $\partial$ _Ωv. Введем прямое и обратное преобразование Фурье для произвольной функции f(x, t)

$\begin{array}{c}\tilde{f}(x,\omega )={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}f(x,t)\exp (-i\omega t)dt},\\ f(x,t)=\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\tilde{f}(x,\omega )\exp (i\omega t)d\omega }.\end{array}$         (4)

Получение временного представления термовязких граничных условий путем формального выполнения обратного преобразования Фурье (4) от выражения (1) сопряжено с определенными трудностями, так как частота входит в граничное условие в дробной степени (первое слагаемое пропорционально ^ω-1/2, а второе – ^ω3/2). В результате во временной постановке термовязких граничных условий появляются дробные производные (см. приложение А):

$\begin{array}{c}\sqrt{\nu }{}_0{}^{C}D{}_t^{-1/2}{\Delta }_{\tau }p(x,t)+\frac{\gamma -1}{c_0^2}\sqrt{\frac{\kappa }{{\rho }_0{c}_p}}{}_0{}^{C}D{}_t^{3/2}p(x,t)+\\ +\frac{\partial p(x,t)}{\partial n}=\mathrm{0,}x\in \partial {\Omega }_v.\end{array}$       (5)

Здесь введено обозначение дробной производной в смысле Капуто:

${}_a{}^{C}D{}_t^{\beta }f(t)=\frac{1}{\Gamma (n-\beta )}{\displaystyle \underset{a}{\overset{t}{\int }}\frac{d^{n}f(\tau )}{d{\tau }^n}\frac{1}{{(t-\tau )}^{\beta +1-n}}d\tau }$,

где n – ближайшая к β целая степень, такая что n > β.

Задачу необходимо дополнить начальными условиями:

$\begin{array}{cc}{p}_0=p(x,0),& {p^{\prime }}_0=\frac{dp}{dt}(x,0).\end{array}$

Отметим, что условие (5), учитывающее термовязкие потери, в настоящей работе используется впервые.

Конечно-элементная формулировка

Пусть имеется некоторое пространство гладких весовых функций H. Рассмотрим слабую постановку задачи (2) для всех w $\in$ H:

$\frac{1}{c_0^2}{\displaystyle \underset{\Omega }{\int }\frac{{\partial }^{2}p(x,t)}{\partial t^2}wdV}-{\displaystyle \underset{\Omega }{\int }\Delta p(x,t)wdV}={\rho }_0{\displaystyle \underset{\Omega }{\int }f(x,t)wdV}.$

Применим первую формулу Грина ко второму члену в левой части

${\displaystyle \underset{\Omega }{\int }\Delta p(x,t)wdV}=-{\displaystyle \underset{\Omega }{\int }\nabla p(x,t)\nabla wdV}+{\displaystyle \underset{\partial \Omega }{\int }\frac{\partial p(x,t)}{\partial n}wdS}$.

Пользуясь (5), преобразуем член с нормальной производной:

$\begin{array}{c}{\displaystyle \underset{\partial \Omega }{\int }\frac{\partial p(x,t)}{\partial n}wdS}=-\sqrt{\nu }{\displaystyle \underset{\partial {\Omega }_{\nu }}{\int }{}_0{}^{C}D{}_t^{-1/2}{\Delta }_{\tau }p(x,t)wdS}-\\ -\frac{\gamma -1}{c_0^2}\sqrt{\frac{\kappa }{{\rho }_0{c}_p}}{\displaystyle \underset{\partial {\Omega }_{\nu }}{\int }{}_0{}^{C}D{}_t^{3/2}p(x,t)wdS}=\\ =\sqrt{\nu }{\displaystyle \underset{\partial {\Omega }_{\nu }}{\int }{}_0{}^{C}D{}_t^{-1/2}{\nabla }_{\tau }p(x,t){\nabla }_{\tau }wdS}-\\ -\frac{\gamma -1}{c_0^2}\sqrt{\frac{\kappa }{{\rho }_0{c}_p}}{\displaystyle \underset{\partial {\Omega }_{\nu }}{\int }{}_0{}^{C}D{}_t^{3/2}p(x,t)wdS}.\end{array}$

Тогда получаем:

$\begin{array}{c}\frac{1}{c_0^2}{\displaystyle \underset{\Omega }{\int }\frac{{\partial }^{2}p(x,t)}{\partial t^2}wdV}+{\displaystyle \underset{\Omega }{\int }\nabla p(x,t)\nabla wdV}-\\ -\sqrt{\nu }{\displaystyle \underset{\partial {\Omega }_{\nu }}{\int }{}_0{}^{C}D{}_t^{-1/2}{\nabla }_{\tau }p(x,t){\nabla }_{\tau }wdS}+\\ +\frac{\gamma -1}{c_0^2}\sqrt{\frac{\kappa }{{\rho }_0{c}_p}}{\displaystyle \underset{\partial {\Omega }_{\nu }}{\int }{}_0{}^{C}D{}_t^{3/2}p(x,t)wdS}={\rho }_0{\displaystyle \underset{\Omega }{\int }f(x,t)wdV}.\end{array}$

Последнее выражение представляет собой выражение для волнового уравнения в слабой постановке. В соответствии с идеей метода конечных элементов, аппроксимируем давление следующим образом [8]:

$p(x,t)={\displaystyle \sum _{k=1}^L{p}^k(t)N_k(x)}$,

где p^k(t) – значения давления в узлах, N_k(x) – так называемые функции формы, L – число узлов дискретизации. Учитывая граничные условия (3), (5) и используя функции формы в качестве весовых функций w (метод Галеркина), получим конечно-элементную формулировку:

$M\frac{d^2}{d{t}^2}p+V{}_0{}^{C}D{}_t^{-1/2}p+T{}_0{}^{C}D{}_t^{3/2}p+Kp=f$,           (6)

Где

$\begin{array}{c}\begin{array}{cc}M=\frac{1}{{\rho }_0{c}_0^2}{\displaystyle \underset{\Omega }{\int }N{N}^{T}dV},& K=\frac{1}{{\rho }_0}{\displaystyle \underset{\Omega }{\int }\nabla N\nabla N^{T}dV},\end{array}\\ \begin{array}{l}\begin{array}{cc}V=-\frac{\sqrt{\nu }}{{\rho }_0}{\displaystyle \underset{\partial {\Omega }_{\nu }}{\int }{\nabla }_{\tau }N{\nabla }_{\tau }{N}^{T}dS},& \end{array}\\ T=\frac{\gamma -1}{{\rho }_0{c}_0^2}\sqrt{\frac{\kappa }{{\rho }_0{c}_p}}{\displaystyle \underset{\partial {\Omega }_{\nu }}{\int }N{N}^{T}dS,}\end{array}\\ f={\displaystyle \underset{\Omega }{\int }f\text{}NdV}.\end{array}$

Здесь p – столбец давлений:

$p={\left(\begin{array}{cccc}{p}^1& p^2& \cdots & p^L\end{array}\right)}^T$,

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

$N={\left(\begin{array}{cccc}{N}_1& N_2& \cdots & N_L\end{array}\right)}^T$,

а f – столбец функций источника. Матрица M называется матрицей массы, K – матрицей жесткости, а появление матриц V, T обусловлено термическими и вязкими эффектами.

Система (6) является системой интегро-дифференциальных уравнений Вольтерры. Перепишем (6) следующим образом:

$\begin{array}{c}\widehat{M}\frac{d}{dt}y=\widehat{K}y+\widehat{V}{\displaystyle \underset{0}{\overset{t}{\int }}G(t-\tau )y(\tau )d\tau }+\\ +\widehat{T}{\displaystyle \underset{0}{\overset{t}{\int }}G(t-\tau )\frac{d}{d\tau }y(\tau )d\tau }+\widehat{f},\end{array}$      (7)

Где

$\begin{array}{c}\widehat{M}=\left(\begin{array}{cc}M& 0\\ 0& I\end{array}\right),\widehat{K}=\left(\begin{array}{cc}0& -K\\ I& 0\end{array}\right),\\ \widehat{V}=\left(\begin{array}{cc}0& -V\\ 0& 0\end{array}\right),\widehat{T}=\left(\begin{array}{cc}-T& 0\\ 0& 0\end{array}\right),\widehat{f}=\left(\begin{array}{c}f\\ 0\end{array}\right),\end{array}$

G(t – τ) – ядро интегральных слагаемых:

$G(t-\tau )=\frac{1}{\sqrt{\pi (t-\tau )}}$,

а y – новый столбец неизвестных:

$y=\left(\begin{array}{c}\frac{dp}{dt}\\ p\end{array}\right)$.

Интегрирование во времени уравнения (7)

Как было отмечено выше, (7) представляет собой систему интегро-дифференциальных уравнений. Будем отдельно аппроксимировать дифференциальную и интегральную часть уравнений.

Производная столбца dy/dt аппроксимируется центральной разностью:

$\frac{dy}{dt}\approx \frac{y_{n+1}-y_{n-1}}{2\Delta t}$,

где Δt – шаг дискретизации по времени, а нижние индексы обозначают текущий временной шаг y_n = y(nΔt).

Для самого первого шага используется аппроксимация разностью справа:

$\frac{dy}{dt}\approx \frac{y_1-y_0}{\Delta t}$.

Введем следующие квадратурные аппроксимации для интегральных слагаемых в (7) [9]:

$\begin{array}{c}{\displaystyle \underset{0}{\overset{t_n}{\int }}G(t_n-\tau )y(\tau )d\tau }\approx {\displaystyle \sum _{j=0}^n{b}_{nj}{y}_j},\\ {\displaystyle \underset{0}{\overset{t_n}{\int }}G(t_n-\tau )\frac{d}{d\tau }y(\tau )d\tau }\approx \frac{1}{2\Delta t}{\displaystyle \sum _{j=1}^n{b}_{nj}(y_{j+1}-y_{j-1})},\end{array}$ (8)

где b_nj дается формулой

$\begin{array}{c}{b}_{nj}\equiv {\displaystyle \underset{t_j}{\overset{t_{j+1}}{\int }}\frac{d\tau }{\sqrt{\pi }\sqrt{t_{n+1}-\tau }}}=\\ =\frac{2}{\sqrt{\pi }}\left(\sqrt{t_{n+1}-t_j}-\sqrt{t_{n+1}-t_{j+1}}\right).\end{array}$

Заметим, что для того, чтобы оценить значения интегралов в (8) в момент времени t_n+1, необходимо знать неизвестные y_n только в предшествующие моменты времени t₀, …, t_n. Отметим также, что, несмотря на сингулярность ядра G(t – τ) при t = τ, квадратурные коэффициенты b_nj везде имеют конечные значения.

Таким образом, приходим к следующей неявной конечно-разностной схеме:

 $\begin{array}{c}{y}_{n+1}={\left(\widehat{M}-b_{nn}\widehat{T}\right)}^{-1}\left(\widehat{M}{y}_{n-1}+2\Delta t\widehat{K}{y}_n\right.+\\ +2\Delta t\widehat{V}{\displaystyle \sum _{j=0}^n{b}_{nj}{y}_j}+\\ \left.+\widehat{T}\left({\displaystyle \sum _{j=1}^{n-1}{b}_{nj}(y_{j+1}-y_{j-1})}-b_{nn}{y}_{n-1}\right)\right).\end{array}$             (9)

Верификация конечно-элементной формулировки

Рассмотрим задачу распространения монохроматической звуковой волны частоты ω в тонкой полубесконечной трубе с квадратным сечением со стороной a (см. рис. 2).

Рис. 2. Геометрия модели.

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

Аналитическое описание

Ввиду простоты условия (1) будем исследовать стационарную задачу. Введем аксиальную координату x вдоль трубы. В таком случае давление подчиняется уравнению Гельмгольца: $\frac{{\partial }^2\tilde{p}(x,\omega )}{\partial x^2}+{\omega }^2\rho (\omega )\kappa (\omega )\tilde{p}(x,\omega )=0$,

где ρ(ω) – комплексная плотность воздуха с учетом вязких эффектов, а κ(ω) – комплексная сжимаемость воздуха с учетом термических эффектов [10]:

 $\begin{array}{c}\rho (\omega )={\rho }_0\left(1+\frac{s^{2}F(z_{\nu })}{1-F(z_V)}\right),\\ \kappa (\omega )=\frac{1}{{\rho }_0{c}_0^2}\left(\gamma -(\gamma -1)\frac{1-F(z_T)}{s^{2}F(z_T)}\right).\end{array}$       (10)

Здесь

$\begin{array}{c}F(z)=\frac{2J_1(z)}{z{J}_0(z)},\\ z_{V,T}=i^{3/2}\sqrt{2}\frac{R}{{\delta }_{V,T}},\begin{array}{cc}R=s^2\frac{a}{2},& s^2=\frac{8}{7}\end{array},\end{array}$

где J₀(z) и J₁(z) – функции Бесселя соответственно нулевого и первого порядка. Формулы (10) выполняются при условиях R > 10^-5 м и Rf ^3/2 < 10⁴ м с^–3/2 [11]. Первое условие означает, что труба является “широкой” в смысле Вестона [12], т.е. ее радиус больше длины свободного пробега. Второе же условие означает, что труба одновременно является “узкой”, т.е. длина звуковой волны больше длины свободного пробега. Число s² является фактором формы для квадратного сечения [13].

Оценим ошибку формул (10). Для этого опишем, из каких соображений она была получена. Здесь повторяются рассуждения из [13]. В случае трубы с круглым сечением формула (10) получается аналитически (для этого надо взять s² = 1, а в качестве R – радиус сечения). В случае трубы с сечением в форме щели также можно получить аналитическое выражение для комплексных плотности и сжимаемости. Однако если в формулах (10) взять в качестве фактора формы s² = 2/3, а в качестве R – полуширину щели, то на высоких частотах будут получаться правильные асимптотические представления для действительной и мнимой частей комплексной плотности, а на низких – для мнимой части. На средних же частотах ошибка довольно мала. Из этого факта делается вывод, что формулы (10) подходят для описания труб с любым постоянным сечением. В [14] проводится моделирование труб с сечениями различной формы с помощью метода конечных элементов. В частности, получаются асимптотические пределы для действительной и мнимой частей комплексной плотности воздуха. Путем подбора фактора формы s², дающего правильную асимптотику для комплексной плотности на низких частотах, можно хорошо аппроксимировать комплексные плотность и сжимаемость для труб с произвольным постоянным сечением на широком диапазоне частот. В то же время, конечно, формулы (10) все же будут давать небольшую ошибку, особенно на средних частотах.

На основании трубы выполняется следующее граничное условие:

$\frac{d\tilde{p}}{dx}(\mathrm{0,}\omega )=\frac{{\rho }_0}{a^2}\tilde{u}(\omega )$,

где $\tilde{u}(\omega )$ – Фурье-образ функции источника.

Давления и их производные в двух точках трубы $\tilde{p}(x_1,\omega )$, $d\tilde{p}/dx(x_1,\omega )$ и $\tilde{p}(x_2,\omega )$, $d\tilde{p}/dx(x_2,\omega )$ связаны следующим матричным соотношением:

$\begin{array}{c}\left(\begin{array}{c}\tilde{p}(x_2,\omega )\\ \frac{d\tilde{p}}{dx}(x_2,\omega )\end{array}\right)=\\ =\left(\begin{array}{cc}\cos \left(k(\omega )(x_2-x_1)\right)& -\frac{Z(\omega )}{{\rho }_0\omega }\sin \left(k(\omega )(x_2-x_1)\right)\\ \frac{{\rho }_0\omega }{Z(\omega )}\sin \left(k(\omega )(x_2-x_1)\right)& \cos \left(k(\omega )(x_2-x_1)\right)\end{array}\right)\times \\ \times \left(\begin{array}{c}\tilde{p}(x_1,\omega )\\ \frac{d\tilde{p}}{dx}(x_1,\omega )\end{array}\right),\end{array}$

где k(ω) и Z(ω) – соответственно комплексные волновое число и импеданс:

$\begin{array}{cc}k(\omega )=\omega \sqrt{\rho (\omega )\kappa (\omega )},& Z(\omega )=\sqrt{\frac{\rho (\omega )}{\kappa (\omega )}}\end{array}$.

С учетом граничного условия получается следующее выражение для давления в точке x = x₀:

$\tilde{p}(x_0,\omega )=\frac{1}{i\omega a^2}\tilde{u}(\omega )Z(\omega )\exp (-ik(\omega )x_0)$.               (11)

Для получения выражения для давления и производной давления во временной области, применим к (11) обратное преобразование Фурье:

$p(x_0,t)=\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\tilde{p}(x_0,\omega )\exp (i\omega t)d\omega }$,

 $\frac{dp}{dt}(x_0,t)=\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}i\omega \tilde{p}(x_0,\omega )\exp (i\omega t)d\omega }$.              (12)

Из (12) и свойства линейности преобразования Фурье можно получить, что при отсутствии термовязких эффектов (т.е. при k(ω) = ω/c₀, Z(ω) = ρ₀c₀) производная давления по времени dp/dt(x₀, t) будет пропорциональна функции источника f(t), сдвинутой на x₀/c₀.

Численное моделирование

Будем моделировать трубу длины L = 120 см с квадратным сечением со стороной a = 2 мм с помощью конечно-элементной формулировки во временной области (6). При таких размерах трубы условие Rf ^3/2 < 10⁴ м с^–3/2 будет справедливо на частотах до 44 кГц. В основании трубы находятся точечные источники объемного ускорения, в качестве функции источника используется следующая функция:

$f(t)=\cos \left(2\pi f_0(t-t_0)\right)\exp \left(-\frac{{(t-t_0)}^2}{2{\sigma }^2}\right)$,

где c₀t₀ = 0.4 м, c₀σ ≈ 0.0858 м, f₀ = 2000 Гц. Данные параметры обеспечивают необходимую локализацию функции источника как во временной, так и в частотной областях.

В момент времени t = 0 амплитуда огибающей функции источника составляет только примерно 0.002% от максимальной, поэтому с определенной точностью можно дополнить (6) начальным условием y₀ = 0. Следует также обратить внимание, что граничные условия уже учтены в (6) с помощью матриц $\widehat{T}$, $\widehat{V}$, а граничные условия Неймана в конечно-элементной формулировке учитываются автоматически. Следует напомнить, что в данной задаче на боковые стенки накладываются термовязкие граничные условия (5), а на торцы – граничные условия Неймана.

Метод конечных элементов реализуется на тетраэдральных сетках различной плотности. Каждая сетка собирается следующим образом: трубка разбивается на равные параллелепипеды, затем каждый параллелепипед разбивается на пять тетраэдров. Характерные размеры сеток (наименьшие стороны параллелепипедов) составляют 2 мм, 1 мм, 0.666… мм и 0.5 мм, число тетраэдров – соответственно 600, 2400, 5400 и 9600. Функции формы линейные. Плотность среды и скорость звука имеют следующие значения:

ρ₀ = 1.2043 кг/м³, c₀ = 343.2 м/с.

Материальные константы имеют следующие значения:

$\begin{array}{c}\nu =1.51\times {10}^{-5}{\text{ м}}^2/\text{с},{\gamma }_a=\mathrm{1.4,}\\ c_p=1.03\text{\hspace{0.33em}кДж/(кг К)}\text{,}\kappa =0.025\text{\hspace{0.33em}Вт/(м К)}.\end{array}$

Моделирование проводилось до момента времени, соответствующего распространению импульса на 1.5 м.

Поскольку квадратурные коэффициенты b_nj стремятся к нулю при n >> j, для численного счета можно взять не все коэффициенты. Анализ показывает, что для самой грубой сетки достаточно взять последние 4000 квадратурных коэффициентов, что для данной задачи является приблизительно половиной от всех используемых в (9) коэффициентов. Для каждой следующей сетки число коэффициентов подбиралось таким образом, чтобы оно соответствовало тому же промежутку времени, что и для самой грубой сетки.

Анализ сходимости схемы по времени показывает, что, несмотря на тот факт, что использование центральных производных дает ошибку порядка O((Δt)²), использование квадратур [9] приводит к ошибке порядка O(Δt), таким образом полученная схема теоретически имеет первый порядок точности. В то же время, как показывает численный счет, при использовании обычных правых производных решение является неустойчивым даже при Δt = 0.01 h_min/c₀, где h_min – длина наименьшего из ребер в сетке. Поэтому, несмотря на кажущуюся избыточность, было принято решение использовать именно центральные производные для конечно-разностного представления производных по времени. При использовании центральных производных решение становится устойчивым при Δt = 0.1 h_min/c₀. Шаг по времени для каждой сетки подбирается в соответствии с этим условием, таким образом, шаг по времени зависит от шага сетки.

Можно также показать, что использование конечно-элементной постановки с линейными функциями формы дает ошибку по пространству порядка $O(h_{\max }^2)$, где h_max – длина наибольшего из ребер в сетке.

На рис. 3 представлена вычисленная с помощью ранее описанной схемы производная давления по времени dp/dt(x₀, t) на расстоянии 50 см от источника и ее сравнение с формулой (12). В качестве значения производной давления берется среднее значение производной по всем узлам сетки, находящимся на сечении трубы в 50 см от источника. Амплитуда импульса нормирована на максимальное значение при отсутствии термовязких эффектов.

Рис. 3. Сравнение численного моделирования с аналитическим решением. Пунктирной линией обозначен импульс на расстоянии 50 см от излучающей поверхности без учета термовязких эффектов, прерывистой линией – с учетом термовязких эффектов, сплошной линией – аналитическое решение (12).

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

Из представленных графиков можно сделать два вывода: во-первых, для тонких труб термовязкие эффекты будут значительными уже на расстояниях порядка 0.5–1 м; во-вторых, результаты численного моделирования достаточно хорошо согласуются с результатами аналитических расчетов (12).

Была также исследована сходимость полученной схемы по норме L₂([0,T]). Вычислялась относительная ошибка численного решения согласно формуле

$\begin{array}{l}\epsilon =\frac{{||{\left(\frac{dp}{dt}\right)}_{\text{числ}}-{\left(\frac{dp}{dt}\right)}_{\text{аналит}}||}_{L_2}}{{||{\left(\frac{dp}{dt}\right)}_{\text{аналит}}||}_{L_2}}\equiv \\ \equiv \frac{\sqrt{{\displaystyle \underset{0}{\overset{T}{\int }}{\left({\left(\frac{dp}{dt}\right)}_{\text{числ}}-{\left(\frac{dp}{dt}\right)}_{\text{аналит}}\right)}^{2}dt}}}{\sqrt{{\displaystyle \underset{0}{\overset{T}{\int }}{\left({\left(\frac{dp}{dt}\right)}_{\text{аналит}}\right)}^{2}dt}}},\end{array}$

где под аналитическим решением понимается формула (12).

Согласно рис. 4, сходимость метода порядка O((Δt)^3/2), а относительная ошибка для сетки с наименьшим шагом составляет примерно 2.5%. Однако на более плотных сетках порядок сходимости заметно ухудшается. Это можно объяснить тем, что численная ошибка на плотных сетках приближается к ошибке самой теоретической формулы (12), что обсуждалось ранее. Косвенно это подтверждается исследованиями сеток на стационарных задачах, где ошибки решения на сетках разной плотности практически не отличаются.

Рис. 4. Зависимость разницы численного и аналитического решения от шага сетки.

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

Заключение

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

Работа выполнена при поддержке гранта РФФИ 19-29-06048. Работа А.И. Королькова была также поддержана Королевским обществом посредством исследовательской стипендии Дороти Ходжкин д-ра Кисиль.

Приложение А основные понятия теории дробных производных

Преобразование Фурье производной функции n-го (целого) порядка дается следующим выражением:

 $\begin{array}{l}\frac{d^n}{d{t}^n}f(t)=\frac{d^n}{d{t}^n}\left[\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\tilde{f}(\omega )\exp (i\omega t)d\omega }\right]=\\ =\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\left[{(i\omega )}^n\tilde{f}(\omega )\right]\exp (i\omega t)d\omega .}\end{array}$               (13)

Определим свертку двух действительных функций f(t) и g(t) следующим образом:

$(f\times g)(t)\equiv {\displaystyle \underset{-\infty }{\overset{t}{\int }}f(\tau )g(\tau -t)d\tau }$.

Пользуясь теоремой Фубини, можно найти преобразование Фурье свертки (теорема о свертке)

$(f\times g)(\omega )=\tilde{f}(\omega )\tilde{g}(\omega )$.              (14)

Поэтому производную свертки можно выразить следующим образом:

 $\frac{d}{dt}(f\times g)=\frac{df}{dt}\times g=f\times \frac{dg}{dt}$.  (15)

Одностороннее преобразование Фурье степенной функции t^-α дается следующим выражением:

${\displaystyle \underset{0}{\overset{\infty }{\int }}{t}^{-\alpha }\exp (-i\omega t)dt={(i\omega )}^{\alpha -1}\Gamma (1-\alpha )}$.     (16)

Здесь Γ – гамма-функция. Соотношение справедливо для α < 1.

С помощью формулы (13) дробная производная f(t) может быть введена следующим образом:

 $\begin{array}{l}\frac{d^{\beta }}{d{t}^{\beta }}f(t)\equiv {\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{(i\omega )}^{\beta }\tilde{f}(\omega )\exp (i\omega t)dt}\equiv \\ \equiv {\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{(i\omega )}^n\tilde{f}(\omega ){(i\omega )}^{\beta -n}\exp (i\omega t)},\end{array}$                       (17)

где n – ближайшая к β целая степень, такая что n > β. Используем (14) и (15), чтобы получить временное представление дробной производной. Выражение (17) можно представить в виде свертки классической производной f(t) и степенного ядра (16)

$\frac{1}{\Gamma (m-\beta )}\frac{1}{t^{\beta +1-n}}$.

Из (15) следует два временных определения дробной производной, различающиеся порядком дифференцирования и интегрирования:

$\frac{d^{\beta }}{d{t}^{\beta }}f(t)=\frac{1}{\Gamma (n-\beta )}\frac{d^n}{d{t}^n}\left(f\times \frac{1}{t^{\beta +1-n}}\right)$

И

$\frac{d^{\beta }}{d{t}^{\beta }}f(t)=\frac{1}{\Gamma (n-\beta )}\left(\frac{d^n}{d{t}^n}f\right)\times \frac{1}{t^{\beta +1-n}}$.

Первое соответствует дробной производной в смысле Римана-Лиувилля [15]:

${}_{a}D{}_t^{\beta }=\frac{1}{\Gamma (n-\beta )}\left(\frac{d^n}{d{t}^n}\right){\displaystyle \underset{a}{\overset{t}{\int }}\frac{f(\tau )}{{(t-\tau )}^{\beta +1-n}}d\tau }$,

а второе – дробной производной в смысле Капуто:

${}_a{}^{C}D{}_t^{\beta }f(t)=\frac{1}{\Gamma (n-\beta )}{\displaystyle \underset{a}{\overset{t}{\int }}\frac{d^{n}f}{d{\tau }^n}\frac{1}{{(t-\tau )}^{\beta +1-n}}d\tau }$.

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