FDTD در OghmaNano

1. مقدمه

Finite-Difference Time-Domain (FDTD) یک روش عددی تمام‌موج برای حل معادلات ماکسول به‌صورت مستقیم در حوزه زمان روی یک شبکه مکانی گسسته است. برخلاف رهگیری پرتو، حل‌گرهای مُد، یا روش‌های transfer-matrix، روش FDTD هیچ فرض هندسی یا مُدی اعمال نمی‌کند. در عوض، میدان‌های الکتریکی و مغناطیسی با گسسته‌سازی صریح معادلات کرل—قانون آمپر و قانون فاراده—در زمان پیش برده می‌شوند، به‌گونه‌ای که امواج الکترومغناطیسی به‌طور طبیعی در دامنه شبیه‌سازی منتشر شوند، پراکنده شوند، تداخل کنند و دچار پراش شوند.

\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]

در روش Finite-Difference Time-Domain (FDTD)، معادلات ماکسول به‌طور مستقیم در حوزه زمان حل می‌شوند، با جایگزین کردن مشتقات مکانی با تفاضل‌های محدود روی یک شبکه پله‌ای (معمولاً شبکه Yee) و پیش‌بردن میدان‌ها در زمان. از آن‌جا که این رویکرد عملاً هیچ فرضی درباره جهت انتشار، ساختار مُدی، یا رفتار میدان اعمال نمی‌کند، یک شبیه‌سازی منفرد FDTD می‌تواند پاسخ پهن‌باند، پراش و تداخل، اثرات میدان نزدیک، هندسه زیرطول‌موج، و پدیده‌های گذرای الکترومغناطیسی را در یک چارچوب ثبت کند. همین نبود فرض‌های ساده‌ساز است که FDTD را تا این حد عمومی می‌سازد، اما همین ویژگی نیز آن را از نظر محاسباتی پرهزینه می‌کند: شبکه مکانی باید کوچک‌ترین مقیاس طولی فیزیکی را تفکیک کند، گام زمانی به‌شدت توسط پایداری عددی محدود می‌شود، و شبیه‌سازی‌های سه‌بعدی حتی برای هندسه‌هایی که از نظر اندازه ظاهراً متعادل هستند نیز خیلی سریع از نظر حافظه و زمان پرهزینه می‌شوند.

2. چه زمانی از FDTD استفاده کنیم و چه زمانی از FDTD استفاده نکنیم

در اغلب مسائل عملی اپتیکی، FDTD نباید نخستین روشی باشد که به سراغ آن می‌روید. برای سامانه‌های نوری در مقیاس بزرگ، فاصله‌های انتشار طولانی، یا موقعیت‌هایی که اپتیک هندسی بر آن‌ها غالب است، FDTD انتخابی بسیار ناکارآمد است. سامانه‌های تصویربرداری، لنزهای ماکروسکوپی، روزنه‌ها، و مجموعه‌های نوری با ابعاد میلی‌متری تا سانتی‌متری بسیار بهتر است با رهگیری پرتو بررسی شوند، در حالی که ساختارهای چندلایه تخت مانند سلول‌های خورشیدی لایه‌نازک، پشته‌های OLED، و پوشش‌ها با روش‌های transfer-matrix بسیار دقیق‌تر و کارآمدتر مدل می‌شوند. در این نواحی، FDTD مانند یک پتک محاسباتی عمل می‌کند: به‌طور چشمگیری زمان اجرا و مصرف حافظه را افزایش می‌دهد، در حالی که بینش فیزیکی اضافه‌ای بسیار اندک یا هیچ به‌دست نمی‌دهد.

در عمل، FDTD باید به‌عنوان یک روش آخرین چاره در نظر گرفته شود، که برای مواردی نگه داشته می‌شود که مدل‌های ساده‌تر واقعاً از کار می‌افتند. نقطه قوت آن در موقعیت‌هایی است که میدان کامل الکترومغناطیسی باید به‌طور صریح حل شود، مانند پراش از شبکه‌های زیرطول‌موج، بلورهای فوتونیکی، متاسطوح، فصل‌مشترک‌های به‌شدت بی‌نظم یا زبر، و اثرات کوپله میدان نزدیک که با اپتیک پرتو یا مدل‌های محیط لایه‌ای قابل ثبت نیستند. حتی در این موارد نیز، FDTD باید به‌صورت محدود و هدفمند به‌کار گرفته شود، با درکی روشن از هزینه عددی، قیود پایداری، و نیازمندی‌های حافظه آن. بنابراین بخش‌های بعدی معادلات به‌روزرسانی گسسته مورد استفاده در OghmaNano را مستقیماً از قوانین آمپر و فاراده استخراج می‌کنند، نه برای تشویق به استفاده بی‌رویه، بلکه برای روشن کردن این‌که FDTD چه کاری می‌تواند انجام دهد و چرا باید با احتیاط از آن استفاده شود.

