-
有限体积法引入:
考虑如下守恒型方程:
\begin{equation}
\frac{\partial u}{\partial t} + \frac{\partial f}{\partial x} = 0
\label{eq_1}
\end{equation}
对第$i$个控制体进行积分:
\begin{equation}
\begin{split}
& \quad\int_{i - 1/2}^{i + 1/2} \frac{\partial u}{\partial t} \mathrm{d}x + \int_{i - 1/2}^{i + 1/2} \frac{\partial f}{\partial x} \mathrm{d}x = 0 \\
\Rightarrow & \quad\frac{\partial}{\partial t}\int_{i - 1 / 2}^{i + 1 / 2}u\mathrm{d}x + f_{i + 1/2} - f_{i - 1/2} = 0 \\
\Rightarrow & \quad\frac{\partial u_i}{\partial t} = \frac{1}{\Delta x}(f_{i - 1/2} - f_{i + 1/2})
\end{split}
\label{eq_2}
\end{equation}
左侧时间微分$\frac{\partial u_i}{\partial t}$采用Runge-Kutta方法处理之, 右侧单元界面处通量值$f_{i - 1/2}$、$f_{i + 1/2}$采用Godunov scheme处理之. -
Runge-Kutta方法:
考虑如下ODE:
\begin{equation}
\frac{\mathrm{d}u_i}{\mathrm{d}t} = g(t, u_i)
\label{eq_3}
\end{equation}
经典RK4给出如下数值解:
\begin{equation}
u_{i}^{n+1} = u_{i}^{n} + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)
\label{eq_4}
\end{equation}
其中:
\begin{gather*}
k_1 = g(t^n, u_i^n) \\
k_2 = g(t^n + \frac{h}{2}, u_i^n + \frac{h}{2}k_1) \\
k_3 = g(t^n + \frac{h}{2}, u_i^n + \frac{h}{2}k_2) \\
k_4 = g(t^n + h, u_i^n + hk_3)
\end{gather*} -
Godunov scheme:
Part Ⅰ: Godunov flux
purpose: upwind scheme
取$u_l$、$u_r$分别为单元界面左右两侧之u值, 则
\begin{equation}
f =
\begin{cases}
\min f(u), & \text{if } u_l < u_r \\
\max f(u), & \text{if } u_l > u_r
\end{cases}
\label{eq_5}
\end{equation}
其中, $u \in [\min\{u_l, u_r\}, \max\{u_l, u_r\}]$.
Part Ⅱ: reconstruction
purpose: ①. monotonicity; ②. high accuracy
考虑单元界面$i + 1/2$, 则
\begin{align}
u_{i + 1/2}^l = u_i + \frac{u_i - u_{i-1}}{2}\cdot \phi (r), & \quad\text{where } r = \frac{u_{i + 1} - u_i}{u_i - u_{i - 1}} \\
u_{i + 1/2}^r = u_{i+1} + \frac{u_{i+1} - u_{i+2}}{2}\cdot \phi (r), & \quad\text{where } r = \frac{u_i - u_{i+1}}{u_{i+1} - u_{i+2}}
\label{eq_6_7}
\end{align}
其中, $\phi (r)$为限制器函数, 取值范围如下:
\begin{equation}
\phi (r)
\begin{cases}
= 0, &\quad\text{if }r \leq 0 \\
\leq 2r, &\quad\text{if } 0 < r < 1 \\
= 1, &\quad\text{if }r = 1 \\
\leq 2, &\quad\text{if } r > 1
\end{cases}
\label{eq_8}
\end{equation}
常见2阶精度限制器函数如下:
\begin{equation}
\phi_{\text{superbee}}(r) =
\begin{cases}
0, &\quad\text{if }r \leq 0 \\
2r, &\quad\text{if }0 \leq r \leq 1/2 \\
1, &\quad\text{if }1/2 \leq r \leq 1 \\
r, &\quad\text{if }1 \leq r \leq 2 \\
2, &\quad\text{if } r \geq 2
\end{cases}
\label{eq_9}
\end{equation}
\begin{equation}
\phi_{\text{minmod}}(r) =
\begin{cases}
0, &\quad\text{if }r \leq 0 \\
r, &\quad\text{if }0 \leq r \leq 1 \\
1, &\quad\text{if }r \geq 1
\end{cases}
\label{eq_10}
\end{equation}
\begin{equation}
\phi_{\text{Van Leer}}(r) =
\begin{cases}
0, &\quad\text{if }r \leq 0 \\
\frac{2r}{1 + r}, &\quad\text{if }r \geq 0
\end{cases}
\label{eq_11}
\end{equation} -
一维Burgers方程初边值问题:
本文拟采用Van Leer限制器, 求解如下Burgers' equation:
\begin{equation*}
\frac{\partial u}{\partial t} + \frac{\partial u^2 / 2}{\partial x} = 0, \quad x \in \Omega = [0, 1], t \in [0, 0.5]
\end{equation*}
初边值条件如下:
\begin{equation*}
\begin{cases}
u(x, 0) = \sin(2\pi x), & \quad x \in \Omega \\
u(x, t) = 0, & \quad x \in \partial\Omega, t \in [0, 0.5]
\end{cases}
\end{equation*} -
代码实现:
View Code1 # Burgers' equation之求解 - 有限体积法 2 # 1. monotonicity 3 # 2. high accuracy 4 5 import numpy 6 from matplotlib import pyplot as plt 7 from matplotlib import animation 8 9 10 # Runge-Kutta方法 11 def RK4(func, t, u, h): 12 ''' 13 func: 导函数func(t, u) 14 t: 当前时刻 15 u: t时刻之u值 16 h: 下一时刻与当前时刻之时间差 17 ''' 18 half_h = h / 2 19 20 k1 = func(t, u) 21 k2 = func(t + half_h / 2, u + half_h * k1) 22 k3 = func(t + half_h / 2, u + half_h * k2) 23 k4 = func(t + h, u + h * k3) 24 25 uNew = u + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4) 26 return uNew 27 28 29 # Burgers' equation求解 30 class BurgersEq(object): 31 32 def __init__(self, nx=100, nt=100): 33 self.__nx = nx # x轴网格数 34 self.__nt = nt # t轴网格数 35 36 self.__init_grid() # 网格初始化 37 38 39 def get_solu(self): 40 UList = list() 41 42 U0 = self.__calc_U0() 43 UList.append(U0) 44 for t in self.__T[:-1]: 45 Uk = self.__calc_Uk(t, U0) 46 UList.append(Uk) 47 U0 = Uk 48 49 return self.__T, self.__X, UList 50 51 52 def __calc_Uk(self, t, UOld): 53 UNew = RK4(self.__calc_dUdT, t, UOld, self.__ht) 54 return UNew 55 56 57 def __calc_dUdT(self, t, UOld): 58 U = self.__fill_U_boundary(UOld) # u边界填充 59 UHL = self.__calc_UHalfLeft(U) # 不含首尾边界 60 UHR = self.__calc_UHalfRight(U) # 不含首尾边界 61 FH = self.__get_GodunovFlux(UHL, UHR) # 不含首尾边界 62 F = self.__fill_F_boundary(FH) # flux边界填充 63 dUdT = (F[:-1] - F[1:]) / self.__hx 64 return dUdT 65 66 67 def __fill_F_boundary(self, FH): 68 tmpFH = numpy.hstack(([0], FH, [0])) 69 return tmpFH 70 71 72 def __get_GodunovFlux(self, UHL, UHR): 73 FH = numpy.zeros(UHL.shape) 74 FHL = UHL ** 2 / 2 75 FHR = UHR ** 2 / 2 76 FH[UHL >= UHR] = numpy.max((FHL, FHR), axis=0)[UHL >= UHR] 77 FH[UHL < UHR] = numpy.min((FHL, FHR), axis=0)[UHL < UHR] 78 FH[(UHL < UHR) & (UHL * UHR < 0)] = 0 79 return FH 80 81 82 def __calc_UHalfRight(self, U): 83 Ui = U[1:-2] 84 UiPlusOne = U[2:-1] 85 UiPlusTwo = U[3:] 86 87 tmpItem = UiPlusOne - UiPlusTwo 88 UHR = list() 89 for idx, ele in enumerate(tmpItem): 90 if ele == 0: 91 UHR.append(UiPlusOne[idx]) 92 else: 93 r = (Ui[idx] - UiPlusOne[idx]) / ele 94 phi = self.__get_limitor(r) 95 UHR.append(UiPlusOne[idx] + ele / 2 * phi) 96 UHR = numpy.array(UHR) 97 return UHR 98 99 100 def __calc_UHalfLeft(self, U): 101 Ui = U[1:-2] 102 UiMinusOne = U[0:-3] 103 UiPlusOne = U[2:-1] 104 105 tmpItem = Ui - UiMinusOne 106 UHL = list() 107 for idx, ele in enumerate(tmpItem): 108 if ele == 0: 109 UHL.append(Ui[idx]) 110 else: 111 r = (UiPlusOne[idx] - Ui[idx]) / ele 112 phi = self.__get_limitor(r) 113 UHL.append(Ui[idx] + ele / 2 * phi) 114 UHL = numpy.array(UHL) 115 return UHL 116 117 118 def __get_limitor(self, r): 119 phi = 0 120 if r >= 0: 121 phi = 2 * r / (1 + r) 122 return phi 123 124 125 def __fill_U_boundary(self, UOld): 126 tmpUOld = numpy.hstack(([0], UOld, [0])) 127 return tmpUOld 128 129 130 def __calc_U0(self): 131 ''' 132 计算第0层时间步数值 133 ''' 134 U0 = self.__calc_u0(self.__X) 135 return U0 136 137 138 def __calc_u0(self, x): 139 u0 = numpy.sin(2 * numpy.pi * x) 140 return u0 141 142 143 def __init_grid(self): 144 xMin, xMax = 0, 1 145 tMin, tMax = 0, 0.5 146 self.__hx = (xMax - xMin) / self.__nx 147 self.__ht = (tMax - tMin) / self.__nt 148 self.__X_interface = numpy.linspace(xMin, xMax, self.__nx + 1) 149 self.__X = (self.__X_interface[:-1] + self.__X_interface[1:]) / 2 150 self.__T = numpy.linspace(tMin, tMax, self.__nt + 1) 151 152 153 class BurgersPlot(object): 154 155 fig = None 156 ax1 = None 157 line = None 158 txt = None 159 T, X, UList = None, None, None 160 161 @classmethod 162 def plot_ani(cls, burgersObj): 163 cls.fig = plt.figure(figsize=(6, 4)) 164 cls.ax1 = plt.subplot() 165 cls.line = cls.ax1.plot([], [], lw=3)[0] 166 cls.txt = cls.ax1.text(x=0.8, y=0.8, s="t = {:.3f}".format(0)) 167 cls.T, cls.X, cls.UList = burgersObj.get_solu() 168 169 ani = animation.FuncAnimation(cls.fig, cls.update, frames=numpy.arange(len(cls.UList)), init_func=cls.init, blit=True, interval=50, repeat=True) 170 ani.save("plot_ani.gif", writer="PillowWriter", fps=20) 171 172 plt.close() 173 174 175 @classmethod 176 def update(cls, frame): 177 cls.line.set_data(cls.X, cls.UList[frame]) 178 cls.txt.set_text("t = {:.3f}".format(cls.T[frame])) 179 return cls.line, cls.txt 180 181 182 @classmethod 183 def init(cls): 184 cls.ax1.set(xlim=(0, 1), ylim=(-1, 1), xlabel="x", ylabel="u") 185 return cls.line, cls.txt 186 187 188 @staticmethod 189 def plot_fig(burgersObj): 190 T, X, UList = burgersObj.get_solu() 191 192 fig = plt.figure(figsize=(6, 4)) 193 ax1 = plt.subplot() 194 for U in UList: 195 ax1.plot(X, U) 196 ax1.set(xlabel="x", ylabel="u") 197 198 fig.savefig("plot_fig.png", dpi=100) 199 # plt.show() 200 plt.close() 201 202 203 204 if __name__ == "__main__": 205 burgersObj = BurgersEq(200, 100) 206 BurgersPlot.plot_ani(burgersObj) 207 # BurgersPlot.plot_fig(burgersObj)
相关文章：