CN105760596A - Two-dimensional vacuum Crank-Nicolson complete matching layer implementation algorithm based on auxiliary differential equation - Google Patents

Abstract

The invention relates to a two-dimensional vacuum Crank-Nicolson complete matching layer implementation algorithm based on an auxiliary differential equation and belongs to the technical field of numerical simulation.The method aims at reducing a vacuum FDTD computational domain and simulating limited memory space of a computer into infinite space.The implementation algorithm is technically characterized in that in the process that a two-dimensional modified Maxwell equation with plural stretching coordinate variables is converted into the time domain finite difference from the frequency domain, a Douglas-Gunn (DG) algorithm is utilized, an iteration equation with coefficients being block tridiagonal matrixes is approximately decomposed into two iteration equations with coefficients being tridiagonal matrixes, wherein the two iteration equations can be efficiently solved, and computational efficiency is obviously improved.The implementation algorithm has the advantages of achieving unconditional stability, increasing the electromagnetic field computational speed and saving memory.

Description

Two-dimentional vacuum C rank-Nicolson completely permutation based on auxiliary differential equation realizes Algorithm

Technical field

The present invention relates to technical field of value simulation, particularly to one based on auxiliary differential equation Two-dimentional vacuum C rank-Nicolson completely permutation realize algorithm.

Background technology

Time-domain finite difference (FDTD) widely should as a kind of electromagnetic method that calculates In the Electromagnetic Simulation of various time domains calculates, such as antenna, radio circuit, optics and half Conductor etc..FDTD has wide applicability, is suitable for parallel computation, calculation procedure versatility Etc. feature.

But, along with going deep into and various more and more wide variety of demands of scientific research, it is calculated Method originally experience Courant Friedrichs Lewy (CFL) the numerical stability conditions restriction lack Fall into more and more obvious.The suffered numerical stability condition of algorithm itself limits: time during calculating Between step-length and spatial mesh size must be fulfilled for CFL constraints, i.e.

$\mathrm{\Δ}t\≤{\left(c_0\sqrt{{\left(\mathrm{\Δ}x\right)}^{-2}+{\left(\mathrm{\Δ}y\right)}^{-2}+{\left(\mathrm{\Δ}z\right)}^{-2}}\right)}^{-1}$

In formula, Δ t is for calculating time step, c₀For the free space light velocity, Δ x, Δ y and Δ z are Three dimensions step-length.In Practical Calculation, spatial spreading step-length and time step relative wavelength and Cycle is the least, so the feelings of inadequate resource will necessarily be occurred when calculating Electrically large size object Condition, the computational efficiency causing FDTD is the lowest.Therefore to eliminate the restriction of CFL condition, The alternating direction implicit (Altemating Direction Implicit, ADI) of unconditional stability FDTD method, local one-dimension (Local One Dimension, LOD) FDTD method and Crane gram Nicholson (Crank-Nicolson, CN) FDTD method is suggested in succession.

Although overcoming to a certain extent for ADI-FDTD algorithm and LOD-FDTD algorithm Stability condition restriction, but the computational accuracy of algorithm is too low, and performance is unsatisfactory, its reason It is owing to, after time step increases, the numerical dispersion caused increases, and then causes the mistake of algorithm Difference is bigger.2004, G.Sun et al. used Crank-Nicolson difference scheme to Max Wei Fangcheng carries out sliding-model control, i.e. CN-FDTD, and algorithm is much larger than in time step value Stability condition (such as 20 times) remains to keep good lasting accuracy, shows the suitableeest By property, and CN-FDTD algorithm is a kind of method of easier unconditional stability, will 2 above required in two kinds of algorithms calculating processes are simplified to 1 calculating process, thus significantly Reducing calculation resources, therefore scholars unanimously think that CN-FDTD has broader development Prospect.

