Discontinuous Galerkin method for steady transport problem
2016-04-04 22:45
225 查看
下面讨论如何使用 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)\)。
1.简介
恒定状态对流方程\[\begin{equation}
a\cdot \nabla \mathbf{u} = f
\end{equation}\]
出现在多种问题中,如海洋模型中求解连续方程计算垂向速度,明渠恒定流动问题等。
2.数值离散
首先需要将方程写为离散格式,以一维问题为例\[\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\]
即始终采用边界左侧节点值计算数值通量。这主要是因为边界条件定义在计算域左侧,在计算时也将由左向右逐个单元进行计算,也就是数据从左向右进行传递。
3.联立求解
与包含时间项的非恒定方程不同,恒定对流方程需要将系数矩阵联立求解。与普通系数矩阵构造不同,这里首先设多元函数 \(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)\)。
4.边界条件
在恒定输运方程中,在给定边界 \(\Gamma_D\) 上为Dirichlet边界 \(u = u_D\)。参考有限元方法,可以采用置大数法,即修改对应系数矩阵 \(A\),与源项 \(f\) 来耦合 \(\Gamma_D\) 上已知解。具体方法请参考有限元边界 Dirichlet 条件处理
5.代码
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
6.计算结果
已知方程精确解为 \(u(x) = -cos(x) + 2\),分别采用不同阶(N=1,2,3)与不同个数(Ne=10,20,40,80)单元进行计算,统计对应的 \(L_1\)、\(L_2\) 误差及收敛速率。6.1.N=1
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 |
6.2.N=2
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 |
6.3.N=3
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 |
相关文章推荐
- 两个Activity之间的观察者模式
- 2016.4.4HTML表格,列表,块,布局
- cocospod 搭建环境和使用
- 分析Linux内核启动过程:从start_kernel到init
- Gson解析数据使用笔记
- bzoj 1812 river
- Alphabet 和 Facebook的Stratospheric Internet纠结于繁文缛节
- Android消息处理机制
- grafana图形化数据展示软件介绍
- Socket使用大全
- MQTT协议学习笔记
- shiro整合redis做缓存
- storm安装
- strcpy 字符串拷贝函数解析
- 二叉树遍历
- SignalR技术
- 借助 Lucene.Net 构建站内搜索引擎(下)
- 20159315《网络攻防实践》第五周学习总结
- 多线程并发测试代码及解决方法
- Content Delivery Network (CDN)