Формулы для двумерного моделирования методом FDTD с использованием граничных условий на основе PML

Иван Мелешко
Моделирование-4
В работе F. Collino, C. Tsogka. Application of the PML absorbing layer model to the linear elastodynamic очень хорошо и подробно расписано использование поглощающих слоев на основе PML (Perfect Matched Layers) с разделенными координатами для моделирования методом FDTD. Но есть одна проблема: уравнения там даны в виде соотношений, и для практического применения нужно выразить из них вычисляемый параметр. Вроде бы дело не сложное. Но из-за их громоздкости и запутанности пропустить ошибку в каком-нибудь индексе очень легко. На поиски же ее можно потратить не один день. Есть неплохая реализация этого метода на фортране. Тем не менее, в виде, пригодном для переноса в программу, мне эти формулы не встречались. Думаю, что тем, кто будет самостоятельно реализовывать PML, они могут пригодиться.

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

T_{yy1}^{i,j*}=T_{yy1}^{i,j}+\frac{\left ( \lambda +2\mu \right )}{\Delta y }\left ( \left (v_{y1}^{i+1,j}-v_{y1}^{i,j} \right )+\left (v_{y2}^{i+1,j}-v_{y2}^{i,j} \right ) \right )\Delta t


T_{zz1}^{i,j*}=T_{zz1}^{i,j}+\frac{ \lambda }{\Delta z }\left ( \left (v_{y1}^{i+1,j}-v_{y1}^{i,j} \right )+\left (v_{y2}^{i+1,j}-v_{y2}^{i,j} \right ) \right )\Delta t


T_{yy2}^{i,j*}=T_{yy2}^{i,j}+\frac{\lambda }{\Delta y }\left ( \left (v_{z1}^{i,j}-v_{z1}^{i,j-1} \right )+\left (v_{z2}^{i,j}-v_{z2}^{i,j-1} \right ) \right )\Delta t


T_{zz2}^{i,j*}=T_{zz2}^{i,j}+\frac{\left ( \lambda +2\mu \right )}{\Delta z }\left ( \left (v_{z1}^{i,j}-v_{z1}^{i,j-1} \right )+\left (v_{z2}^{i,j}-v_{z2}^{i,j-1} \right ) \right )\Delta t


T_{yz1}^{i,j*}=T_{yz1}^{i,j}+\frac{\mu}{\Delta z }\left ( \left (v_{y1}^{i,j}-v_{y1}^{i-1,j} \right )+\left (v_{y2}^{i,j}-v_{y2}^{i-1,j} \right ) \right )\Delta t


T_{yz2}^{i,j*}=T_{yz2}^{i,j}+\frac{\mu}{\Delta y }\left ( \left (v_{z1}^{i,j+1}-v_{z1}^{i,j} \right )+\left (v_{z2}^{i,j+1}-v_{z2}^{i,j} \right ) \right )\Delta t

v_{y1}^{i,j*}=v_{y1}^{i,j}+\frac{\left ( \left (T_{yy1}^{i,j}-T_{yy1}^{i-1,j} \right )+\left (T_{yy2}^{i,j}-T_{yy2}^{i-1,j} \right ) \right )\Delta t}{\rho\Delta y}


v_{y2}^{i,j*}=v_{y2}^{i,j}+\frac{\left ( \left (T_{yz1}^{i,j}-T_{yz1}^{i,j-1} \right )+\left (T_{yz2}^{i,j}-T_{yz2}^{i,j-1} \right ) \right )\Delta t}{\rho\Delta z}


v_{z1}^{i,j*}=v_{z1}^{i,j}+\frac{\left ( \left (T_{yz1}^{i+1,j}-T_{yz1}^{i,j} \right )+\left (T_{yz2}^{i+1,j}-T_{yz2}^{i,j} \right ) \right )\Delta t}{\rho\Delta y}


v_{z2}^{i,j*}=v_{z2}^{i,j}+\frac{\left ( \left (T_{zz1}^{i,j+1}-T_{zz1}^{i,j} \right )+\left (T_{zz2}^{i,j+1}-T_{zz2}^{i,j} \right ) \right )\Delta t}{\rho\Delta z}

Для участков в зоне PML:

T_{yy1}^{i,j*}=\frac{T_{yy1}^{i,j}\cdot \left ( \frac{1}{\Delta t} - d_{k}^{\frac{1}{2}}\right )+\frac{\left ( \lambda +2\mu \right )}{\Delta y }\left ( \left (v_{y1}^{i+1,j}-v_{y1}^{i,j} \right )+\left (v_{y2}^{i+1,j}-v_{y2}^{i,j} \right ) \right )\Delta t}{\left ( \frac{1}{\Delta t} + d_{k}^{\frac{1}{2}}\right )}


T_{zz1}^{i,j*}=\frac{T_{zz1}^{i,j}\cdot \left ( \frac{1}{\Delta t} - d_{k}^{\frac{1}{2}}\right )+\frac{ \lambda }{\Delta z }\left ( \left (v_{y1}^{i+1,j}-v_{y1}^{i,j} \right )+\left (v_{y2}^{i+1,j}-v_{y2}^{i,j} \right ) \right )\Delta t}{\left ( \frac{1}{\Delta t} + d_{k}^{\frac{1}{2}}\right )}