Due to the restriction in calculator memory space, numerical computations can only be carried out in limited region, Open or the problem such as the electromagnetic radiation in semi-open region and scattering to be able to simulation, calculating district Absorbing boundary condition is must be provided with, in order to limited mesh space simulation at the cutoff boundary in territory Open infinite space, solves the Electromagnetic Wave Propagation in arbitrary medium and various electromagnetic problem. The completely permutation (Perfectly Matched Layer, PML) proposed by Berenger is mesh The absorbing boundary condition that front application is wider, PML is it is to be understood that pass through in FDTD region A kind of special media layer, the wave impedance of this layer of medium and adjacent media wave resistance are set at cutoff boundary Resist and mate completely, so that incidence wave areflexia ground enters PML layer through interface, PML Layer is lossy dielectric, finally by electro-magnetic wave absorption.The most conventional PML absorbing boundary is main There are stretching coordinate transform completely permutation (SC-PML) and unit anisotropy completely permutation (UPML)。

Summary of the invention

It is an object of the invention to for FDTD algorithm by lacking that CFL stability condition is limited Fall into, improve and block the computational efficiency of PML algorithm of two dimension vacuum and assimilation effect and propose SC-PML algorithm based on Sub-ODE method and CN-FDTD.This algorithm combines Douglas-Gunn solves thought, can significantly improve asking of CN-FDTD-PML algorithm Solve efficiency.

Two-dimentional vacuum C rank-Nicolson completely permutation based on auxiliary differential equation realizes calculating Method, comprises the following steps:

Step 1: Maxwell equation in frequency domain is modified to the Mike with stretching coordinate operator This Wei Fangcheng, and represent in rectangular coordinate system；

Step 2: according to the mapping transformation relation of frequency domain and time domain, by two in rectangular coordinate system Keep in repair positive Maxwell equation and transform to time-domain representation, be simultaneously based on Sub-ODE method Auxiliary variable is set；

Step 3: time domain expanded form based on Crank-Nicolson Finite Difference Time Domain, Two dimension Maxwell equation in the rectangular coordinate system of forms of time and space is launched into Fdtd Method Form, also auxiliary differential equation is transformed to the form of Fdtd Method simultaneously；

Step 4: Finite Difference-Time Domain form-separating is organized into the form solved, result produces one group The coupling implicit equation in electric field and magnetic field, collated after to obtain coefficient matrix be block tridiagonal matrix The implicit iterative equation of form；

Step 5: use Crank-Nicolson Douglas-Gunn (CNDG) method, will step Electric field iterative equation obtained by rapid 4 be approximately decomposed into can with Efficient Solution, coefficient be three right Two iterative equations of angular moment battle array；

Step 6: utilize the iterative equation obtained by step 5 to solve electric field value, then will solve The electric field value gone out is updated in the iterative equation in magnetic field, solves magnetic-field component, by solve Electric field value and magnetic field value are updated in the iterative equation of auxiliary variable, solve the value of auxiliary variable；

Repetition step 6, thus iterative in time.

Use CNDG method can with effectively by coefficient matrix for block tridiagonal matrix form Electric field iterative equation is approximately decomposed into can with Efficient Solution, coefficient as triple diagonal matrix two Iterative equation, reduces computation complexity, improves computational efficiency, has guidance meaning to FDTD algorithm Justice.

Accompanying drawing illustrates:

Fig. 1 is FB(flow block) of the present invention；

Relative to reflection error figure when Fig. 2 is CFLN difference of the present invention.

Detailed description of the invention:

The purport of the present invention is to propose a kind of two-dimentional vacuum based on auxiliary differential equation Crank-Nicolson completely permutation realizes algorithm, utilizes Douglas-Gunn to solve thought pole The earth improves Electromagnetic Calculation speed.

Embodiment of the present invention is described further in detail below in conjunction with the accompanying drawings.

Fig. 1 is flow chart of the present invention, implements step as follows:

Step 1: Maxwell equation in frequency domain is modified to the Mike with stretching coordinate operator This Wei Fangcheng, and Maxwell equation revised in frequency domain is represented in rectangular coordinate system, The dissemination in a vacuum of TM ripple can be described as

$-{\mathrm{j\ωH}}_x=c_0\·S_y^{-1}\frac{\∂E_z}{\∂y}---\left(1\right)$

