带有复杂界面的二维热传导方程的数值模拟-无辅助函数
科研笔记 数值模拟
1. 问题模型
给定一个由模具区域和铸件区域组成的二维区域 Ω⊂R2, 记模具区域为 Ω+, 铸件区域为 Ω−, 两者的分界面为 Γ. 这里假定 Γ 可由一个水平集函数 ϕ(x) 表示, 则有 Ω=Ω−∪Γ∪Ω+, 其中:
Γ=Ω+=Ω−={x∈Ω|ϕ(x)=0}{x∈Ω|ϕ(x)>0}{x∈Ω|ϕ(x)<0}
记
| 符号 |
意义 |
| T(x,t) |
区域 Ω 上的温度分布变化函数 |
| T0=T(x,0) |
区域 Ω 上的初始温度分布函数 |
| T0,+ |
模具的初始温度 |
| T0,− |
铸件的初始温度 |
| T+ |
T 在模具区域 Ω+ 上的限制 |
| T− |
T 在铸件区域 Ω− 上的限制 |
| q(x,t) |
区域 Ω 上的热源分布变化函数 |
| q+ |
q 在 Ω+ 上的限制 |
| q−− |
q 在 Ω− 上的限制 |
| c(x,t)=sρ |
比热和密度的乘积函数, 实际为单位质量材料升高一度需要的热量 |
| c+ |
c 在 Ω+ 的限制 |
| c− |
c 在 Ω− 的限制 |
| κ(T) |
材料的导热系数, 它是关于材料温度 T 的函数 |
| κ+=κ(T+) |
κ 在 Ω+ 上的限制 |
| κ−=κ(T−) |
κ 在 Ω− 上的限制 |
| a |
外边界∂Ω 的对流系数 |
| b |
区域外界的温度 |
| H1(Ω) |
|
| H1(Ω−∪Ω) |
|
考虑定义在 Ω 上的如下热传导模型:
⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪c∂T∂t=∇⋅(κ(T)∇T)+q−κ(T)∂T∂n=a(T−b)T(x,0)=T0 in Ω on ∂Ω in Ω
2. 变分形式
T在界面 Γ 的函数值间断, 这里记:
[T]Γ=T+−T−
但 T 在 Γ 上的通量连续， 即：
κ−∂T−∂n=κ+∂T+∂n=β(T+−T−)
这里
n 表
Γ 上由
Ω− 指向
Ω+ 的单位法线向量。
⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪c+∂T+∂t−∇⋅(κ+∇T+)=q+T+(x,0)=T0,+c−∂T−∂t−∇⋅(κ−∇T−)=q−T−(x,0)=T0,−κ−∂T−∂n=κ+∂T+∂n=β(T+−T−)−κ+∂T+∂n+=a(T+−b) in Ω+ in Ω+ in Ω− in Ω− in Γ on ∂Ω
取 v∈H1(Ω), 分别在 Ω+ 和 Ω− 上考虑上述模型的热传导变分形式, 可得到如下的偶合连续变分问题:
(c+∂T+∂t,v)Ω++(κ+∇T+,∇v)Ω+=(c−∂T−∂t,v)Ω−+(κ−∇T−,∇v)Ω−=(κ−∂T−∂n,v)Γ=(q+,v)Ω+−(κ+∂T+∂n,v)Γ−(aT+,v)∂Ω+(ab,v)∂Ω(q−,v)Ω−+(κ−∂T−∂n,v)Γ(κ+∂T+∂n,v)Γ
(c+∂T+∂t,v)Ω++(κ+∇T+,∇v)Ω+=(c−∂T−∂t,v)Ω−+(κ−∇T−,∇v)Ω−=(q+,v)Ω+−(βT+,v)Γ+(βT−,v)Γ−(aT+,v)∂Ω+(ab,v)∂Ω(q−,v)Ω−+(βT+,v)Γ−(βT−,v)Γ
3. 界面拟合网格上的线性有限元离散
下面考虑上述连续变分问题在三角形界面拟合网格T上的线性有限元的离散形式, 引入线性有限元空间 Vh(Ωh), 这里 Ωh 是 T 覆盖的区域对Ω 的逼近.
因为是界面拟合网格, 所以 T 存在一条由网格边组成的折线 Γh (是Γ的逼近), 把T分成两部分 T− 和 T+, 它们覆盖的区域分别记为 Ωh,− (Ω− 的逼近) 和 Ωh,+(Ω+的逼近).
设 T 有 N 个节点, T− 中有 N− 个节点, T+ 中有 N+ 个节点, Γh 上的有 NΓ 个节点, 则
N=N−+N+−NΓ
引入定义在 T 的 N 维线性有限元空间
Vh(Ωh)=span{ϕi,i=1,⋯,N}
把 Vh 分别限制在 Ωh,− 和 Ωh,+, 可得相应的 N− 维空间 Vh,− 和 N+ 维空间 Vh,+, 引入如下的有限元函数:
Tn+,h=∑i=1N+T+,inϕi,Tn+1+,h=∑i=1N+Tn+1+,iϕi∈VhTn−,h=∑i=1N−Tn−,iϕi,Tn+1−,h=∑i=1N−Tn+1−,iϕi∈Vh
其中上标 n 和 n+1 分别表示第 n 个和第 n+1 个时间层, 记 δn 为第 n 个时间层到第 n+1 个时间层的步长.
记 M− 和 M+ 为为铸件区域和模具区域上质量矩阵, A− 和 A+ 为铸件区域和模具区域上刚度矩阵, C− 和 C+ 分别为铸件区域和模具区域在边界 Γ 上的通量矩阵, C 为外边界上的由通量条件引入的矩阵, U+ 和 U− 分别为铸件区域和模具区域上的解向量, F− 和 F+ 为为铸件区域和模具区域上热源向量(在本模型中都取零向量), F0 由外边界通量条件引入的常向量.
下面给出不同的离散格式的数学形式.
3.1 显格式
M+Un+1+=M−Un+1−=M+Un+−δn(A++C++C)Un++δnC−Un−+δnF0M−Un−−δn(A−+C−)Un−+C+Un+
3.2 隐格式
M+Un+1++δn(A++C++C)Un+1+−δnC−Un+1−=M−Un+1−+δn(A−+C−)Un+1−−δnC+Un+1+=M+Un++δnF0M−Un−
3.3 Crank-Nicolson
