colin/fdtd-coremem
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Advance the E field by using Ampere's circuital law MACROSCOPIC equation

Browse Source
This commit is contained in:
Colin
2020-08-27 15:23:52 -07:00
parent f948820edf
commit 57258eb6cc
2 changed files with 34 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files

55
src/lib.rs
View File

@@ -25,7 +25,7 @@ pub mod consts {
/// Vacuum Permeability
pub const MU0: f64 = 1.2566370621219e-6; // H/m
// Vacuum Permittivity
// pub const Eps0: f64 = 8.854187812813e-12 // F⋅m−1
pub const EPS0: f64 = 8.854187812813e-12; // F⋅m−1
pub mod real {
use crate::R64;
pub fn C() -> R64 {
@@ -34,6 +34,9 @@ pub mod consts {
pub fn C2() -> R64 {
C() * C()
}
pub fn EPS0() -> R64 {
super::EPS0.into()
}
pub fn MU0() -> R64 {
super::MU0.into()
}
@@ -172,18 +175,23 @@ impl<M: Material> Cell<M> {
consts::MU0 * (self.hz() + self.mat.mz())
}
fn bz_to_hz(&self, bz: f64) -> f64 {
// B = mu0*(H + M) => H = B/mu0 - M
bz/consts::MU0 - self.mat.mz()
}
fn impulse_bz(&mut self, delta_bz: f64) {
self.state.hz = self.mat.step_h(&self.state, delta_bz).into();
}
fn step_h(mut self, right: &Self, down: &Self, _delta_t: R64) -> Self {
// ```tex
// Maxwell's equation: $\nabla x E = -dB/dt$
// Maxwell-Faraday equation: $\nabla x E = -dB/dt$
// Expand: $dE_y/dx - dE_x/dy = -dB_z/dt$
// Rearrange: $dB_z/dt = dE_x/dy - dE_y/dx$
// Discretize: $(\Delta B_z)/(\Delta t) = (\Delta E_x)/(\Delta y) - (\Delta E_y)/(\Delta x)$
//
// light travels C meters per second, and it we're only stepping half a timestep here (the other half when we advance E),
// light travels C meters per second, and we're only stepping half a timestep here (the other half when we advance E),
// therefore $(\Delta x)/(2 \Delta t) = (\Delta y)/(2 \Delta t) = c$ if we use SI units.
// Simplify: $\Delta B_z = (\Delta E_x)/(2c) - (\Delta E_y)/(2c)$
// ```
@@ -204,39 +212,40 @@ impl<M: Material> Cell<M> {
/// feature_size = how many units apart is the center of each adjacent cell on the grid.
fn step_e(self, left: &Self, up: &Self, delta_t: R64, feature_size: R64) -> Self {
// ```tex
// Ampere's circuital law with Maxwell's addition, in SI units:
// $\nabla x B = \mu_0 (J + \epsilon_0 dE/dt)$ where J = current density = $\sigma E$, $\sigma$ being a material parameter
// Expand: $\nabla x B = \mu_0 \sigma E + \mu_0 \epsilon_0 dE/dt$
// Substitute: $\nabla x B = \mu_0 \sigma E + 1/c^2 dE/dt$ where $c = 1/\sqrt{\mu_0 \epsilon_0}$ is the speed of light
// Rearrange: $dE/dt = c^2 (\nabla x B - \mu_0 \sigma E)$
// Expand: $dE_x/dt = c^2 (dB_z/dy - \mu_0 \sigma E_x)$ (1); $dE_y/dt = c^2 (-dB_z/dx - \mu_0 \sigma E_y)$ (2)
// Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"):
// $\nabla x H = J_f + dD/dt$ where $J_f$ = current density = $\sigma E$, $\sigma$ being a material parameter ("conductivity")
// Note that $D = \epsilon_0 E + P$, but we don't simulate any material where $P \ne 0$.
// Expand $D$: $\nabla x H = J_f + \epsilon_0 dE/dt$
// Expand $J$: $\nabla x H = \sigma E + \epsilon_0 dE/dt$
// Rearrange: $dE/dt = 1/\epsilon_0 (\nabla x H - \sigma E)$
// Expand: $dE_x/dt = 1/\epsilon_0 (dH_z/dy - \sigma E_x)$ (1); $dE_y/dt = 1/\epsilon_0 (-dH_z/dx - \sigma E_y)$ (2)
//
// Consider (1): let $E_p$ be $E_x$ at $T-\Delta t$ and $E_n$ be $E_x$ at $T+\Delta t$.
// Linear expansion about $t=T$, and discretized:
// $(E_n-E_p)/(2\Delta t) = c^2(\Delta B_z/\Delta y - \mu_0\sigma(E_n+E_p)/2)$
// Normalize: $E_n - E_p = 2\Delta{t} c^2 \Delta{B_z}/\Delta{y} - c^2 \mu_0 \sigma \Delta{t} (E_n + E_p)$
// Rearrange: $E_n(1 + c^2 \mu_0 \sigma \Delta{t}) = E_p(1 - c^2 \mu_0 \sigma \Delta{t}) + 2\Delta{t} c^2 \Delta{B_z}/\Delta{y}$
// $(E_n-E_p)/(2\Delta{t}) = 1/\epsilon_0(\Delta{H_z}/\Delta{y} - \sigma(E_n+E_p)/2)$
// Normalize: $E_n - E_p = 2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{y} - \sigma/\epsilon_0 \Delta{t} (E_n + E_p)$
// Rearrange: $E_n(1 + \sigma/\epsilon_0 \Delta{t}) = E_p(1 - \sigma/\epsilon_0 \Delta{t}) + 2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{y}$
// Then $E_n$ (i.e. the x value of $E$ after this step) is trivially solved
//
// Consider (2): let $E_p$ be $E_y$ at $T-\Delta t$ and $E_n$ be $E_y$ at $T+\Delta t$.
// Linear expansion about $t=T$, and discretized:
// $(E_n-E_p)/(2\Delta t) = c^2(-\Delta B_z/\Delta x - \mu_0\sigma(E_n+E_p)/2)$
// Normalize: $E_n - E_p = -2\Delta{t} c^2 \Delta{B_z}/\Delta{x} - c^2 \mu_0 \sigma \Delta{t} (E_n + E_p)$
// Rearrange: $E_n(1 + c^2 \mu_0 \sigma \Delta{t}) = E_p(1 - c^2 \mu_0 \sigma \Delta{t}) - 2\Delta{t} c^2 \Delta{B_z}/\Delta{x}$
// $(E_n-E_p)/(2\Delta{t}) = 1/\epsilon_0(-\Delta{H_z}/\Delta{x} - \sigma(E_n+E_p)/2)$
// Normalize: $E_n - E_p = -2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{x} - \sigma/\epsilon_0 \Delta{t} (E_n + E_p)$
// Rearrange: $E_n(1 + \sigma/\epsilon_0 \Delta{t}) = E_p(1 - \sigma/\epsilon_0 \Delta{t}) - 2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{x}$
// Then $E_n$ (i.e. the y value of $E$ after this step) is trivially solved
// ```
use consts::real::{C, C2, HALF, MU0, ONE, TWO};
use consts::real::{C, C2, EPS0, HALF, ONE, TWO};
let sigma = self.mat.conductivity();
let sigma: R64 = self.mat.conductivity().into();
let delta_bz_y = self.bz() - up.bz();
let ex_rhs = self.state.ex*(ONE() - C2()*MU0()*sigma*delta_t) + TWO()*delta_t*C2()*delta_bz_y/feature_size;
let ex_next = ex_rhs / (ONE() + C2()*MU0()*sigma*delta_t);
let delta_hz_y = self.hz() - self.bz_to_hz(up.bz());
let ex_rhs = self.state.ex*(ONE() - sigma/EPS0()*delta_t) + TWO()*delta_t/EPS0()*delta_hz_y/feature_size;
let ex_next = ex_rhs / (ONE() + sigma/EPS0()*delta_t);
let delta_bz_x = self.bz() - left.bz();
let ey_rhs = self.state.ey*(ONE() - C2()*MU0()*sigma*delta_t) - TWO()*delta_t*C2()*delta_bz_x/feature_size;
let ey_next = ey_rhs / (ONE() + C2()*MU0()*sigma*delta_t);
let delta_hz_x = self.hz() - self.bz_to_hz(left.bz());
let ey_rhs = self.state.ey*(ONE() - sigma/EPS0()*delta_t) - TWO()*delta_t/EPS0()*delta_hz_x/feature_size;
let ey_next = ey_rhs / (ONE() + sigma/EPS0()*delta_t);
Cell {
state: CellState {

3
src/render.rs
View File

@@ -53,8 +53,9 @@ impl ColorTermRenderer {
let r = norm_color(cell.mat.mz());
let b = (55.0*cell.mat().conductivity()).min(255.0) as u8;
//let b = 0;
//let g = norm_color(cell.ex());
//let b = norm_color(cell.ey());
//let g = 0;
//let g = norm_color(cell.ex());
//let g = norm_color(curl(cell.ex(), cell.ey()));
let g = norm_color((cell.bz() * 3.0e8).into());
//let g = norm_color(cell.ey().into());