Application of C-Bézier and H-Bézier Basis in Solving Heat Conduction Problems
-
摘要:
运用C-Bézier和H-Bézier基函数构造有限元方法中的试探函数和测试函数空间, 并分析了误差。通过数值算例, 数值解的误差精度较Lagrange基函数提升了1到3个数量级, 说明这两类基函数在模拟特定类型的热传导问题时具有更好的逼近效果。
-
关键词:
- 热传导方程 /
- 有限元方法 /
- C-Bézier基函数 /
- H-Bézier基函数
Abstract:C-Bézier and H-Bézier basis functions are applied to the space of trial function and test function in the finite element method and the errors are analyzed. By means of numerical examples, the error accuracy of numerical solution is improved by one to three orders of magnitude compared with Lagrange basis function, it is shown that these two kinds of basis functions have better approximations in simulating specific heat conduction problems.
-
开放科学(资源服务)标识码(OSID):
0. 引言
在自然界中, 热传导、电磁势和流体等问题都可用偏微分方程来描述, 它们的解一般很难用解析公式表达出来, 运用数值方法来计算它们的近似解是十分重要的。有限元方法是现代科学与工程计算领域中求解偏微分方程的一个主要工具。1943年, 数学家COURANT[1]提出了有限元方法;20世纪50年代, 有限元方法是作为解决固体力学问题的方法出现的;1960年, CLOUGH[2]在处理平面弹性问题时, 第一次提出了“有限元方法”这一名称, 并成功应用于飞机结构的分析;20世纪60年代, 冯康发表的论文《基于变分原理的差分格式》[3], 标志着有限元方法在我国的诞生。基于传统的有限元方法, 学者们也提出了许多新的数值计算方法, 例如: 混合有限元方法[4]、自适应有限元法[5]、DG方法[6]、弱Galerkin有限元法[7]。
有限元基函数的选取会对数值解的精确度产生影响。为解决高阶多项式曲线的表示问题, 文献[8-9]等对混合代数和三角多项式空间进行了扩展, 在空间T=span{1, t, …, tp-2, sin t, cos t}上定义了p次C-Bézier基函数以及空间T=span{sinh t, cosh t, tq-2, …, t, 1}上的q次H-Bézier基函数。
如何得到一个精确度高且计算量较小的方法一直是有限元方法研究中的一个重要课题。HUGES等[10]将NURBS基函数应用于有限元分析,BHATTI等[11]提出在求解偏微分方程中使用Bernstein基函数。本文考虑一维热传导方程初边值问题的数学模型, 在有限元方法中利用C-Bézier和H-Bézier基函数构造试探函数和测试函数空间, 通过数值实验来验证理论结果的正确性与方法的可行性。
1. 一维热传导方程
热传导问题是描述在某个区域内物体的温度分布规律。给定一根均匀且各向同性的细杆, 长为L(L≥0), 假设细杆内部有热源且与周边介质有热交换, f0(x, t)为热源强度, C为热容, ρ为密度, k为传导率, u=u(x, t)描述在t时刻x位置材料的温度, 则Dirichlet边界条件下的一维热传导方程初边值问题为[12-13]:
{ut−∇⋅(c∇u)=f,x∈Ω,t>0,u(x,0)=ψ(x),x∈Ω,t=0,u(0,t)=ga,u(L,t)=gb, (1) 其中: Ω=(0,L),c=kρC为热扩散项系数, f=f0(x, t)/C, 初始温度为ψ(x), 边界条件为ga和gb。
一维热传导方程初边值问题的变分形式为: 求u=u(x, t)∈H1(Ω)(t为参数), 满足
{(ut,v)+a(u,v)=(f,v),u(x,0)=ψ(x), (2) 对任意的v∈H01(Ω)都成立。
在H1(Ω)中取一n维子空间Vh, 令Vh0表示Vh中在边界处值为0的函数构成的空间, 则半离散的Gakerkin方法为:求uh(x, t)∈Vh, 满足
{(uht,vh)+a(uh,vh)=(f,vh),∀vh∈Vh0,uh(0)=ψh。 (3) 进一步利用差分法对时间t离散, 得到一维热传导方程的全离散化Galerkin方法。
2. C-Bézier和H-Bézier基函数
定义1 在空间T=span{1, ω, …, ωp-2, sin ω, cos ω}上的p次C-Bézier基函数{Cip(ω)}i=0p为:
Cp0(ω)=1−∫ω0δp−10Cp−10(s)ds,Cpi(ω)=∫ω0δp−1i−1Cp−1i−1(s)ds−∫ω0δp−1iCp−1i(s)ds,(i=1,2,⋯,p−1),Cpp(ω)=∫ω0δp−1p−1Cp−1p−1(s)ds, (4) 其中: p≥2, ω∈[0, α], α∈(0, π],
C10(ω)=sin(α−ω)sinα,C11(ω)=sinωsinα,δpi=(∫α0Cpi(ω)dω)−1, α为形状参数[8]。
定义2 在空间T=span{1, ω, …, ωq-2, sinh ω, cosh ω}上的q次H-Bézier基函数{Hiq(ω)}i=0q为:
Hq0(ω)=1−∫ω0δq−1oHq−10(s)ds,Hqi(ω)=∫ω0δq−1i−1Hq−1i−1(s)ds−∫ω0δq−1iHq−1i(s)ds,(i=1,2,⋯,q−1),Hqq(ω)=∫ω0δq−1q−1Hq−1q−1(s)ds, (5) 其中: q≥2, ω∈[0, α], α∈(0, +∞),
H10(ω)=sinh(α−ω)sinhα,H11(ω)=sinhωsinhα,δqi=(∫α0Hqi(ω)dω)−1, α为形状参数[9]。
C-Bézier和H-Bézier基函数都具有单位分解、端点插值、升阶等性质, 具体性质及证明见参考文献[8-9]。
定义3 称参数曲线段
P(ω)=p∑i=0PiCpi(ω),ω∈[0,α], (6) 为一条p次C-Bézier曲线, 其中α∈(0, π], Cip(ω)(i=0, 1, …, p)为p次C-Bézier基函数, Pi(i=0, 1, …, p)为控制顶点, α称为全局形状参数[8]。
图 1(a)为α=π4时的二次C-Bézier曲线, 控制顶点为P0=(0, 2), P1=(0, 0), P2=(2, 0), 图 1(b)为α=π时的三次的C-Bézier曲线, 对应的控制顶点为P0=(0, 0), P1=(2, 4), P2=(4, 4), P3=(6, 0)。
定义4 称参数曲线段
Q(ω)=q∑i=0QiHqi(ω),ω∈[0,α], (7) 为一条q次H-Bézier曲线, 其中α∈(0, ∞), Hiq(ω)(i=0, 1, …, q)为q次H-Bézier基函数, Qi(i=0, 1, …, q)为控制顶点, α称为全局形状参数[9]。
图 2(a)为α=2时的二次H-Bézier曲线, 控制顶点为Q0=(4, 0), Q1=(0, 2), Q2=(4, 4), 图 2(b)为α=3时的三次H-Bézier曲线, 其控制顶点为Q0=(0,0),Q1=(32,5),Q2=(3,6),Q3=(5,2)。
C-Bézier和H-Bézier基函数引入了形状参数α, α可以调节曲线的形状和几何性质。以C-Bézier基函数为例, 调节区域为Bézier曲线到α=π时的C-Bézier曲线之间。曲线首端和末端的位置及切矢方向不随形状参数α的改变而改变。同时在多段C-Bézier曲线进行G1拼接时, α改变不会影响曲线之间原有的连续性, 为曲线曲面造型带来了便利性和灵活性。
3. 以C-Bézier和H-Bézier为基函数的有限元方法
将区间[0, L]均匀剖分为N个单元, 步长h=LN, 设节点为xi=(i-1)h(i=1, 2, …, N+1), En=[xn, xn+1](n=1, …, N)表示网格单元, Nb表示有限元节点个数。令
Vh=span{φj}Nbj=1∈H1(Ω) 为C-Bézier或H-Bézier基函数构成的测试函数空间。
当p=2时, 二次C-Bézier基函数为:
φ1=1−cos(α−ω)1−cosα,φ2=−1−cosω+cosα−cos(α−ω)1−cosα,φ3=1−cosω1−cosα, (8) 其中ω∈[0, α], α∈(0, π]。
当q=2时, 二次H-Bézier基函数为:
φ1=1−cosh(α−ω)1−coshα,φ2=−1−coshω+coshα−cosh(α−ω)1−coshα,φ3=1−coshω1−coshα, (9) 其中: ω∈[0, α], α∈(0, +∞)。图 3(a)为α=π6、p=2时的C-Bézier基函数, 图 3(b)为α=π、q=2时的H-Bézier基函数。
做变量替换m=ωα, 则m∈[0, 1], ω=mα, 通过仿射变换, 可得任一单元[xn, xn+1]上的二次C-Bézier和H-Bézier基函数。以C-Bézier基函数为例:
φn1=1−cos(xn+1−xhα)1−cosα,φn2=−1−cos(x−xnhα)+cosα−cos(xn+1−xhα)1−cosα,φn3=1−cos(x−xnhα)1−cosα, (10) 其中x∈[xn, xn+1], α∈(0, π]。
设uh(x,t)=Nb∑j=1uj(t)φj, 并取vh=φi(i=1, …, Nb), 得到关于u1(t), …, uNb(t)的常微分方程组:
Nb∑j=1u′j(t)∫L0(φjφi)dx+Nb∑j=1uj(t)∫L0(c∇φj⋅∇φi)dx=∫L0fφi dx,(i=1,⋯,Nb)。 (11) 用向量表示为:
\boldsymbol{M} \boldsymbol{X}^{\prime}(t)+\boldsymbol{A}(t) \boldsymbol{X}(t)=\boldsymbol{b}(t), (12) 进一步用差分法对时间t离散, 得到全离散的Galerkin方法。
4. 稳定性及误差分析
全离散格式是由常微分方程组(11)离散化得到的, 因此需要讨论它的稳定性。使用Crank-Nicolson格式对时间t进行离散:
\begin{aligned} & \frac{1}{\tau}\left(u_{h}^{n+1}-u_{h}^{n}, v_{h}\right)+\frac{1}{2} a\left(u_{h}^{n+1}+u_{h}^{n}, v_{h}\right)= \\ & \quad\left(f, v_{h}\right), \end{aligned} 其中: 上标n表示在t=tn=nτ处的近似。取
v_{h}=u_{h}^{n+1}-u_{h}^{n}, f=0, 得
\begin{aligned} & \frac{1}{\tau}\left\|u_{h}^{n+1}-u_{h}^{n}\right\|^{2}+ \\ & \quad \frac{1}{2} a\left(u_{h}^{n+1}+u_{h}^{n}, u_{h}^{n+1}-u_{h}^{n}\right)=0。 \end{aligned} 利用a(·, ·)的对称性得
\left|u_{h}^{n+1}\right|_{1} \leqslant\left|u_{h}^{n}\right|_{1} \leqslant \cdots \leqslant\left|u_{h}^{0}\right|_{1}, 由Poincaré不等式可以证得稳定性。
由Céa引理[14]可得, 有限元的误差估计转化为插值误差估计。当一维热传导方程半离散时, 即仅对空间离散, L2范数误差估计为
\left\|u-u_{h}\right\| \leqslant C_{1} h^{2}\|u\|_{2}, 其中:C1为常数, h为步长, 利用Aubin-Nitche技巧可证明之。
引理1 设u和uh分别为方程(1)和(3)的解, 对任意的t∈[0, +∞), 误差估计为
\begin{aligned} & \left\|u_{h}(t)-u(t)\right\| \leqslant \\ & \quad C_{2} h^{2}\left(\|\psi\|_{2}+\int_{0}^{t}\left\|u_{t}\right\|_{2} \mathrm{~d} x\right), \end{aligned} 其中:C2为常数。
引理1的证明参见文献[14], 对时间t进行离散, 由引理1可以得到一维热传导方程全离散化的误差估计。
5. 数值算例
运用Lagrange、C-Bézier和H-Bézier基函数的有限元法求解热传导方程并比较数值解的精度。在下述算例中, 均使用二次的基函数构成有限元方法中的试探函数与测试函数空间, 有限元节点相同, 所用网格一致。
例1 考虑一根均匀且长度为1 m的细杆, 初始温度分布为0 ℃, 热扩散系数c=2, 热容为1, 细杆内部有热源, 在t时刻, 热源强度f0(x, t)=xcos(xt)+2t2sin(xt), 细杆左端保持温度0 ℃, 右端温度为sin t ℃, 运用有限元方法求解细杆内部温度分布, 并分析误差。
其一维热传导方程的初边值问题为:
\left\{\begin{array}{l} u_{t}-\nabla \cdot(2 \nabla u)=x \cos (x t)+2 t^{2} \sin (x t), \\ \quad 0<x<1, t>0, \\ u(x, 0)=0, 0<x<1, t=0, \\ u(0, t)=0, u(1, t)=\sin t, t>0。 \end{array}\right. (13) 方程(13)的精确解为u=sin(xt), 运用有限元方法, 在t=2时刻得到在节点处的数值解误差见表 1, 图 4(a)为步长h=\frac{1}{4}, t=2时的误差曲线图, 细杆内温度随时间的变化情况如图 4(b)所示。
表 1 方程(13)在t=2时刻节点处的无穷范误差Table 1. Infinite norm errors at nodes of Eq.(13) for t=2Lagrange基函数 C-Bézier基函数 h ‖u-uh‖∞ α ‖u-uh‖∞ 1/4 2.829 4e-05 π/6 5.134 9e-06 1/8 1.800 4e-06 π/12 3.221 3e-07 1/16 1.136 9e-07 π/24 2.016 9e-08 1/32 7.159 2e-08 π/48 1.261 1e-09 由例1可见, 运用C-Bézier基函数构成测试函数及试探函数空间, 在步长不变的情况下, 通过选取适当的形状参数α, 此时\alpha=\frac{2}{3} \pi h, 得到的数值解的精度比Lagrange型基函数高一个数量级, 说明C-Bézier基函数的有限元法对热传导过程的数值模拟效果更好。
例2 考虑一根均匀且长为1 m的细杆, 初始温度分布为0 ℃, 热扩散系数c和热容都为1。细杆内部有热源, 在t时刻, 热源强度f0(x, t)=\cosh \left(\frac{x}{2}\right)\left(2-\frac{t}{2}\right), 细杆左端温度是2t ℃, 右端为2 t \cosh \left(\frac{1}{2}\right) ℃, 运用有限元方法求解细杆内部温度分布, 并分析误差。
其一维热传导方程边值问题为
\left\{\begin{array}{l} u_{t}-\nabla \cdot(\nabla u)=\cosh \left(\frac{x}{2}\right)\left(2-\frac{t}{2}\right), \\ \quad 0<x<1, t>0, \\ u(x, 0)=0, 0<x<1, t=0, \\ u(0, t)=2 t, u(1, t)=2 t \cosh \left(\frac{1}{2}\right), t>0。 \end{array}\right. (14) 方程(14)的精确解为u=2 t \cosh \left(\frac{x}{2}\right), 运用有限元方法, 在t=1时刻得到在节点处的数值解误差见表 2, 图 5(a)为h=\frac{1}{4}, t=1时的误差图像, 细杆内温度的分布情况如图 5(b)所示。
表 2 方程(14)在t=1时刻在节点处的无穷范误差Table 2. Infinite norm errors at nodes of Eq.(14) for t=1Lagrange基函数 H-Bézier基函数 h ‖u-uh‖∞ α ‖u-uh‖∞ 1/4 2.445 3e-07 1/8 2.744 5e-13 1/8 1.643 9e-08 1/16 6.483 0e-14 1/16 1.071 0e-09 1/32 1.465 5e-14 1/32 6.841 1e-11 1/64 8.304 5e-14 由例2可见, 运用H-Bézier基函数构成测试函数及试探函数空间, 在步长不变的情况下, 通过选取适当的形状参数α, 此时\alpha=\frac{1}{2} h, 较传统Lagrange基函数,当h取1/4或1/8时,精度可提高6个数量级;当h取1/16时,精度可提高5个数量级;当h取1/32时,精度可提高3个数量级;综上,得到的数值解的精度比Lagrange型基函数高3个数量级以上, 说明H-Bézier基函数的有限元法对精确解的逼近效果更好。
例3 考虑一根均匀且长为1 m的细杆, 初始温度分布为0 ℃, 热扩散系数c和热容都为1。细杆内部有热源, 在t时刻, 热源强度f0(x, t)=tsin(πx)(2+π2t), 细杆两端温度为0 ℃。运用有限元方法求解细杆内部温度分布, 并分析误差。其一维热传导方程的初边值问题为:
\left\{\begin{array}{l} u_{t}-\nabla \cdot(\nabla u)=t \sin (\pi x)\left(2+\pi^{2} t\right), \\ \quad 0<x<1, t>0, \\ u(x, 0)=0, 0<x<1, t=0, \\ u(0, t)=0, u(1, t)=0, t>0。 \end{array}\right. (15) 方程(15)的精确解为u=t2sin(πx), 运用有限元方法, 在t=1时刻得到在节点处的数值解误差见表 3, 图 6(a)为h=\frac{1}{4},t=1时的误差图像, 细杆内温度的分布情况如图 6(b)所示。
表 3 方程(15)在t=1时刻在节点处的无穷范误差Table 3. Infinite norm errors at nodes of Eq.(15) for t=1Lagrange基函数 C-Bézier基函数 h ‖u-uh‖∞ α ‖u-uh‖∞ 1/4 9.818 5e-05 π/4 2.147 0e-08 1/8 6.312 6e-06 π/8 8.104 0e-11 1/16 3.972 9e-07 π/16 1.252 3e-13 1/32 2.487 3e-08 π/32 9.059 4e-14 由例3可见, 运用C-Bézier基函数构成测试函数及试探函数空间, 在步长不变的情况下, 通过选取适当的形状参数α, 此时α=πh, 得到的数值解的精度比Lagrange型基函数高三个数量级以上, 说明C-Bézier基函数的有限元法对热传导过程的模拟效果更好。
6. 总结
介绍了C-Bézier和H-Bézier基函数的定义及良好性质, 并用其构造有限元方法中的测试函数与试探函数空间, 通过例子说明利用这两种基函数求出的数值解与精确解的误差更小, 精度更高, 证明了理论结果。
形状参数α的选取会影响数值解的精确度, 同时在进行误差估计时, 常数C1、C2的确定也非常重要。下一步将继续研究如何选取最优的形状参数α, 使得数值解与精确解之间的误差达到最小, 并探讨如何利用简单的方法确定常数C1、C2进行误差估计。
-
表 1 方程(13)在t=2时刻节点处的无穷范误差
Table 1 Infinite norm errors at nodes of Eq.(13) for t=2
Lagrange基函数 C-Bézier基函数 h ‖u-uh‖∞ α ‖u-uh‖∞ 1/4 2.829 4e-05 π/6 5.134 9e-06 1/8 1.800 4e-06 π/12 3.221 3e-07 1/16 1.136 9e-07 π/24 2.016 9e-08 1/32 7.159 2e-08 π/48 1.261 1e-09 表 2 方程(14)在t=1时刻在节点处的无穷范误差
Table 2 Infinite norm errors at nodes of Eq.(14) for t=1
Lagrange基函数 H-Bézier基函数 h ‖u-uh‖∞ α ‖u-uh‖∞ 1/4 2.445 3e-07 1/8 2.744 5e-13 1/8 1.643 9e-08 1/16 6.483 0e-14 1/16 1.071 0e-09 1/32 1.465 5e-14 1/32 6.841 1e-11 1/64 8.304 5e-14 表 3 方程(15)在t=1时刻在节点处的无穷范误差
Table 3 Infinite norm errors at nodes of Eq.(15) for t=1
Lagrange基函数 C-Bézier基函数 h ‖u-uh‖∞ α ‖u-uh‖∞ 1/4 9.818 5e-05 π/4 2.147 0e-08 1/8 6.312 6e-06 π/8 8.104 0e-11 1/16 3.972 9e-07 π/16 1.252 3e-13 1/32 2.487 3e-08 π/32 9.059 4e-14 -
[1] COURANT R. Variational methods for the solution of problems of equilibrium and vibrations[J]. Bulletin of the American Mathematical Society, 1943, 49: 1-23.
[2] CLOUGH R. The finite element method in plane stress analysis[C]//Proceedings of the 2nd Conference on Electronic Computation of American Society of Civil Engineers, 1960: 345-378.
[3] 冯康. 基于变分原理的差分格式[J]. 应用数学与计算数学, 1965, 2(4): 238-262. https://www.cnki.com.cn/Article/CJFDTOTAL-JSWL200204003.htm FENG Kang. Difference scheme based on variational principle[J]. Applied Mathematics & Computational Mathematics, 1965, 2(4): 238-262. https://www.cnki.com.cn/Article/CJFDTOTAL-JSWL200204003.htm
[4] BABUŠKA I. Error-bounds for finite element method[J]. Numerische Mathematik, 1971, 16(4): 322-333. doi: 10.1007/BF02165003
[5] LI F, YI N Y. Analysis of a goal-oriented adaptive two-grid finite-element algorithm for semilinear elliptic problems[J]. Computational and Applied Mathematics, 2022, 41(3): 108. doi: 10.1007/s40314-022-01815-4
[6] TASSI P A, RHEBERGEN S, VIONNET C A. A discontinuous Galerkin finite element model for river bed evolution under shallow flows[J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(33/40): 2930-2947.
[7] AL-TAWEEL A, WANG X S, YE X, et al. A stabilizer free weak Galerkin finite element method with super closeness of order two[J]. Numerical Methods for Partial Differential Equations, 2021, 37(2): 1012-1029. doi: 10.1002/num.22564
[8] CHEN Q Y, WANG G Z. A class of Bézier-like curves[J]. Computer Aided Geometric Design, 2003, 20(1): 29-39. doi: 10.1016/S0167-8396(03)00003-7
[9] LYU Y G, WANG G Z, YANG X N. Uniform hyperbolic polynomial B-spline curves[J]. Computer Aided Geometric Design, 2002, 19(6): 379-393. doi: 10.1016/S0167-8396(02)00092-4
[10] HUGHES T J R, COTTRELL J A, BAZILEVS Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39/40/41): 4135-4195.
[11] BHATTI M I, BRACKEN P. Solutions of differential equations in a Bernstein polynomial basis[J]. Journal of Computational and Applied Mathematics, 2007, 205(1): 272-280. doi: 10.1016/j.cam.2006.05.002
[12] MEZGHANI B, TOUNSI F, MASMOUDI M. Development of an accurate heat conduction model for micromachined convective accelerometers[J]. Microsystem Technologies, 2015, 21(2): 345-353. doi: 10.1007/s00542-014-2079-x
[13] LIU Q, MING P J, ZHAO H Y, et al. A high order control volume finite element method for transient heat conduction analysis of multilayer functionally graded materials with mixed grids[J]. Journal of Thermal Science, 2020, 29(1): 144-158. doi: 10.1007/s11630-019-1167-8
[14] THOMÉE V. Galerkin finite element methods for parabolic problems[M]. Heidelberg: Springer-Verlag, 1984: 186-187.
[15] BADIA S, VERDUGO F, MARTÍN A F. The aggregated unfitted finite element method for elliptic problems[J]. Computer Methods in Applied Mechanics and Engineering, 2018, 336: 533-553. doi: 10.1016/j.cma.2018.03.022
[16] BAUMANN C E, ODEN J T. A discontinuous hp finite element method for convection-diffusion problems[J]. Computer Methods in Applied Mechanics and Engineering, 1999, 175(3/4): 311-341.
[17] CHEN Zhiming, TUO Rui, ZHANG Wenlong. A finite element method for elliptic problems with observational boundary data[J]. Computational Mathematics, 2020, 38(2): 355-374. doi: 10.4208/jcm.1810-m2017-0168
[18] CAO S H, CHEN L, HUANG X H. Error analysis of a decoupled finite element method for quad-curl problems[J]. Journal of Scientific Computing, 2021, 90(1): 29.
-
期刊类型引用(1)
1. 孙兰银,庞琨琨. C-Bézier基函数在稳态线弹性方程求解中的应用. 信阳师范学院学报(自然科学版). 2024(02): 197-202 . 本站查看
其他类型引用(1)