带有复杂界面的二维热传导方程的数值模拟
科研笔记
数值模拟
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 |
区域外界的温度 |
考虑定义在 Ω 上的如下热传导模型:
⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪c∂T∂t=∇⋅(κ(T)∇T)+q−κ(T)∂T∂n=a(T−b)T(x,0)=T0 in Ω on ∂Ω in Ω
2. 变分形式
分别在 Ω+ 和 Ω− 上考虑上述模型的热传导变分形式:
⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪c+∂T+∂t=∇⋅(κ+∇T+)+q+c−∂T−∂t=∇⋅(κ−∇T−)+q−T(x,0)=T0,+T(x,0)=T0,−−κ(T+)∂T+∂n+=a(T+−b) in Ω+ in Ω− in Ω+ in Ω− on ∂Ω
T在界面 Γ 的函数值和通量跳跃:
[T]Γ=[κTn]Γ=T+−T−κ+∂T+∂n−κ−∂T−∂n
其中
n 为
Γ 上由
Ω− 指向
Ω+ 的单位法线向量.
⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪c+∂T+∂t=∇⋅(κ+∇T+)+q+c−∂T−∂t=∇⋅(κ−∇T−)+q−T(x,0)=T0,+T(x,0)=T0,−[T]Γ=T+−T−=p0−κ(T+)∂T+∂n+=a(T+−b) in Ω+ in Ω− in Ω+ in Ω− on Γ on ∂Ω
取 v∈H1(Ω), 做变分可得:
(c+∂T+∂t,v)Ω++(κ+∇T+,∇v)Ω+=(c−∂T−∂t,v)Ω−+(κ−∇T−,∇v)Ω−=(q+,v)Ω++(κ+∂T+∂n+,v)Γ+(κ+∂T+∂n+,v)∂Ω(q−,v)Ω−+(κ−∂T−∂n−,v)Γ
引入函数 ϕ−:Ω−→R, 且 ϕ−∈H1(Ω−) 和 ϕ−|Γ=p0. 记 ϕ 为 ϕ− 在Ω 的零扩展. 构造新的函数 u=T+ϕ, 则 u 满足:
(c+∂u∂t,v)Ω++(κ+∇u,∇v)Ω+=(c−∂u∂t,v)Ω−+(κ−∇u,∇v)Ω−=(q+,v)Ω++(κ+∂u∂n+,v)Γ+(κ+∂u∂n+,v)∂Ω(q−,v)Ω−+(κ−∂u∂n−,v)Γ+(c−∂ϕ∂t,v)Ω−+(κ−∇ϕ,∇v)Ω−
上面两式相加, 可得
=(c∂u∂t,v)Ω+(κ∇u,∇v)Ω(q,v)Ω+(κ+∂u∂n+,v)∂Ω++(κ−∂u∂n−,v)Γ+(c−∂ϕ∂t,v)Ω−+(κ−∇ϕ,∇v)Ω−
显格式
隐格式
=1τ(cun+1−cun,v)Ω+(κ∇un+1,∇v)Ω(q,v)Ω+(κ+∂un+1∂n+,v)∂Ω++(κ−∂un+1∂n−,v)Γ+1τ(c−ϕn+1−c−ϕn,v)Ω−+(κ−∇ϕn+1,∇v)Ω−
整理可得:
=1τ(cun+1,v)Ω+(κ∇un+1,∇v)Ω−(κ+∂un+1∂n+,v)∂Ω+−(κ−∂un+1∂n−,v)Γ−1τ(c−ϕn+1,v)Ω−−(κ−∇ϕn+1,∇v)Ω−1τ(cun,v)Ω+(qn+1,v)Ω+1τ(c−ϕn,v)Ω−
un+1=un=ϕn=ϕn+1=∑i=1Nun+1iψi∑i=1Nuniψi∑ϕniψi∑ϕn+1iψi
1τMUn+1+AUn+1−(C++C−)Un+1