${\mathrm{j\ωH}}_y=c_0\·S_x^{-1}\frac{\∂E_z}{\∂x}---\left(2\right)$

${\mathrm{j\ωE}}_z=c_0\·S_x^{-1}\frac{\∂H_y}{\∂x}-c_0\·S_y^{-1}\frac{\∂H_x}{\∂y}---\left(3\right)$

In formula, S_η(η=x, y) for stretching coordinate variable, can be expressed as

$S_{\mathrm{\η}}^{-1}=\frac{1}{1+{\mathrm{\σ}}_{\mathrm{\η}}/\left({\mathrm{j\ω\ϵ}}_0\right)}=1-\frac{{\mathrm{\σ}}_{\mathrm{\η}}/{\mathrm{\ϵ}}_0}{j\mathrm{\ω}+{\mathrm{\σ}}_{\mathrm{\η}}/{\mathrm{\ϵ}}_0}---\left(4\right)$

Step 2: according to the mapping transformation relation of frequency domain and time domain, by two in rectangular coordinate system Keep in repair positive Maxwell equation and transform to time-domain representation, be simultaneously based on Sub-ODE method Auxiliary variable is set, i.e.

$-\frac{\∂}{\∂t}{H}_x=c\·\frac{\∂E_z}{\∂y}-g_{xy}---\left(5\right)$

$\frac{\∂}{\∂t}{H}_y=c\·\frac{\∂E_z}{\∂x}-g_{yx}---\left(6\right)$

$\frac{\∂}{\∂t}{E}_z=c\·\frac{\∂H_y}{\∂x}-c\·\frac{\∂H_x}{\∂y}-f_{zx}+f_{zy}---\left(7\right)$

In formula, f_zx、f_zy、g_xyAnd g_yxFor auxiliary variable.

Step 3: time domain expanded form based on Crank-Nicolson Finite Difference Time Domain, Two dimension Maxwell equation in the rectangular coordinate system of forms of time and space is launched into Fdtd Method Form, is also transformed to the form of Fdtd Method simultaneously by time domain auxiliary differential equation, utilizes CN item is by formula (5)-(7) discretization, and can obtain discrete equation is

$H_{x_{i,j+1/2}}^{n+1}=H_{x_{i,j+1/2}}^n-{\mathrm{\χ}}_y\·{\mathrm{\Γ}}_y\left(E_{z_{i,j}}^n\right)+{\mathrm{\Δt}}_h\·(g_{{\mathrm{xy}}_{i,j+1/2}}^{n+1}+g_{{\mathrm{xy}}_{i,j+1/2}}^n)---\left(8\right)$

$H_{y_{i+1/2,j}}^{n+1}=H_{y_{i+1/2,j}}^n+{\mathrm{\χ}}_x\·{\mathrm{\Γ}}_x\left(E_{z_{i,j}}^n\right)-{\mathrm{\Δt}}_h\·(g_{{\mathrm{yx}}_{i+1/2,j}}^{n+1}+g_{{\mathrm{yx}}_{i+1/2,j}}^n)---\left(9\right)$

$\begin{array}{c}{E}_{z_{i,j}}^{n+1}=E_{z_{i,j}}^n+{\mathrm{\χ}}_x\·{\mathrm{\Γ}}_x\left(H_{y_{i+1/2,j}}^n\right)-{\mathrm{\χ}}_y\·{\mathrm{\Γ}}_y\left(H_{x_{i,j+1/2}}^n\right)\\ -{\mathrm{\Δt}}_h\·\left(f_{{\mathrm{zx}}_{i,j}}^{n+1}+f_{{\mathrm{zx}}_{i,j}}^n\right)+{\mathrm{\Δt}}_h\·\left(f_{{\mathrm{zy}}_{i,j}}^{n+1}+f_{{\mathrm{zy}}_{i,j}}^n\right)\end{array}---\left(10\right)$

In formula, Δ t_h=Δ t/2, χ_η=c Δ t/ (2 Δ η), Δ η (η=x, y) is space cell size, In order to clear, Г_η[*] is the shorthand in CN method, as