3. پس‌زمینه نظری

این بخش از راهنما با هدف توصیف کامل کد FDTD همراه با استخراج‌های مفصل برای کمک به فهم/یافتن خطاها نوشته شده است.

قانون آمپر به‌صورت \[\sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t} = \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ H_{x} & H_{y} & H_{z} \end{vmatrix}\]

که می‌توان آن را به‌صورت زیر بسط داد

\[\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} = \frac{\partial H_{z}}{\partial y}-\frac{\partial H_{y}}{\partial z}\]

\[\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\]

\[\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y}\]

برای حالت \(\frac{\partial}{\partial y}=0\)

\[\begin{split} &\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} =-\frac{\partial H_{y}}{\partial z}\\ &\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\\ &\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x} \end{split}\]

برای \(E_{x}\) \[\begin{split} &\sigma E_{x} + \epsilon \frac{\partial E_{x}}{\partial t} =-\frac{\partial H_{y}}{\partial z}\\ &\sigma \frac{E_{x}^{t+1}[]+E_{x}^{t}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]-E_{x}^{t}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}\\ &\sigma \frac{E_{x}^{t+1}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ &\sigma \frac{E_{x}^{t+1}[]}{2} + \epsilon \frac{E_{x}^{t+1}[]}{\Delta t} = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ & \frac{\sigma \Delta t + 2 \epsilon }{ 2 \Delta t}E_{x}^{t+1}[] = -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t}\\ & E_{x}^{t+1}[] = \left ( -\frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{x}^{t}[]}{2}+\epsilon \frac{E_{x}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon} \end{split}\]

برای \(E_{y}\) \[\begin{split} &\sigma E_{y} + \epsilon \frac{\partial E_{y}}{\partial t} = -\frac{\partial H_{z}}{\partial x}+\frac{\partial H_{x}}{\partial z}\\ &\sigma \frac{E_{y}^{t+1}[]+E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t+1}[]-E_{y}^{t}[]}{\Delta t} = -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}\\ &\sigma \frac{E_{y}^{t+1}[]}{2} + \epsilon \frac{E_{y}^{t+1}[]}{\Delta t} = -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t}[]}{\Delta t}\\ &E_{y}^{t+1}[] = \left ( -\frac{H_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}+\frac{H_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z}-\sigma \frac{E_{y}^{t}[]}{2} + \epsilon \frac{E_{y}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon} \end{split}\]

برای \(E_{z}\) \[\begin{split} &\sigma E_{z} + \epsilon \frac{\partial E_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x}\\ &\sigma \frac{E_{z}^{t+1}[]+E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t+1}[]-E_{z}^{t}[]}{\Delta t} = \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}\\ &\sigma \frac{E_{z}^{t+1}[]}{2} + \epsilon \frac{E_{z}^{t+1}[]}{\Delta t} = \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\sigma \frac{E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t}[]}{\Delta t}\\ &E_{z}^{t+1}[]= \left ( \frac{H_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-H_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\sigma \frac{E_{z}^{t}[]}{2} + \epsilon \frac{E_{z}^{t}[]}{\Delta t} \right ) \frac{2 \Delta t}{\sigma \Delta t + 2 \epsilon}\\ \end{split}\]

قانون فاراده به‌صورت \[-\sigma_{m} \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t} = \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ E_{x} & E_{y} & E_{z} \end{vmatrix}\]

که می‌توان آن را بسط داد تا به‌صورت زیر به‌دست آید:

\[-\sigma_{m} H_{x} - \mu \frac{\partial H_{x}}{\partial t} = \frac{\partial E_{z}}{\partial y}-\frac{\partial E_{y}}{\partial z}\]

\[-\sigma_{m} H_{y} - \mu \frac{\partial H_{y}}{\partial t} = -\frac{\partial E_{z}}{\partial x}+\frac{\partial E_{x}}{\partial z}\]

\[-\sigma_{m} H_{z} - \mu \frac{\partial H_{z}}{\partial t} = \frac{\partial E_{y}}{\partial x}-\frac{\partial E_{x}}{\partial y}\]

با \(\sigma_m=0\) و \(\frac{\partial}{\partial y}=0\)

\[\begin{split} &\frac{\partial H_{x}}{\partial t} = \frac{1}{\mu} \left ( \frac{\partial E_{y}}{\partial z} \right )\\ &\frac{\partial H_{y}}{\partial t} = \frac{1}{\mu} \left ( \frac{\partial E_{z}}{\partial x}-\frac{\partial E_{x}}{\partial z} \right )\\ &\frac{\partial H_{z}}{\partial t} = - \frac{1}{\mu} \left ( \frac{\partial E_{y}}{\partial x} \right ) \end{split}\]

که با گسسته‌سازی به‌صورت زیر درمی‌آید

