下面讨论如何使用 Discontinuous Galerkin 求解恒定对流问题。
恒定状态对流方程
\[\begin{equation}
a\cdot \nabla \mathbf{u} = f
\end{equation}\]
出现在多种问题中,如海洋模型中求解连续方程计算垂向速度,明渠恒定流动问题等。
首先需要将方程写为离散格式,以一维问题为例
\[\begin{eqnarray}
\begin{aligned}
& \frac{\partial u}{\partial x} = sin(x) \quad x\in[0, 2\pi] \cr
& u(0) = 1
\end{aligned}
\end{eqnarray}\]
将方程乘以实验函数(test function)并在控制单元内积分,利用分部积分将方程转化为
\[\begin{equation}
\oint_{\Omega}(l_i \cdot u) n_x dx - \int_{\Omega} \frac{\partial l_i}{\partial x} u dx = \int_{\Omega} l_i f_h dx
\end{equation}\]
用基函数线性组合近似变量 \(u \approx u_h = \sum_{j=1}^P l_j u_j\),源项 \(f_h\) 也采用基函数线性近似 \(f_h = \sum_{j=1}^P l_j f_j\) 结果为
\[\begin{equation}
\left( \oint_{\Omega}(l_i \cdot l_j) dx \right) u_j \cdot n_x - \left( \int_{\Omega} \frac{\partial l_i}{\partial x} l_j dx \right) u_j = \left( \int_{\Omega} l_i l_j dx \right) f_j
\end{equation}\]
写成矩阵形式为
\[\begin{equation}
J_s M_e \hat{u}_n - J M Dr \cdot u = J M \cdot f
\end{equation}\]
其中边界通量 \(\hat{u}_n\) 采用如下方式计算
\[\hat{u}_n = u_L \cdot n_x
\]
即始终采用边界左侧节点值计算数值通量。这主要是因为边界条件定义在计算域左侧,在计算时也将由左向右逐个单元进行计算,也就是数据从左向右进行传递。
与包含时间项的非恒定方程不同,恒定对流方程需要将系数矩阵联立求解。
与普通系数矩阵构造不同,这里首先设多元函数 \(L:\mathbb{R}^{Np} \to \mathbb{R}^{Np}\),
\[L(\mathbf{u}) = J_s M_e \hat{u}_n - J M Dr \cdot \mathbf{u} = (L_1, L_2, \cdots, L_{Np})^T
\]
最终目标是寻找变量 \(\mathbf{u}_0\),使得等式 \(L(\mathbf{u}_0) = JM \cdot f\) 成立。
这里令 \(\mathbf{e}_i\) 表示第\(i\)个分量为单位1,其余分量为0的解。假设函数满足线性关系:若
\[\begin{equation}
u_0 = \sum_{i=1}^{Np} u_i \mathbf{e}_i
\end{equation}\]
那么对应的函数有
\[\begin{equation}
L(\mathbf{u}_0) = \sum_{i=1}^{Np} u_i L(\mathbf{e}_i)
\end{equation}\]
那么我们通过构造系数矩阵 \(A = \left( L(\mathbf{e}_1), L(\mathbf{e}_2), \cdots, L(\mathbf{e}_{Np}) \right)\),联立 \(Au_i = JM \cdot f\),便可得到最终未知解 \(u_0 = \left( u_1, u_2, \cdots, u_{Np} \right)\)。
在恒定输运方程中,在给定边界 \(\Gamma_D\) 上为Dirichlet边界 \(u = u_D\)。参考有限元方法,可以采用置大数法
,即修改对应系数矩阵 \(A\),与源项 \(f\) 来耦合 \(\Gamma_D\) 上已知解。具体方法请参考有限元边界 Dirichlet 条件处理
SteadyConvectionDriver
负责构造计算所需的网格及标准线单元系数矩阵(刚度矩阵,质量矩阵等),方程源项 \(f\) 也在此给定。
function SteadyConvectionDriver
% solving steady convection problem by DGM
% \nabla u = sin(x)
%
x1 = 0; x2 = 2*pi; % domain
N = 1; nElement = 20;
[~, VX, ~, EToV] = Utilities.Mesh.MeshGen1D(x1, x2, nElement);
BC = [2,1];
%
line = StdRegions.Line(N);
mesh = MultiRegions.RegionLineBC(line, EToV, VX, BC);
f = sin(mesh.x);
u = SteadyConvectionSolver(mesh, f);
plot(mesh.x(:), u(:), 'b', mesh.x(:), cos(mesh.x(:)), 'r');
end% func
SteadyConvectionSolver
负责求解方程组,根据所给边界条件,采用置大数法修改对应系数矩阵及右端项系数,
function u = SteadyConvectionSolver(mesh, f)
% set up and solve the equation system
% Input:
% mesh - mesh object
% f - source term
% Output:
% u - unknown variable
f = mesh.J.*(mesh.Shape.M*f);
%% set up and solve global matrix coeffcient
% get system global matrix coefficient
A = SteadyConvectionCoeffMatrix(mesh);
% boundary condition
u0 = 1; M = 1e8;
A(1, 1) = M; f(1) = u0*M;
solvec = A\f(:);
u = reshape(solvec, size(mesh.x) );
end% func
SteadyConvectionCoeffMatrix
负责构造系数矩阵,每次计算单位变量 \(\mathbf{e}_{i}\) 对应的多元函数值 \(L(\mathbf{e}_{i})\)
function A = SteadyConvectionCoeffMatrix(mesh)
% set up symmetric matrix
A = zeros(mesh.nNode, mesh.nNode);
g = zeros(size(mesh.x));
% Build matrix -- one column at a time
for i = 1:mesh.nNode
g(i) = 1;
Avec = SteadyConvectionRHS(mesh, g);
A(:, i) = Avec(:);
g(i) = 0;
end% for
end% func
SteadyConvectionRHS
计算函数 \(L(\mathbf{u})\),根据信息传递方向,在边界处数值通量采用迎风格式计算。
function rhs = SteadyConvectionRHS(mesh, u)
% right hands of equation
%
us = zeros( size(u(mesh.vmapM)) );
us(1, :) = u(mesh.vmapP(1, :));
us(2, :) = u(mesh.vmapM(2, :));
us = us.*mesh.nx;
rhs = (mesh.Shape.Mef * us) - mesh.J.*( mesh.rx .*( mesh.Shape.Dr'*(mesh.Shape.M*u) ));
end% func
已知方程精确解为 \(u(x) = -cos(x) + 2\),分别采用不同阶(N=1,2,3)与不同个数(Ne=10,20,40,80)单元进行计算,统计对应的 \(L_1\)、\(L_2\) 误差及收敛速率。
Ne
L1
rate
10
0.029536
\
20
0.007977
1.888594
40
0.002040
1.967584
80
0.000513
1.991309
Ne
L2
rate
10
0.037214
\
20
0.009872
1.914490
40
0.002506
1.978170
80
0.000629
1.994526
Ne
L1
rate
10
0.001776
\
20
0.000172
3.368547
40
0.000019
3.170322
80
0.000002
3.081502
Ne
L2
rate
10
0.002184
\
20
0.000233
3.227791
40
0.000028
3.073642
80
0.000003
3.020009
Ne
L1
rate
10
0.000175
\
20
0.000011
3.938733
40
0.000001
3.975718
80
0.000000
3.855050
Ne
L2
rate
10
0.000189
\
20
0.000012
3.946121
40
0.000001
3.978587
80
0.000000
3.871224
手机扫一扫
移动阅读更方便
你可能感兴趣的文章