О вычислении амплитуды рассеяния включения в сложном поле

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

ПОСТАНОВКА ПРОБЛЕМЫ

Пусть идеальное жидкое пространство, в котором находится включение, полагается однородным с плотностью р и скоростью c . В пространстве подразумевается наличие суммы падающей u n ( x ) и рассеянной от включения u s ( x ) гармонических волн с фактором e ^-,mt, который далее опускается:

u ( x ) = u ( x ) + u s ( x ).                (1)

Включение расположено в области D, внутри которой находится начало координат. Рассеянное поле в волновой зоне, как известно, описывается выражением eikx          - us(x) = u„(x)-rr + o(|xr )’        (2)

x где k = to I c — волновое число; x = (x, y, z) = = (R, 0, ф) — точка наблюдения соответственно в декартовых и сферических координатах;

x x = -- = (sinOcosф,sinOsinф,cosO) — вектор на x единичной сфере, направленный в точку наблюдения; u„ (x) — амплитуда рассеяния включения, характеризующая угловое распределение поля рассеяния включения при его облучении полем u. (x). Обычно выделяется случай плоской падающей волны u (x) = e'kd'x.      Здесь inc d = (sin a cos в, sin a sin в, cos a) — единичный вектор, совпадающий по направлению с волновым вектором плоской волны. В этом случае фиксируется зависимость рассеянной волны us (x) = = us (x, d) от направления распространения падающей плоской волны единичной амплитуды и нулевой фазы и, следовательно, аналогичная зависимость ар u„ (x) = u„ (x, d). Очевидно, что вектор x в функции u„ (x, d) привязан к точке расположения включения. В силу линейности задачи очевидно, что в случае, когда падающая волна представляет собой совокупность плоских волн со спектром g(d)

u_inc ( x ) = ∫ g ( d ) e^{ik d ⋅ x} d s ( d ),             (3)

Ω ^'

ар вычисляется следующим образом (см., например, [2]):

u _∞ ( x ˆ) = ∫ g ( d ) u _∞ ( x ˆ, d )d s ( d ).         (4)

Ω ^'

Здесь Ω ^' — область определения спектра падающей волны; d s ( d ) = sin α d α d β — элемент поверхности единичной сферы.

Известно [3], что если функция u_inc ( x ) является решением однородного уравнения Гельмгольца во всем пространстве x ∈ R ³ (так называемое целое решение), то в ее представлении по плоским волнам (3) участвуют только однородные плоские волны и интегрирование в (3), (4) осуществляется по поверхности единичной сферы Ω ^' =Ω (область видимости). Если же u_inc ( x ) является решением однородного уравнения Гельмгольца только в бесконечной подобласти R ³ и удовлетворяет условиям излучения Зоммерфельда, то интегрирование в (3), (4) происходит как в области видимости (однородные плоские волны), так и вне ее (неоднородные плоские волны).

РЕШЕНИЕ ЗАДАЧИ

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

При построении алгоритма будем исходить из самого сложного случая: первичное падающее поле представляет собой излученное решение.

Воспользуемся выражением для излученного решения в виде разложения по плоским волнам, приведенного в [3, 4, 5], для первичного падающего поля u_inc ( x )

∞∞ uinc(x) =

2 π

-∞ -∞

u ˆ( k , k )

x ^, y e ik d ⋅ x _{d k d k = k}       xy

π

- i ∞

2 π    2

= ^ik   d β    u ˆ( d ) eⁱ

2 π ₀₀

ikR (sin α sin θ cos( β - ϕ ) + cos α cos θ )

× sin α d α .

Здесь uˆ(kx ,ky ) — спектр плоских волн в первичном поле. В (5) подразумевается, что источник поля сосредоточен в области D , включающей в себя начало координат. Пусть включение находится в точке x=x0=(x0,y0,z0)=(R0,θ0,ϕ0). Тогда, согласно (5), включение облучают плоские волны с амплитудой

g ( d ) = ^ik u ˆ( d ) e^{ik d ⋅ x 0} .                 (6)

2 π

С учетом этого запишем (4) в виде интеграла по сферическим углам типа (5):

u _∞ ( θ , ϕ , R ₀, θ ₀, ϕ ₀) =

π

- i ∞

2 π 2

=          u ˆ( α , β ) u _∞ ( θ , ϕ , α , β ) ×

2 π ₀₀

_× e ikR ₀ (sin α sin θ ₀ cos( β - ϕ ₀) + cos α cos θ ₀) s_in α _d α _{d β .}

Между выражением (7), где включение смещено в точку x = x ₀ , и выражениями (2), (4), где включение находится в начале координат, противоречий нет, т. к. углы θ и ϕ в результирующей амплитуде рассеяния u _∞ привязаны, по определению функции u _∞ ( x ˆ, d ) и согласно (4), и согласно (7), к системе координат с началом в точке расположения включения.

