Asymptotic study of flows induced by oscillations of cylindrical bodies

Full Text

Исследования аэрогидродинамики колебательного движения тел в жидкости имеют большое значение для широкого спектра прикладных направлений, таких как морское строительство [1, 2], авиационно-космическая инженерия [3, 4], атомно-силовая микроскопия [5, 6], альтернативная энергетика [7], микрофлюидика [8, 9], разработка биомиметических пропульсивныхсистем [10–12] и др.

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

При больших и средних амплитудах колебаний, когда около тел реализуются вихревые отрывные режимы обтекания, изучение аэрогидродинамики главным образом основывается на экспериментах и численных решениях (см., напр., [6, 9, 10, 13]), а также на моделях теории идеальной жидкости [14, 15]. При малых амплитудах колебаний, когда вязкие эффекты проявляются главным образом в пограничном слое около тела, полную картину гидродинамических процессов можно получить, используя асимптотический подход.

Первые результаты в этой области пол учены еще Стоксом [16] в XIX в. Им, в частности, в первом приближении было построено аналитическое решение для случая малоамплитудных колебаний круглого цилиндра. Развитие асимптотических исследований этой задачи проводилось далее в многочисленных работах, напр. [17–19, 28, 29].

Переход к некруглой форме сечения существенно осложняет решение задачи в рамках асимптотического подхода. Исследования в этой области не так многочисленны. Обобщение решения Стокса на цилиндрическое тело с произвольной формой сечения проводилось в работе Така [20], им был предложен метод расчета гидродинамической силы, основанный на численном методе граничных элементов. Развитие этот подход получил в работах [21, 22]. Также следует отметить аналитические результаты для цилиндрических тел с разной формой сечения, полученные в исследованиях [23, 24] для случая высокочастотных малоамплитудных колебаний.

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

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

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

$u_0=U_0\cos \left(\omega t\right)(\cos \alpha ,\sin \alpha ),$

где U₀ – характерная амплитуда скорости колебаний; w – частота колебаний; a – угол между осью колебаний и осью Ox.

Решение гидродинамической задачи будем проводить в подвижной криволинейной ортогональной системе координат $(r,\theta )$ жестко связанной с телом, которая задается конформнным отображением

$z\left(\zeta \right)=x+iy,\zeta =r+i\theta ,\frac{\partial x}{\partial r}=\frac{\partial y}{\partial \theta },\frac{\partial y}{\partial r}=-\frac{\partial x}{\partial \theta }.$

Оно отображает полуполосу $r_0\le r\le \infty \mathrm{,0}\le \theta \le 2\pi$ во внешность сечения цилиндра в физической области. При этом r = r₀ переходит в замкнутую кривую, определяющую профиль сечения цилиндра, бесконечно удаленная точка перейдет в бесконечно удаленную $z(\infty )=\infty$ и $\partial z/\partial \zeta >0.$ Коэффициенты Ламе (масштабные коэффициенты) определяются следующим образом:

$\lambda ={\lambda }_r={\lambda }_{\theta },{\lambda }^2=\frac{\partial z}{\partial \zeta }\frac{\partial \overline{z}}{\partial \overline{\zeta }}=\frac{\partial z}{\partial r}\frac{\partial \overline{z}}{\partial r}=\frac{\partial z}{\partial \theta }\frac{\partial \overline{z}}{\partial \theta }.$

Нормируя пространственные координаты на R, время на w^-1, функцию тока y на U₀R, завихренность w на w, запишем уравнения, описывающие движение вязкой несжимаемой жидкости около цилиндра, в терминах завихренность – функция тока:

$\frac{\partial w}{\partial t}+\frac{k}{{\lambda }^2}\frac{\partial \left(\psi ,w\right)}{\partial \left(r,\theta \right)}={\beta }^{-1}\Delta w,w=\Delta \psi ,$

$\frac{\partial (\psi ,w)}{\partial (r,\theta )}=\frac{\partial \psi }{\partial r}\frac{\partial w}{\partial \theta }-\frac{\partial \psi }{\partial \theta }\frac{\partial w}{\partial r},\Delta =\frac{1}{{\lambda }^2}\left(\frac{{\partial }^2}{\partial r^2}+\frac{{\partial }^2}{\partial {\theta }^2}\right).$ (1)

Здесь безразмерные параметры

$k=U_0/\left(R\omega \right),\beta =R^2\omega /\nu$

