计算科学 - 半无限平面上平面波散射的计算 - 吾爱随笔录

半无限平面上平面波散射的计算

计算科学 pde Python 数值分析

2021-12-11 11:23:26

我试图编写简单的数学运算来绘制由入射平面波在半无限平板上建立的总场，可以在此处找到。

总结一下：

ϕ_{s} (r, θ) = e^{i k r} / 2 (w [e^{i π / 4} \sqrt{2 k r} \sin (\frac{θ + Θ}{2})] + w [e^{i π / 4} \sqrt{2 k r} \sin (\frac{θ - Θ}{2})]),

$\phi_s(r,\theta ) = e^{ikr}/2 \left( w\left[ e^{i\pi/4}\sqrt{2kr}\sin(\frac{\theta + \Theta}{2}) \right]+ w\left[ e^{i\pi/4}\sqrt{2kr}\sin(\frac{\theta - \Theta}{2}) \right]\right) \enspace,$

s 根据文章的极坐标中的散射势。表示入射平面波的角度，是波数。最后 $\Theta$ $k$

w (z) = e^{- z^{2}} (1 - e r f (z))

$w(z) = e^{-z^2}(1-\mathrm{erf}(z)) \enspace$

如果您运行代码，您会看到我没有得到图片中的预期波场。

我的问题如下：是我在执行过程中犯了错误还是给出的公式不正确？

我的计算：

他们的计算：

我的代码如下：

import numpy as np
from scipy.special import erf
import matplotlib.pyplot as plt


"""Plane wave scattering by a plane wave incident on a seminfinite plane
"""

global T
global k
T = np.pi/4.0 # Incident plane wave angle
k = 10        # Wave number


def phi_scattered(r,t):

    W = lambda z: np.exp(-z**2)*(1-erf(z))
    term1 = W(np.exp(1j*np.pi/4.0)*np.sqrt(2*k*r)*np.sin((t+T)/2.0))
    term2 = W(np.exp(1j*np.pi/4.0)*np.sqrt(2*k*r)*np.sin((t-T)/2.0))
    res =  np.exp(1j*k*r)/2.0*(term1 + term2)

    return res

def phi_incident(r,t):
    return np.exp(1j*k*r*np.cos(t-T))



x1 = np.linspace(-50,50,100)
x2 = np.linspace(-50,50,100)
X1, X2 = np.meshgrid(x1,x2)


R = np.sqrt(X1**2+X2**2)
THETA = np.arctan2(X2,X1)

phi = phi_scattered(R,THETA) + phi_incident(R,THETA)
plt.figure()
plt.imshow(np.real(phi), vmin = np.min(np.real(phi)), vmax = np.max(np.real(phi)))
plt.colorbar()
plt.show()

1个回答

我不能直接回答你的问题，但是：

1）使用contour而不是imshow，否则你的情节是x翻转的。

2) 使用与论文中相同的轴限制（即从-20 到 20）。

3）正如其他人提到的，使用与论文中相同的参数（即k=1）。

4) 以与论文中相同的方式设计情节也使情节比较更容易。

我已经解决了这些问题，代码如下：

import numpy as np
from scipy.special import erf
import matplotlib.pyplot as plt

"""Plane wave scattering by a plane wave incident on a semi-infinite plane"""

T = np.pi / 4.0  # Incident plane wave angle
k = 1.0  # 10  # Wave number

def phi_scattered(r,t):
    W = lambda z: np.exp(-z**2)*(1-erf(z))
    term1 = W(np.exp(1j*np.pi/4.0)*np.sqrt(2*k*r)*np.sin((t+T)/2.0))
    term2 = W(np.exp(1j*np.pi/4.0)*np.sqrt(2*k*r)*np.sin((t-T)/2.0))
    res =  np.exp(1j*k*r)/2.0*(term1 + term2)

    return res

def phi_incident(r,t):
    return np.exp(1j*k*r*np.cos(t - T))

x1 = np.linspace(-20, 20, 100)
x2 = np.linspace(-20, 20, 100)
X1, X2 = np.meshgrid(x1,x2)

R = np.sqrt(X1**2 + X2**2)
THETA = np.arctan2(X2, X1)

phi = phi_scattered(R, THETA) + phi_incident(R, THETA)

plt.figure()
plt.contourf(x1, x2, np.real(phi), cmap=plt.cm.gist_rainbow)
plt.colorbar()
plt.contour(x1, x2, np.real(phi), colors='k', linestyles='solid')
plt.axis('image')
plt.savefig('scattering.png', dpi=300)
plt.show()

结果是：

您可以看到它沿 x 轴翻转。

其它你可能感兴趣的问题

上一篇大型稀疏矩阵的动态/无矩阵 SVD 下一篇C中的克罗内克积和基收缩（即BATranspose [B]）？