Функция u _∞ ( θ , ϕ , R ₀, θ ₀, ϕ ₀) в (7) зависит от R ₀ , θ ₀ , ϕ ₀ как от параметров, характеризующих координаты включения. Как видно из (7), функция u _∞ ( θ , ϕ , R ₀, θ ₀, ϕ ₀) отвечает всем признакам излученного решения вида (5) по переменным R ₀ , θ ₀ и ϕ ₀ (излученное решение со спектром u ˆ( α , β ) u _∞ ( θ , ϕ , α , β )). Это означает, что при фиксированных значениях θ и ϕ функция u _∞ ( θ , ϕ , R ₀, θ ₀, ϕ ₀) удовлетворяет однородному уравнению Гельмгольца везде, где справедливо представление (5) для первичного поля по плоским волнам. А это в свою очередь серьезно облегчает вычисление интеграла (7). А именно вместо того, чтобы производить прямые вычисления или пользоваться общим, но трудоемким методом перевала, можно воспользоваться изящной теоремой, полученной Аткинсоном (Atkinson F.V.) и Уилкоксом (Wilcox C.H.) [6, с. 84], которая в данном случае записывается так:

u ∞ ( θ , ϕ , R 0 , θ 0 , ϕ 0 ) = ^{eikR 0} ∑ ^{∞ Fn ( θ 0, ϕ} _n ⁰⁾ ,        (8)

R _{0 n = 0} ( kR ₀)

где функции F_n ( θ ₀, ϕ ₀) определяются из следующего рекуррентного соотношения:

F n =

^ 1 d d

---(sin 0_O ---) +

( sin 0 0 50 0      0 5 0 0

+ -k A + n (n — 1) I^, sin² 0 0 дф ₀            J 2 in

a F о ( O о , Ф о ) равно

F M^) = U ( O o , ф ,) и „ ( О , ф , О о , ф о).      (10)

При и „ ( 0 , ф , О _о, ф ₀) ^ 1 равенство (8) сходится абсолютно и равномерно при R ₀ >  R ', его можно дифференцировать по R ₀, 0 _о, ф ₀ почленно любое число раз и полученные после этого ряды также сходятся абсолютно и равномерно [6]. Здесь R ' — радиус наименьшего шара, включающего в себя область D источника первичной волны. При функции и „ ( 0 , ф , 0 ₀, ф ₀), отличной от константы, этот вопрос требует отдельного рассмотрения, и наверняка можно сказать только то, что асимптотический ряд (8) по обратным степеням волнового параметра kR ₀ определяется рекурсией (9), (10) единственным образом в силу единственности представления асимптотическими рядами и удовлетворения функции и „ ( 0 , ф , R₀, 0 ₀, ф ₀) однородному уравнению Гельмгольца в указанной области. А следовательно, и в этом случае ряд (8) позволяет вычислить совокупную амплитуду рассеяния и _м ( 0 , ф , R ₀, 0 ₀, ф ₀) с заданной точностью.

В случае, когда первичное поле является целым решением, вычисление интеграла (4), собственного в этом случае, является исходно существенно более простой задачей и может быть проведено точно или приближенно. Остановимся на асимптотике интеграла (4) в случае, когда падающие плоские волны являются только уходящими, т. е. имеют угол падения а е [0, п / 2]. Для этого случая интеграл (7) преобразуется к виду и „ (0,ф, R 0,00, ф0) =

п ik 2п

= — П и ( ав ) и „ ( 0 , ф , а , в) X

^{2 п} 0 0

X e ik 0 (sin а sin 0 0 cos( в - ф 0 ) + cos а cos 0 0 ) sin а d a d p .

Гельмгольца, для совокупной амплитуды рассеяния и „ ( 0 , ф , R ₀, 0 ₀, ф ₀) остается справедливым

Л представление (8)-(10). При 00 = 0 или 00 = — асимптотика (12) для (11) нарушается, и пользоваться представлением (8)–(10) неправомерно.

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

Рассмотрим применение (8)–(10) на некоторых примерах.

Пример 1. Случай сферической падающей волны. Рассмотрим пример расчета радиационного давления на сферическую частицу, находящуюся в сферической волне. Как известно (см., например, [7]), амплитуда рассеяния сферической частицы при падении плоской волны в длинноволновом приближении равна сумме монопольного и дипольного моментов и„ (x, d) = A + A] cos у =

= A 0 + A 1 (sin 0 sin а cos( ф - 3 ) + cos 0 cos а ).   (13)

Здесь у — угол между векторами x и d ; A 0 и A ₁ — постоянные коэффициенты, зависящие от частоты и свойств включения и среды. Пусть источник сферической волны расположен в начале координат. Центр рассеивателя находится в точке x ₀ = (0,0, z ₀), где z ₀ >  0. Для определения и „ ( x ) воспользуемся выражениями (3), (4). Представление сферической волны при z >  0 по плоским волнам в этом случае имеет вид (см., например, [3, 4]):

eikx и (x) = ТТ = inc

п

м

• 7 2 п 2

_ ik Г Г eik(sin а cos px+sin а sin py+z cos а) sin ada dp.(14)