определяют безразмерную амплитуду колебаний и безразмерную частоту колебаний (где n – кинематическая вязкость жидкости) цилиндра соответственно, далее будем считать $k\ll 1.$

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

$r=r_0:\psi =\frac{\partial \psi }{\partial r}=\mathrm{0,}$ (2)

$r\to \infty :\psi ~m_{\infty }{e}^r\sin \left(\theta -\alpha \right)\cos t,m_{\infty }={\left(\frac{\partial z}{\partial \zeta }{e}^{-\zeta }\right)}_{\zeta =\infty }.$

2. Асимптотическое разложение

Для решения задачи будем использовать метод формальных асимптотических разложений. Полагая k малой величиной, представим функции y и w в виде ряда по степеням k:

$\psi =k^0{\psi }_0+k^1{\psi }_1+\mathrm{...},w=k^0{w}_0+k^1{w}_1+\mathrm{...}$ (3)

Подставляя (3) в (1) и (2) и собирая члены при одинаковых степенях k, получим серию линейных подзадач.

В главном члене решение описывает нестационарное стоксовское (осцилляционное) течение жидкости около цилиндра:

$\frac{\partial w_0}{\partial t}={\beta }^{-1}\Delta w_0,w_0=\Delta {\psi }_0,$ (4)

$r=r_0:{\psi }_0=\frac{\partial {\psi }_0}{\partial r}=\mathrm{0,}$

$r\to \infty :{\psi }_0=m_{\infty }{e}^r\sin \left(\theta -\alpha \right)\cos t+O\left(1\right).$

В соответствии с имеющимися граничными условиями решение подзадачи в этом члене можно представить в следующем виде:

${\psi }_0=\mathrm{Re}al\left(\left({\psi }_0^1+i{\psi }_0^2\right)e^{it}\right)={\psi }_0^1\cos t-{\psi }_0^2\sin t,$ (5)

$w_0=\mathrm{Re}al\left(\left(w_0^1+i{w}_0^2\right)e^{it}\right)=w_0^1\cos t-w_0^2\sin t.$

Подставляя (5) в (4), перепишем эту подзадачу как