\[\begin{split} & H_{x}^{t+1} = \frac{1}{\mu} \left ( \frac{E_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z} \right ) \Delta t + H_{x}^{t}[]\\ & H_{y}^{t+1} = \frac{1}{\mu} \left ( \frac{E_{z}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{z}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x}-\frac{E_{x}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{x}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta z} \right ) \Delta t+ H_{y}^{t}[]\\ & H_{z}^{t+1} = \frac{1}{\mu} \left ( - \frac{E_{y}^{t+\frac{1}{2}}[\frac{1}{2}]-E_{y}^{t+\frac{1}{2}}[-\frac{1}{2}]}{\Delta x} \right ) \Delta t + H_{x}^{z}[] \end{split}\]

4. حل عددی روی رایانه

معادلاتی که در بالا استخراج شدند تحول زمانی و مکانی پیوسته میدان‌های الکتریکی و مغناطیسی را توصیف می‌کنند. برای حل آن‌ها روی یک رایانه، هم فضا و هم زمان گسسته می‌شوند. دامنه شبیه‌سازی به یک شبکه منظم از سلول‌ها تقسیم می‌شود و مؤلفه‌های میدان در مکان‌های مکانی گسسته و گام‌های زمانی گسسته ذخیره می‌شوند. در فرمول‌بندی استاندارد FDTD، میدان‌های الکتریکی و مغناطیسی هم در فضا و هم در زمان به‌صورت پله‌ای قرار می‌گیرند، به‌گونه‌ای که هر مؤلفه میدان دقیقاً در جایی ارزیابی شود که به‌طور طبیعی توسط معادلات کرل مورد نیاز است.

در عمل، این به آن معناست که مؤلفه‌های میدان الکتریکی در گام‌های زمانی صحیح (\(t, t+\Delta t, t+2\Delta t\)) ذخیره می‌شوند، در حالی که مؤلفه‌های میدان مغناطیسی در گام‌های زمانی نیمه‌صحیح (\(t+\tfrac{1}{2}\Delta t, t+\tfrac{3}{2}\Delta t\)) ذخیره می‌شوند. مشتقات مکانی مانند \(\partial / \partial x\) و \(\partial / \partial z\) با تفاضل‌های محدود مرکزی بین نقاط شبکه همسایه جایگزین می‌شوند. این آرایش پله‌ای به یک طرح به‌روزرسانی صریح و ساده منجر می‌شود که در آن میدان‌ها در زمان به‌صورت leapfrog از روی یکدیگر عبور می‌کنند.

در هر گام زمانی، ابتدا میدان مغناطیسی با استفاده از فرم گسسته قانون فاراده به‌روزرسانی می‌شود، به‌طوری که کرل میدان الکتریکی از مقادیر فعلی میدان الکتریکی محاسبه می‌گردد. پس از آن‌که میدان مغناطیسی به‌اندازه نیم گام زمانی پیش برده شد، میدان الکتریکی با استفاده از فرم گسسته قانون آمپر، همراه با پارامترهای ماده مانند گذردهی، تراوایی، و رسانندگی، به‌روزرسانی می‌شود. هیچ حل ماتریسی سراسری لازم نیست: هر مؤلفه میدان در هر سلول شبکه فقط با استفاده از همسایه‌های محلی خود از گام زمانی قبلی به‌روزرسانی می‌شود.

این رویه گام‌برداری زمانی صریح تا رسیدن به زمان شبیه‌سازی مورد نظر تکرار می‌شود. چون این روش در حوزه زمان است، یک شبیه‌سازی منفرد به‌طور طبیعی حاوی اطلاعاتی در یک بازه فرکانسی پهن است؛ کمیت‌های حوزه فرکانس مانند طیف‌ها یا توزیع‌های میدان حالت پایدار معمولاً با ثبت مقادیر میدان به‌عنوان تابعی از زمان و سپس پس‌پردازش آن‌ها با استفاده از تبدیل‌های فوریه به‌دست می‌آیند. منابع، مرزها، و لایه‌های جذب‌کننده (مانند perfectly matched layers) به‌طور مستقیم در معادلات به‌روزرسانی در مکان‌های مناسب شبکه گنجانده می‌شوند.

سادگی معادلات به‌روزرسانی، FDTD را از نظر مفهومی ساده و به‌شدت موازی‌پذیر می‌سازد، اما این روش هم‌زمان قیود عددی سخت‌گیرانه‌ای نیز اعمال می‌کند. گام زمانی باید یک شرط پایداری را ارضا کند که به فاصله شبکه مکانی و خواص ماده بستگی دارد، و ردپای حافظه مستقیماً با تعداد سلول‌های شبکه و مؤلفه‌های میدان ذخیره‌شده متناسب است. این قیود در نهایت اندازه و مدت شبیه‌سازی‌هایی را که در عمل می‌توان انجام داد محدود می‌کنند و انتخاب دقیق تفکیک‌پذیری شبکه، اندازه دامنه، و پیچیدگی مدل را ضروری می‌سازند.