T_{yy2}^{i,j*}=\frac{T_{yy2}^{i,j}\cdot \left ( \frac{1}{\Delta t} - d_{k}\right )+\frac{\lambda }{\Delta y }\left ( \left (v_{z1}^{i,j}-v_{z1}^{i,j-1} \right )+\left (v_{z2}^{i,j}-v_{z2}^{i,j-1} \right ) \right )\Delta t}{\left ( \frac{1}{\Delta t} + d_{k}\right )}


T_{zz2}^{i,j*}=\frac{T_{zz2}^{i,j}\cdot \left ( \frac{1}{\Delta t} - d_{k}\right )+\frac{\left ( \lambda +2\mu \right )}{\Delta z }\left ( \left (v_{z1}^{i,j}-v_{z1}^{i,j-1} \right )+\left (v_{z2}^{i,j}-v_{z2}^{i,j-1} \right ) \right )\Delta t}{\left ( \frac{1}{\Delta t} + d_{k}\right )}


T_{yz1}^{i,j*}=\frac{T_{yz1}^{i,j}\cdot \left ( \frac{1}{\Delta t} - d_{k}\right )+\frac{\mu}{\Delta z }\left ( \left (v_{y1}^{i,j}-v_{y1}^{i-1,j} \right )+\left (v_{y2}^{i,j}-v_{y2}^{i-1,j} \right ) \right )\Delta t}{\left ( \frac{1}{\Delta t} + d_{k}\right )}


T_{yz2}^{i,j*}=\frac{T_{yz2}^{i,j}\cdot \left ( \frac{1}{\Delta t} - d_{k}^{\frac{1}{2}}\right )+\frac{\mu}{\Delta y }\left ( \left (v_{z1}^{i,j+1}-v_{z1}^{i,j} \right )+\left (v_{z2}^{i,j+1}-v_{z2}^{i,j} \right ) \right )\Delta t}{\left ( \frac{1}{\Delta t} + d_{k}^{\frac{1}{2}}\right )}

v_{y1}^{i,j*}=\frac{v_{y1}^{i,j}\cdot \left ( \frac{1}{\Delta t} - d_{k}\right )+\frac{\left ( \left (T_{yy1}^{i,j}-T_{yy1}^{i-1,j} \right )+\left (T_{yy2}^{i,j}-T_{yy2}^{i-1,j} \right ) \right )\Delta t}{\rho\Delta y}}{\left ( \frac{1}{\Delta t} + d_{k}\right )}


v_{y2}^{i,j*}=\frac{v_{y2}^{i,j}\cdot \left ( \frac{1}{\Delta t} - d_{k}\right )+\frac{\left ( \left (T_{yz1}^{i,j}-T_{yz1}^{i,j-1} \right )+\left (T_{yz2}^{i,j}-T_{yz2}^{i,j-1} \right ) \right )\Delta t}{\rho\Delta z}}{\left ( \frac{1}{\Delta t} + d_{k}\right )}


v_{z1}^{i,j*}=\frac{v_{z1}^{i,j}\cdot \left ( \frac{1}{\Delta t} - d_{k}^{\frac{1}{2}}\right )+\frac{\left ( \left (T_{yz1}^{i+1,j}-T_{yz1}^{i,j} \right )+\left (T_{yz2}^{i+1,j}-T_{yz2}^{i,j} \right ) \right )\Delta t}{\rho\Delta y}}{\left ( \frac{1}{\Delta t} + d_{k}^{\frac{1}{2}}\right )}


v_{z2}^{i,j*}=\frac{v_{z2}^{i,j}\cdot \left ( \frac{1}{\Delta t} - d_{k}^{\frac{1}{2}}\right )+\frac{\left ( \left (T_{zz1}^{i,j+1}-T_{zz1}^{i,j} \right )+\left (T_{zz2}^{i,j+1}-T_{zz2}^{i,j} \right ) \right )\Delta t}{\rho\Delta z}}{\left ( \frac{1}{\Delta t} + d_{k}^{\frac{1}{2}}\right )}

Где коэффициенты PML определяются:

d_{0}=\frac{3\cdot C_{p}\cdot\ln \left ( \frac{1}{R_{coef}} \right )}{2\cdot h}


d_{i}=d_{0}\cdot\left ( \frac{h-x_{i}}{h} \right )^{2}


d_{i}^{\frac{1}{2}}=d_{0}\cdot\left ( \frac{h-x_{i}-\frac{\Delta x}{2}}{h}\right )^{2}


Где:
C_{p} скорость продольной волны
R_{coef} желаемый коэффициент затухания
h толщина PML слоя
x_{i} координата точки относительно PML слоя

Запись опубликована в рубрике wave equation, моделирование с метками , . Добавьте в закладки постоянную ссылку.

Добавить комментарий