$\left\{\begin{array}{r}\Delta {\psi }_0^{\mathrm{1,2}}=w_0^{\mathrm{1,2}}\\ -w_0^2={\beta }^{-1}\Delta w_0^1\\ w_0^1={\beta }^{-1}\Delta w_0^2\end{array}\right.,$ (6)

$r=r_0:{\psi }_0^{\mathrm{1,2}}=\frac{\partial {\psi }_0^{\mathrm{1,2}}}{\partial r}=\mathrm{0,}$

$r\to \infty :{\psi }_0^1~m_{\infty }{e}^r\sin \left(\theta -\alpha \right),{\psi }_0^2=O\left(1\right).$

Далее ее решение будет определяться численно.

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

$\Delta {w_1}^s=\frac{\beta }{{\lambda }^2}{\frac{\partial \left({\psi }_0,w_0\right)}{\partial (r,\theta )}}^s=\frac{\beta }{\left(2{\lambda }^2\right)}\frac{\partial \left({\psi }_0^1,w_0^1\right)}{\partial (r,\theta )}+\frac{\partial \left({\psi }_0^2,w_0^2\right)}{\partial (r,\theta )},$ (7)

$\Delta {{\psi }_1}^s={w_1}^s,$

$r=r_0:{\psi }_1^s=\frac{\partial {{\psi }_1}^s}{\partial r}=0.$

Для симметричных тел на бесконечности будем использовать следующие условия:

$r\to \infty :{\psi }_1^s=O\left(1\right).$

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

$r\to \infty :{\psi }_1^s=O\left(e^r\right).$

Решение стационарной задачи далее также будет определяться численно.

3. Численное решение

Численное решение подзадач (6), (7) проводится с помощью метода конечных разностей. Решение строится в ограниченной области $r_0\le r\le r_{\infty },0\le \theta \le 2\pi$ на равномерных (в пространстве конформного отображения) сетках с $n_r\times n_{\theta }$ узлов. Дискретные значения переменных определяются в узлах расчетной сетки.

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

Остановимся более подробно на дискретизации граничных условий (ГУ). На границе цилиндра ГУ преобразуются к следующему виду:

${\tilde{\psi }}_0^{\mathrm{1,2}}\left(r_0\right)=\mathrm{0,}{\tilde{w}}_0^{\mathrm{1,2}}\left(r_0\right)=\frac{1}{\lambda {\left(r_0\right)}^2}\frac{2}{h^2}{\tilde{\psi }}_0^{\mathrm{1,2}}\left(r_0+h_r\right),$

где тильда обозначает приближенную сеточную функцию; h_r – шаг сетки по координате r. Условие для завихренности здесь отражает широкоизвестное условие Тома 2-го порядка точности [25].

На внешней границ е расчетной области в первой подзадаче (6) ГУ записываются в следующем виде:

${\tilde{\psi }}_0^1\left(r_{\infty }\right)+\frac{{\partial }^{\left(up2\right)}{\tilde{\psi }}_0^1}{\partial r}\left(r_{\infty }\right)=2m_{\infty }{e}^{r_{\infty }}\sin \left(\theta -\alpha \right),$

${\tilde{\psi }}_0^2\left(r_{\infty }\right)+\frac{{\partial }^{\left(up2\right)}{\tilde{\psi }}_0^2}{\partial r}\left(r_{\infty }\right)=\mathrm{0,}\frac{{\partial }^{\left(up2\right)}{\tilde{w}}_0^{\mathrm{1,2}}}{\partial r}\left(r_{\infty }\right)=0.$

Индекс (up2) у символа частной производной указывает здесь на использование направленной разности второго порядка точности для аппроксимации соответствующей производной. Представленные условия призваны смоделировать правильное осцилляционное поведение течения на удаленной (на конечное расстояние) границе. Дополнительные замечания по критериям выбора соответствующих условий даны в разделе 5.1.

Во второй подзадаче дискретные граничные условия на внешней границе задаются для симметричных тел как

$e^{2r_{\infty }}\frac{{\partial }^{\left(up\right)}{\tilde{\psi }}_1^s}{\partial r}\left(r_{\infty }\right)-e^{2\left(r_{\infty }-h_r\right)}\frac{{\partial }^{\left(up\right)}{\tilde{\psi }}_1^s}{\partial r}\left(r_{\infty }-h_r\right)=\mathrm{0,}$ (8)

$2{\tilde{w}}_1^s\left(r_{\infty }\right)+\frac{{\partial }^{\left(up2\right)}{\tilde{w}}_1^s}{\partial r}\left(r_{\infty }\right)=\mathrm{0,}$

для несимметричных тел как

$e^{-r_{\infty }}\frac{{\partial }^{\left(up\right)}{\tilde{\psi }}_1^s}{\partial r}\left(r_{\infty }\right)-e^{-\left(r_{\infty }-h_r\right)}\frac{{\partial }^{\left(up\right)}{\tilde{\psi }}_1^s}{\partial r}\left(r_{\infty }-h_r\right)=\mathrm{0,}$ (9)

$2{\tilde{w}}_1^s\left(r_{\infty }\right)+\frac{{\partial }^{\left(up2\right)}{\tilde{w}}_1^s}{\partial r}\left(r_{\infty }\right)=0.$

Здесь аналогично индекс (up) указывает на использование направленной разности 1-го порядка точности.

Решение результирующих систем линейных алгебраических уравнений проводилось с помощь геометрического многосеточного метода. Использовался V-цикл, количество пре- и поститераций сглаживания в котором варьировалось от 3 до 6 в зависимости от сетки и геометрии цилиндра. В качестве метода сглаживания использовался метод релаксации. Решение считалось найденным при достижение невязки 10^-8.

Программный код был реализован на языке С⁺⁺ на базе кода, созданного в работе [26]. Благодаря использованию многосеточного метода алгоритм решения имеет высокую скорость сходимости, характерное время полного расчета (в двух приближениях) одного случая на сетке с 2.6 ×10⁵ узлов составляет от 30 с до 5 мин (в зависимости от формы тела) на одном ядре мобильного процессора Intel Core i7-1165G7.

4. Определение гидродинамической силы, действующей на цилиндр

Одной из важнейших характеристик, исследуемых течений, является гидродинамическая сила $f=(f_x,f_y),$ действующая на цилиндр. Ее компоненты (обезразмерненные на ) можно определить по следующим формулам:

$f_x=\frac{1}{2\pi \beta }\left({\int }_0^{2\pi }\omega {\frac{\partial x}{\partial \theta }}_{r=r_0}d\theta +{\int }_0^{2\pi }y{\frac{\partial \omega }{\partial r}}_{r=r_0}d\theta \right),$

$f_y=\frac{1}{2\pi \beta }\left({\int }_0^{2\pi }\omega {\frac{\partial y}{\partial \theta }}_{r=r_0}d\theta -{\int }_0^{2\pi }x{\frac{\partial \omega }{\partial r}}_{r=r_0}d\theta \right).$

В соответствии с используемым разложением сила представляется в виде суммы слагаемых

$\begin{array}{c}{f}_x\approx {f_x}^1\cos t-{f_x}^2\sin t+k{f_x}^s,\\ f_y\approx {f_y}^1\cos t-{f_y}^2\sin t+k{f_y}^s,\end{array}$

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

${f_x}^{\mathrm{1,2,}s}=\frac{1}{2\pi \beta }\left(\begin{array}{c}\sum _i{\tilde{w}}_{\mathrm{0,1}}^{\mathrm{1,2,}s}\left(r_0,i{h}_{\theta }\right)\frac{{\partial }^{\left(c\right)}x}{\partial \theta }\left(r_0,i{h}_{\theta }\right)h_{\theta }+\\ +\sum _{i}y\left(r_0,i{h}_{\theta }\right)\frac{{\partial }^{\left(h\right)}{\tilde{w}}_{\mathrm{0,1}}^{\mathrm{1,2,}s}}{\partial r}\left(r_0,i{h}_{\theta }\right)h_{\theta }\end{array}\right),$

${f_y}^{\mathrm{1,2,}s}=\frac{1}{2\pi \beta }\left(\begin{array}{c}\sum _i{\tilde{w}}_{\mathrm{0,1}}^{\mathrm{1,2,}s}\left(r_0,i{h}_{\theta }\right)\frac{{\partial }^{\left(c\right)}y}{\partial \theta }\left(r_0,i{h}_{\theta }\right)h_{\theta }-\\ -\sum _{i}x\left(r_0,i{h}_{\theta }\right)\frac{{\partial }^{\left(h\right)}{\tilde{w}}_{\mathrm{0,1}}^{\mathrm{1,2,}s}}{\partial r}\left(r_0,i{h}_{\theta }\right)h_{\theta }\end{array}\right).$

Здесь индекс (c) указывает на использование центрально-разностной аппроксимации, а индекс (h) – на использование следующей аппроксимации 2-го порядка точности с компактным шаблоном:

$\frac{{\partial }^{\left(h\right)}\tilde{w}}{\partial r}\left(r_0,i{h}_{\theta }\right)=\frac{\tilde{w}\left(r_0+h\right)-\tilde{w}\left(r_0\right)}{h_r}+$

$+\frac{h_r}{2}\left(-g\left(r_0,i{h}_{\theta }\right)\lambda {\left(r_0\right)}^2+\frac{-\tilde{w}\left(r_0,\left(i-1\right)h_{\theta }\right)+2\tilde{w}\left(r_0,i{h}_{\theta }\right)-\tilde{w}\left(r_0,\left(i+1\right)h_{\theta }\right)}{{h_{\theta }}^2}\right),$

при построении которой используются исходные уравнения для завихренности (6), (7), записанные в виде

$\Delta w=g.$

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

Традиционно в численных и экспериментальных работах для аппроксимации нестационарной силы (f_a), действующий на тело вдоль оси колебаний при установившемся д вижении, используют представление Морисона [27]:

$f_{\alpha }\approx -\frac{4}{3{\pi }^2}k{C}_D\cos t+\frac{1}{2}{C}_M\sin t,$

где C_D называется коэффициентом сопротивления или демпфирования, а C_M называется коэффициентом инерциальных сил. В нашем случае коэффициенты Морисона могут быть вычислены как

$\begin{array}{c}{C}_D=-\frac{3{\pi }^2}{4k}\left(f_x^1\cos \alpha +f_y^1\sin \alpha \right),\\ C_M=-2\left(f_x^2\cos \alpha +f_y^2\sin \alpha \right).\end{array}$

Заметим также, что вычисленная в используемой подвижной системе координат сила f отлична от силы f ^* в неподвижной системе координат на величину, называемую силой Крылова–Фруда $f^{FK}$ $(f=f^{*}+f^{FK}).$ Она равна произведению ускорения на площадь S поперечного сечения тела:

$f^{FK}=(\cos \alpha ,\sin \alpha )\sin tS.$

5. Результаты

5.1. Колебания круглого цилиндра

В качестве первого примера рассмотрим колебания круглого цилиндра в жидкости. Преобразование координат для круглого цилиндра задается как

$x=e^r\cos \theta ,y=e^r\sin \theta ;$

${\lambda }^2=e^{2r},r_0=0.$

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

 

Рис. 1. Характерная структура сетки (n = 100×100) в окрестности круглого цилиндра (в физической плоскости).

 

На базе этого примера проведем непосредственную апробацию модели, сравнивая результаты с известным аналитическим решением. В главном члене оно было полученное впервые в знаменитой работе Стокса [16]. Его можно записать в следующем виде [28]:

$\begin{array}{c}{\psi }_0=\mathrm{Re}al\left(({\psi }_0^1+i{\psi }_0^2)e^{it}\right)\sin \theta ,\\ {\psi }_0^1+i{\psi }_0^2=e^r-e^{-r}-S_x{e}^{-r}-z^{-2}{W}_1,\end{array}$

$z=\sqrt{i\beta },W_1=-2z\frac{K_1\left(z{e}^r\right)}{K_0\left(z\right)},S_x\left(z\right)=\frac{2K_1\left(z\right)}{z{K}_0\left(z\right)},$

здесь K_j — модифицированные функции Бесселя второго рода.

Стационарную часть решения в следующем приближении также можно выразить аналитически [19], [28]:

${\psi }_1^s=\sin 2\theta \left(\Psi \left(r\right)-\Psi \left(0\right)e^{-2r}-\left({\Psi }^{\text{'}}\left(0\right)+2\Psi \left(0\right)\right)\left(1-e^{-2r}\right)/2\right),$

$\begin{array}{c}\Psi \left(r\right)=-\frac{e^{4r}}{48}\frac{\beta \pi }{2}{\int }_1^{\infty }G\left(e^{r+x}\right){\left(e^{2x}-1\right)}^{3}dx,\\ G\left(r\right)=\frac{\beta \pi }{4}\mathrm{Im}ag\left(\frac{\overline{S}{\Phi }_0}{r^2}+2S{\Phi }_2\overline{{\Phi }_0}-S{\Phi }_2\right),\end{array}$

${\Phi }_k\left(r\right)=\frac{K_k\left(zr\right)}{K_0\left(z\right)},S=\frac{K_2\left(z\right)}{K_0\left(z\right)}.$

Структура аналитических решений учитывалась при выборе ГУ, представленных в разд. 3, на внешней границе расчетной области.

На рис. 2 представлена относительная погрешность $er{r}_i^j=({\psi }_i^j-{\tilde{\psi }}_i^j)/{\psi }_i^j$ численного решения, полученного на разных сетках при b = 314 в нестационарном и стационарном слагаемых. Как можно видеть, наибольшее значение относительной погрешности в главном члене достигается в пограничном слое около тела. На самой грубой сетке с $n=256\times 256$ узлами максимальная погрешность достигает приблизительно 0.82%, на наилучшей сетке с $n=1024\times 1024$ узлами – чуть менее 0.01%.

 

Рис. 2. Относительная погрешность численного решения ~y01 (а) и ~y0s (б) при b = 314, q = p/4, вычисленного на сетках с n = nr × nq узлов с r∞ = 3: 1 – n = 512 × 512, 2 – n = 1024 × 1024, 3 – n = 2566 × 256.

 

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

Как следствие, в этой области при расчете относительной погрешности возникает деление на ноль, что и объясняет ее резкий скачек. Другое локальное значение максимума погрешности находится в стационарном погранслое у поверхности цилиндра, на лучшей сетке оно не превышает 2%. Во всей остальной области относительная погрешность численного решения в среднем не превышает 0.6% на самой грубой сетке и 0.04% на самой точной.

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

$\begin{array}{c}{C}_D=-\frac{3{\pi }^2}{4k}\mathrm{Im}ag\left(S_x\right)+O\left(k\right),\\ C_M=2+2\mathrm{Re}al\left(S_x\right)+O\left(k^2\right).\end{array}$ (10)

В табл. 1 представлены значения коэффициентов разложения Морисона вычисленные на сетках с количеством узлов $n_1=256\times 256$, $n_2=512\times 512$, $n_3=1024\times 1024$ для разных значений b. Там же представлены значения $C_D,C_M$ в приближении $\beta \to \infty$ [17] (их можно получить из (10), раскладывая S_x в ряд по степеням b):

$\beta \to \infty :\begin{array}{c}{C}_D=\frac{3{\pi }^2}{4k}\left(\sqrt{2}{\beta }^{-1/2}+{\beta }^{-1}\right)+O\left({\beta }^{-3/2}\right)+O\left(k\right),\\ C_M=2+2\sqrt{2}{\beta }^{-1/2}+O\left({\beta }^{-3/2}\right)+O\left(k^2\right);\end{array}$ (11)

и данные прямого численного моделирования [24].

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

Отметим хорошую согласованность рассматриваемой асимптотики (полученной в предположении k << 1) с данными прямого численного моделирования: при относительно малых амплитудах (k = 0.16, 0.22) разница между оценками C_D не превышает 6%.

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

5.2. Эллиптический цилиндр

В качестве второго примера рассмотрим колебания эллиптического цилиндра в жидкости. Преобразование координат для этого случая задается следующим образом:

$\begin{array}{c}x=\sqrt{1-b^2}\cosh r\cos \theta ,\\ y=\sqrt{1-b^2}\sinh r\sin \theta ,b=\frac{R_y}{R_x}<\mathrm{1,}\end{array}$

$r_0=\ln \left(\frac{\sqrt{1+b}}{\sqrt{1-b}}\right),\lambda =\left(\sqrt{1-b^2}\right)\sqrt{{\cosh }^{2}r-{\cos }^2\theta },$

где R_x, R_y – длины полуосей сечения цилиндра, расположенных вдоль соответствующих осей Ox, Oy в физической плоскости. Размер большей полуоси R_x в дальнейшем будем принимать в качестве характерного размера профиля (R = R_x).

Так как при малых значениях b геометрия тела существенно отличается от рассмотренного случая круглого цилиндра (рис. 3), проведем дополнительный анализ сеточной сходимости для b = 0.1.

 

Рис. 3. Характерная структура сетки (n = 100 × 100) в физической плоскости в окрестности эллиптического цилиндра при b = 0.1.

 

В табл. 2 представлены результаты вычисления коэффициентов силы в расчетных областях разного размера (для разных значений r_∞) и на сетках с разным количеством узлов.

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

Проведем анализ зависимости коэффициентов C_M и C_D гидродинамической силы от отношения полуосей эллипса, переходя от круглого цилиндра при b = 1 к бесконечно тонкой пластине при b = 0. Графики соответствующих зависимостей, полученные в настоящей работе, представлены на рис. 4 для b = 772, a = p/2 одновременно с результами прямого численного моделирования и графиками асимптотических зависимостей, выведенных для случая b → ∞ в [24]:

$C_D^{\left(3\right)}=\frac{3{\pi }^2}{8k}{\left(1+b\right)}^2\left\{{\beta }^{-1/2}\frac{1}{\sqrt{2}}\left(A_1{\cos }^2\alpha +B_1{\sin }^2\alpha \right)+{\beta }^{-1}\left[{\cos }^2\alpha \left(-\frac{1}{2}+\sum _{j-1}^{\infty }j{A}_j^2\right)+{\sin }^2\alpha \left(-\frac{1}{2b^2}+\sum _{j-1}^{\infty }j{B}_j^2\right)\right]\right\},$

$C_M^{\left(3\right)}=\left(1+b\right)\left(b{\cos }^2\alpha +{\sin }^2\alpha \right)+{\beta }^{-1/2}\frac{{\left(1+b\right)}^2}{\sqrt{2}}\left(A_1{\cos }^2\alpha +B_1{\sin }^2\alpha \right),$

$\begin{array}{c}{A}_j=\frac{1}{\pi }{\int }_{-\pi }^{\pi }\frac{\sin \theta \sin \left(j\theta \right)}{\lambda \sqrt{{\displaystyle {\cosh }^2{r}_0-{\cos }^2\theta }}}d\theta ,\\ B_j=\frac{1}{\pi }{\int }_{-\pi }^{\pi }\frac{\cos \theta \cos \left(j\theta \right)}{\lambda \sqrt{{\displaystyle {\cosh }^2{r}_0-{\cos }^2\theta }}}d\theta .\end{array}$ (12)

 

Рис. 4. Изменение коэффициентов CM (а) и CD (б, в) в зависимости от соотношения полуосей b для эллиптического цилиндра при b = 314, k = 0.032, k = 0.064. Сравнение малоамптитудной асимптотики (результаты настоящей работы) с высокочастотным приближением (12) и результатами прямого численного моделирования: 1 – численное моделирование при k = 0.064; 2 – асимптотика b → ∞, k → 0; 3 – численное моделирование при k = 0.032; 4 – асимптотика k → 0.

 

Здесь индекс (03) указывает на удержание трех членов при разложении по степеням b.

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

Значения C_D при b > 0.3, полученные в разных приближениях (в том числе по данным численного моделирования), также хорошо совпадают. В области b → 0 слагаемые ряда (12) для C_D попеременно стремятся то к плюс бесконечности, то к минус бесконечности, что ставит вопрос о границах применимости формулы при малых b. Как можно видеть, высокочастотная асимптотика (12) для рассматриваемого b хорошо совпадает с асимптотическим решением, полученным в нашей работе, начиная с b ≈ 0.2.

В табл. 3 приведено сравнение C_D полученного в ходе численного асимптотического решения для произвольных b с высокочастной асимптотикой (12).

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

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

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

На рис. 5, 6 представлены характерные структуры вторичных течений при разных значениях b для b = 4000 и b = 134 При b = 1 для всех рассматриваемых частот можно наблюдать один и тот же режим течения, в литературе он иногда называется закритическим [28].

 

Рис. 5. Картины вторичных стационарных течений. Линии тока y1s = const при b = 4000, a = p/2: (а–е) – b = 1, 0.5, 0.3, 0.25, 0.2, 0.1.

 

Подобный режим течения был впервые теоретически описан Шлихтингом в 1932 г. [29]. Он реализуется при всех b > 37.064 [28] и характеризуется формированием двух циркуляционных зон около цилиндра. Внутренняя циркуляционная зона состоит из четырех симметрично расположенных ячеек, внешняя граница которых образует окружность. Ее радиус зависит от параметра b. Внешнее циркуляционное течение на удалении от тела представляет собой квадрауполь

$b=\mathrm{1,}r\to \infty :{\psi }_1^s~\left(1-\frac{1}{e^{2r}}\right)\sin 2\left(\theta -\alpha \right).$

При b = 4000 аналогичный режим можно наблюдать при b = 0.5 и b = 0.3. С уменьшением b центры внутренних циркуляционных ячеек смещаются вдоль поверхности цилиндра ближе к оси, перпендикулярной оси колебаний. При b = 0.25 (b = 4000) можно наблюдать значительное увеличение размера ячеек внутренней зоны. При b = 0.2 (b = 4000) они становятся бесконечно большого размера, ядро каждой из ячеек разделяется на два: один из центров остается около тела, второй устремляется к бесконечности. Внешняя циркуляционная зона при этом прижимается к оси колебания, превращаясь в циркуляционные ячейки конечных размеров. При b = 0.1 (b = 4000) эти ячейки полностью исчезают, также исчезают и ядра ячеек бывшей внутренней циркуляционной зоны, расположенные около цилиндра.

Изменение структуры течения при b = 314 происходит похожим образом, только изменяются границы перехода между разными режимами. Трансформации закритического режима начинаются уже при b = 0.5. Также разделение ячеек внешней циркуляционной зоны происходит здесь уже после полного исчезновения внешней циркуляционной зоны. До разделения, при b = 0.4 течение приобретает структуру похожую на докритический режим, наблюдаемый у круглого цилиндра при

5.3. Профиль Жуковского

В качестве заключительного примера рассмотрим колебания цилиндрического тела с сечением в форме руля Жуковского. Преобразование координат для описания геометрии этого тела задается как

$x=m_{\infty }\left(\left(e^r\cos \theta -m\right)+\frac{{(1-m)}^2\left(e^r\cos \theta -m\right)}{e^{2r}-2m{e}^r\cos \theta +m^2}\right);$

$y=m_{\infty }\left(\left(e^r\sin \theta \right)-\frac{{(1-m)}^2\left(e^r\sin \theta \right)}{e^{2r}-2m{e}^r\cos \theta +m^2}\right);$

$m_{\infty }\frac{1+m}{2}m\le \mathrm{1,}$

${\lambda }^2=e^{2r}\frac{\left(e^{2r}-2(2m-1)e^r\cos \theta +{(2m-1)}^2\right)\left(e^{2r}-2e^r\cos \theta +1\right)}{{\left(e^{2r}-2m{e}^r\cos \theta +m^2\right)}^2}{m}_{\infty }^2.$

Здесь x, y, l² записаны в виде вещественных функций от вещественных переменных в той форме, в которой они использовались в программной реализации. Для получения классического профиля с острой кромкой нужно выбрать r₀ = 0, при этом сама кромка становится особой точкой в которой, в частности, l² обращается в ноль. Для проведения расчетов профиль скруглялся, и значение r₀ выбиралось равным 0.02. Форма профиля и структура расчетной сетки в физической области для этого случая представлены на рис. 7.

 

Рис. 6. Картины вторичных стационарных течений. Линии тока y1s = const при b = 314, a = p/2: (а–е) – b = 1, 0.6, 0.5, 0.4, 0.28, 0.1.

 

Рис. 7. Характерная структура сетки (n = 100 × 100) в окрестности симметричного профиля Жуковского при m = 0.5, r0 = 0.02.

 

Далее рассмотрим колебания этого цилиндрического тела вдоль оси симметрии (при a = 0). Проверка на сеточную сходимость и влияние размеров расчетной области проводились практически в каждом расчетном случае, в качестве основной была выбрана расчетная сетка с n = 512 × 512 узлами.

Результаты расчета силовых коэффицентов Морисона по решению в главном члене представлены в табл. 4. Для сравнения с круглым цилиндром (m = 1), значения коэффициентов (и значения b) рассчитаны для случая, когда максимальная толщина крылового профиля равна двум. Как можно видеть значения коэффициента инерциальных сил последовательно снижается при уменьшении m, значение коэффициента сопротивления наоборот возрастает как следствие увеличение площади поверхности тела.

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

Второй вариант в рассматриваемом приближении приведет к появлению логарифмической особенности на бесконечности при решении уравнения Стокса [12] (в используемых координатах речь идет о члене re^r в функции тока), поэтому далее будем искать решение, при котором тело движется с крейсерской скоростью u^st при нулевой средней силе.

Для построения такого решения вместо условия (8), применявшегося ранее для симметричных тел, будем использовать на удаленной границе условие (9). Оно не ограничивает появление направленного потока, но запрещает возрастание скорости на бесконечности. При этом u^st находится из решения задачи.

Характерные структуры вторичных течений, возникающих при колебаниях профиля Жуковского, представлены на рис. 8. Как можно видеть, большинство из них похожи на те, что формируются у круглого и эллиптического цилиндров. При малых числах b (см. b = 20) около колеблющегося формируется одна циркуляционная зона, состоящая из четырех вихрей (как в докритическом режиме у круглого цилиндра), которые из-за формы тела и присутствия направленного потока расположены несимметрично относительно оси Oy. При больших значениях b (см. b = 500, 7000, 4 × 10⁵) хорошо различимы уже две циркуляционные зоны: внутренняя, прижатая к телу, и внешняя, как в закретическом режиме. Случаи b = 50, 100 можно отнести к переходным режимам.

 

Рис. 8. Вторичные стационарные течения около кылового профиля. Линии тока для m = 0.3 (a) b = 20, ust = 0.0069k, (б) b = 50, ust = 0.0068k, (в) $\mathit{\beta }$ = 100, u^st = 0.0065k, (г) $\mathit{\beta }$ = 300, u^st = 0.003k, (д) $\mathit{\beta }$ = 7000, u^st = –0.025k, (е) $\mathit{\beta }$ = 4 ×105, u^st = 0.23k; $\mathit{\beta }$ вычислено по длине хорды.

 

Направление движения профиля с ростом b меняется несколько раз (значения скорости при разных b представлены в подписи к рис. 8). При b = 20, 50, 100 профиль движется острым торцом вперед. Затем, с появлением полноценной внешней циркуляционной зоны (см. b = 500, 7000), направление движения меняется на противоположное. В пределе больших b профиль вновь движется острым концом вперед.

В то же время следует отметить, что скорость стационарного движения во всех режимах является весьма незначительной. С учетом условия $k\ll 1$ она по крайней мере на два порядка ниже амплитуды скорости колебаний.

Заключение

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

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

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

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

Финансирование

Работа выполнена при поддержке гранта Российского научного фонда 22-79-10033 (разд. 1–4, 5.3) и гранта Президента Российской Федерации для государственной поддержки молодых российских ученых – кандидатов наук МК-3245.2022.1.1 (разд. 5.1, 5.2).

About the authors

A. N. Nuriev

Kazan Federal University

Author for correspondence.
Email: Artem.Nuriev@kpfu.ru
Russian Federation, Kazan

O. N. Zaitseva

Kazan Federal University

Email: olga_fdpi@mail.ru
Russian Federation, Kazan

A. M. Kamalutdinov

Kazan Federal University

Email: amkamalutdinov@kpfu.ru
Russian Federation, Kazan

E. E. Bogdanovich

Kazan Federal University

Email: helenbogdanovich03@gmail.com
Russian Federation, Kazan

A. R. Baimuratova

Kazan Federal University

Email: angelina.baimuratova@yandex.ru
Russian Federation, Kazan

References

Supplementary files