${\mathrm{\Γ}}_y\left(E_{z_{i,j}}^n\right)=E_{z_{i,j+1}}^{n+1}-E_{z_{i,j}}^{n+1}+E_{z_{i,j+1}}^n-E_{z_{i,j}}^n---\left(11\right)$

Step 4: the form of Fdtd Method is organized into the form solved, result produces one Group electric field and the coupled wave equation in magnetic field, this is one group of implicit equation, this prescription journey is decoupled, whole The coefficient electric field implicit iterative equation that the left side is block tridiagonal matrix form is obtained after reason

$\begin{array}{c}\[1-\left(D_{2x}+D_{2y}\right)\]E_{z_{i,j}}^{n+1}=\[1+\left(D_{2x}+D_{2y}\right)\]E_{z_{i,j}}^n\\ +2\left({\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ex}}_i}\right)\left(H_{y_{i+1/2,j}}^n-H_{y_{i-1/2,j}}^n\right)-{\mathrm{\Δt}}_h\left(1+r_{0{\mathrm{ex}}_i}\right)f_{{\mathrm{zx}}_{i,j}}^n\\ -2\left({\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ey}}_j}\right)\left(H_{x_{i,j+1/2}}^n-H_{x_{i,j-1/2}}^n\right)+{\mathrm{\Δt}}_h\left(1+r_{0{\mathrm{ey}}_j}\right)f_{{\mathrm{zy}}_{i,j}}^n\\ -\left({\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ex}}_i}\right){\mathrm{\Δt}}_h\[\left(1+r_{0{\mathrm{hx}}_{i+1/2}}\right)g_{{\mathrm{yx}}_{i+1/2,j}}^n-\left(1+r_{0{\mathrm{hx}}_{i+1/2}}\right)g_{{\mathrm{yx}}_{i-1/2,j}}^n\]\\ -\left({\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ey}}_j}\right){\mathrm{\Δt}}_h\[\left(1+r_{0{\mathrm{hy}}_{j+1/2}}\right)g_{{\mathrm{xy}}_{i,j+1/2}}^n-\left(1+r_{0{\mathrm{hy}}_{j-1/2}}\right)g_{{\mathrm{xy}}_{i,j-1/2}}^n\]\end{array}---\left(12\right)$

In formula, D_2xWith D_2yCan be defined as

$\begin{array}{c}{D}_{2x}{E}_{z_{i,j}}^{n+1}=\left({\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ex}}_i}\right)\left({\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{hx}}_{i,j+1/2}}\right)E_{z_{i+1,j}}^{n+1}\\ -\left({\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ex}}_i}\right)\left({\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{hx}}_{i+1/2}}+{\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{hx}}_{i-1/2}}\right)E_{z_{i,j}}^{n+1}\\ +\left({\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ex}}_i}\right)\left({\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{hx}}_{i-1/2}}\right)E_{z_{i-1,j}}^{n+1}\end{array}---\left(13\right)$

$\begin{array}{c}{D}_{2y}{E}_{z_{i,j}}^{n+1}=\left({\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ey}}_j}\right)\left({\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{hy}}_{j+1/2}}\right)E_{z_{i,j+1}}^{n+1}\\ -\left({\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ey}}_j}\right)\left({\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{hy}}_{j+1/2}}+{\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{hy}}_{j-1/2}}\right)E_{z_{i,j}}^{n+1}\\ +\left({\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ey}}_j}\right)\left({\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{hy}}_{j-1/2}}\right)E_{z_{i,j-1}}^{n+1}\end{array}---\left(14\right)$

Step 5: formula (12) remains a more complicated matrix, it is still desirable to the biggest Amount of calculation, uses CNDG method, adds respectively on the left side of formula (12) and the rightWithArrange

$\begin{array}{c}(1-D_{2x})(1-D_{2y})E_{z_{i,j}}^{n+1}=(1+D_{2x})(1+D_{2y})E_{z_{i,j}}^n\\ +f\left(H_x^n,H_y^n,f_{zx}^n,f_{zy}^n,g_{xy}^n,g_{yx}^n\right)\end{array}---\left(15\right)$