п 00

В точке х ₀ = (0,0, z ₀) поле равно

п

Согласно [4], при 0 < 00 < — интеграл (11) име- ет асимптотику eik x0

U‘c (0,0, z0) = ^J = ikR 0

и „ ( 0 , ф , R 0 , 0 0 , ф 0 ) ~ F 0 ф 0) — ,       (12)

R

п

1 м

2 п 2

_ ik  Г   Г   e kz 0 cos а

2/7" ^   ^

^^п 0  0

e ikz 0 sin a d a d p _--- .

z ₀

где F ₀( 0 ₀, ф ₀) определяется из (10). А в силу того, что (11) удовлетворяет однородному уравнению

Сравнивая последнее выражение с (3), находим

g(d) = —e‘z°cosa. С учетом этого, подставляя (13) 2r в (4), находим

k ²

4 r²

2 п I

- i го

J J

—k a sin a e ²

cos ax

u „ ( x ) =

П

~ — i го

2 r -

= ib I J ( A + r °   °

+A1 (sin 0 sin a cos(ф - в) + cos 0 cos a)) x x eikz°cosa sin a da de =

z >  °.

x e

ik ( x sin a cos в + y sin a sin в + z cos a )

sin a d a d e ,

Пусть включение (13) находится в точке

x ° = ( R ° , 0 ° , ф ° ). Тогда имеем:

u го ( 0 , ф , R ° , 0 ° , ф °) =

_{A 0} e ikz 0 z ₀

^ e^fc °            ^

+ A 1 —- ( kz ° + i ) I cos 0 .

I ^kz °           )

п

1 го

k2f

= —2 I  I ( A o + A j (sin 0 sin a cos( ф - в ) +

^{4 r} о 0

Вычислим u _го ( 0 ) по соотношениям (8)-(1°).

Имеем: R ° = z ° , 0 ° = °, u ( 0 ° , ф ° ) = 1,

.   —-k ² a ² sin ² a

+ cos0cosa))e 2 x x cos a elkkt°(sin0°sinacos(e-ф°)+cosacos0°) sin a da de.    (18)

FM, ф б ) = u ^{( 0} ° , ф б ) u го ^{( 0 ,} ф ^{, 0} °^, ф °⁾ =

= A ° + A , cos 0 cos 0 ° , F ( 0 ° , Ф ° ) = iA 1 cos 0 cos 0 ° ; F n ( 0 ° , ф ° ) ^ °, n = 2,3,...

Последний интеграл, как видно, уже совсем не прост, поэтому прибегнем к описываемому алгоритму. В силу цилиндрической симметрии падающего поля ограничимся случаем ф ° = °. Тогда, согласно (7) и (10), имеем

Подставляя (16), (17) в (8) и полагая 0 ° = ° , получаем в точности выражение (15).

Пример 2. Случай гауссова пучка. В разобранном выше примере сферической волны тривиально вычислялся интеграл (4) и преимущества предлагаемого алгоритма были не столь очевидны. В рассматриваемом ниже случае гауссова пучка преимущества окажутся существеннее.

Под гауссовым пучком здесь понимается плоская волна e ' ^kz , при z = ° пропускаемая через маску с функцией прозрачности ^ ( x , y ), совпадающей с гауссовой кривой

F ° ( 0 ° , ф ° ) =

k

=--- ( A ° + A j (sin 0 ° sin 0 cos ф + cos 0 ° cos 0 ) ) x

2 r i

—1 k 2 a 2 sin2 0° x e 2        cos0°.                              (19)

Из (9) получаем с точностью до третьей степени множителей cos 0 ° и sin 0 °

^ ( x , У ) =

2п^х^у

x ²

к — k²a ² sin ² 0

F 1 ( 0 ° , ф ° ) - ² x 4 r

x (2cos 0 ° A ° + A j (4cos ² 0 ° - 2sin ² 0 ° )cos 0 +

+ A j (5cos 0 ° sin 0 ° - ctg 0 ° cos ² 0 ° )sin 0 cos ф ).   (2°)

Примем далее, что a x = a y = a . Тогда нетрудно показать, что при z >  ° падающее поле определяется выражением (см., например, [8])

го го

U nc ^{( x} , У , ^z ) = J J 'Н ^kx , ^ky ) е ⁽ k^xx ^{+ 1 у У +} )^{d k} x ^d ky , ^z >  °.

-го -го

Подставляя (19), (20) в (8), получаем с точностью до O ( k ^{- 2}R ° ^{- 3} ) значение интеграла (18).

Отметим, что 0 ° = ° является особой точкой дифференциального оператора в (9), поэтому при вычислении F_n ( 0 ° = °), n = 2,3,... возможны особенности, которых можно избежать поворотом системы координат.

Здесь V ( k_x , k_y ) — фурье-образ функции ^ ( x , y ). Переходя к сферической системе координат, имеем:

U nc ^{( x} , y , ^z ) =

ВЫВОДЫ

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

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