In formula,For remaining n moment field amount and the writing a Chinese character in simplified form of auxiliary variable Form.Formula (15) can be tried to achieve by following two formulas

$(1-D_{2x})E_{z_{i,j}}^{*}=(1+D_{2x}+2D_{2y})E_{z_{i,j}}^n+f(H_x^n,H_y^n,f_{zx}^n,f_{zy}^n,g_{xy}^n,g_{yx}^n)---\left(16\right)$

$(1-D_{2y})E_{z_{i,j}}^{n+1}=E_{z_{i,j}}^{*}-D_{2y}{E}_{z_{i,j}}^n---\left(17\right)$

Step 6: utilize the iterative equation obtained by step 5 to solve electric field value, then will solve The electric field value gone out is updated in the iterative equation in magnetic field, solves magnetic-field component, by solve Electric field value and magnetic field value are updated in the iterative equation of auxiliary variable, solve the value of auxiliary variable；

Repetition step 6, thus iterative in time.

Relative to reflection error figure when Fig. 2 is CFLN difference of the present invention, the extracting method in order to verify, Inventive algorithm is programmed, obtains result shown in Fig. 2 by Computer Simulation, wherein,In formulaIt is that conventional FDTD algorithm ensure that numerical stability Discrete interval maximum time of property.As seen from Figure 2, the absorbability of CNDG-SC-PML Can not change with the increase of CNLN, illustrate that this algorithm has unconditional stability, imitative True process required time is shorter compared with traditional algorithm simulation time.

The foregoing is only presently preferred embodiments of the present invention, be not limiting as the present invention, all at this Within bright spirit and principle, any modification, equivalent substitution and improvement etc. made, all should wrap Within being contained in protection scope of the present invention.

Claims (4)

1. two-dimentional vacuum C rank-Nicolson completely permutation based on auxiliary differential equation realizes Algorithm, comprises the following steps:

Step 1: Maxwell equation in frequency domain is modified to the Maxwell with stretching coordinate operator Equation, and represent in rectangular coordinate system；

Step 2: according to the mapping transformation relation of frequency domain and time domain, by two maintenances in rectangular coordinate system Positive Maxwell equation transforms to time-domain representation, is simultaneously based on Sub-ODE method and arranges Auxiliary variable；

Step 3: time domain expanded form based on Crank-Nicolson Finite Difference Time Domain, by time In the rectangular coordinate system of territory form, two dimension Maxwell equation is launched into the shape of Fdtd Method Formula, is also transformed to the form of Fdtd Method simultaneously by auxiliary differential equation；

Step 4: Finite Difference-Time Domain form-separating is organized into the form solved, result produces one group of electric field With the coupling implicit equation in magnetic field, collated after to obtain coefficient matrix be block tridiagonal matrix form Implicit iterative equation；

Step 5: use Crank-Nicolson Douglas-Gunn method, by the electricity obtained by step 4 Iterative equation be approximately decomposed into can with Efficient Solution, coefficient as triple diagonal matrix two repeatedly For equation；

Step 6: utilize the iterative equation obtained by step 5 to solve electric field value, then will solve Electric field value is updated in the iterative equation in magnetic field, solves magnetic-field component, the electric field that will solve Value and magnetic field value are updated in the iterative equation of auxiliary variable, solve the value of auxiliary variable；

Repetition step 6, thus iterative in time.

2. according to based on auxiliary differential equation the two-dimentional vacuum described in right 1 Crank-Nicolson completely permutation realizes algorithm, it is characterised in that: step 2, by revise Maxwell equation transforms to time domain, is simultaneously based on Sub-ODE method and arranges auxiliary variable f_zx、f_zy、g_xyAnd g_yx。

3. according to based on auxiliary differential equation the two-dimentional vacuum described in right 1 Crank-Nicolson completely permutation realizes algorithm, it is characterised in that: step 3, based on The time domain expanded form of Crank-Nicolson Finite Difference Time Domain

$H_{x_{i,j+1/2}}^{n+1}=H_{x_{i,j+1/2}}^n-{\mathrm{\χ}}_y\·{\mathrm{\Γ}}_y\left(E_{z_{i,j}}^n\right)+{\mathrm{\Δt}}_h\·(g_{{\mathrm{xy}}_{i,j+1/2}}^{n+1}+g_{{\mathrm{xy}}_{i,j+1/2}}^n)---\left(1\right)$

$H_{y_{i+1/2,j}}^{n+1}=H_{y_{i+1/2,j}}^n+{\mathrm{\χ}}_x\·{\mathrm{\Γ}}_x\left(E_{z_{i,j}}^n\right)-{\mathrm{\Δt}}_h\·(g_{{\mathrm{yx}}_{i+1/2,j}}^{n+1}+g_{{\mathrm{yx}}_{i+1/2,j}}^n)---\left(2\right)$

$\begin{array}{c}{E}_{z_{i,j}}^{n+1}=E_{z_{i,j}}^n+{\mathrm{\χ}}_x\·{\mathrm{\Γ}}_x\left(H_{y_{i+1/2,j}}^n\right)-{\mathrm{\χ}}_y\·{\mathrm{\Γ}}_y\left(H_{x_{i,j+1/2}}^n\right)\\ -{\mathrm{\Δt}}_h\·(f_{{\mathrm{zx}}_{i,j}}^{n+1}+f_{{\mathrm{zx}}_{i,j}}^n)+{\mathrm{\Δt}}_h\·(f_{{\mathrm{zy}}_{i,j}}^{n+1}+f_{{\mathrm{zy}}_{i,j}}^n)\end{array}---\left(3\right)$

In formula, Γ_η[*] is the shorthand in CN method, as

${\mathrm{\Γ}}_y\left(E_{z_{i,j}}^n\right)=E_{z_{i,j+1}}^{n+1}-E_{z_{i,j}}^{n+1}+E_{z_{i,j+1}}^n-E_{z_{i,j}}^n---\left(4\right)$

4. according to based on auxiliary differential equation the two-dimentional vacuum described in right 1 Crank-Nicolson completely permutation realizes algorithm, it is characterised in that: step 4, electric field is divided The iterative equation of amount arranges, can be to obtain the coefficient electricity that the left side is block tridiagonal matrix form Field implicit iterative equation

$\begin{array}{c}\[1-(D_{2x}+D_{2y})\]E_{z_{i,j}}^{n+1}=\[1+(D_{2x}+D_{2y})\]E_{z_{i,j}}^n\\ +2({\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ex}}_i})(H_{y_{i+1/2,j}}^n-H_{y_{i-1/2,j}}^n)-{\mathrm{\Δt}}_h(1+r_{0{\mathrm{ex}}_i})f_{{\mathrm{zx}}_{i,j}}^n\\ -2({\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ey}}_j})(H_{x_{i,j+1/2}}^n-H_{x_{i,j-1/2}}^n)+{\mathrm{\Δt}}_h(1+r_{0{\mathrm{ey}}_j})f_{{\mathrm{zy}}_{i,j}}^n\\ -({\mathrm{\χ}}_x-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ex}}_i}){\mathrm{\Δt}}_h\[(1+r_{0{\mathrm{hx}}_{i+1/2}})g_{{\mathrm{yx}}_{i+1/2,j}}^n-(1+r_{0{\mathrm{hx}}_{i-1/2}})g_{{\mathrm{yx}}_{i-1/2,j}}^n\]\\ -({\mathrm{\χ}}_y-{\mathrm{\Δt}}_h{r}_{1{\mathrm{ey}}_j}){\mathrm{\Δt}}_h\[(1+r_{0{\mathrm{hy}}_{j+1/2}})g_{{\mathrm{xy}}_{i,j+1/2}}^n-(1+r_{0{\mathrm{hy}}_{j-1/2}})g_{{\mathrm{xy}}_{i,j-1/2}}^n\]\end{array}---\left(5\right)$