引言

在过去的几十年中，单螺杆挤出机一直是聚合物加工最重要的设备之一，由于设计和操作的需要，对其流体流动的数值模拟研究一直是热点问题^[1]。

单螺杆挤出机内的熔体流动的数学模型一般为非牛顿流体的流动与热传导的耦合的非线性偏微分方程组定解问题^[2]，而且求解区域随时间变化，结构复杂。对此模型的数值求解通常用有限元^[3]等方法或简化为更低维模型求解，但要较精确高效地求解非常困难，主要困难来源于求解区域的复杂性。

单螺杆挤出机流体流动区域一般由机筒及螺杆组成，其中机筒内的螺杆还按一定的角速度旋转。为了克服求解区域的复杂带来的困难，文献[4]和[5]利用不同的变换(类似于拉直)将螺槽区域变换为柱体区域，从而简化了问题的计算，但难以应用于包含漏流层情形。

一般单螺杆挤出机都有漏流层，其流体流动区域的截面是一个二维双连通区域，通常其外边界曲线为圆周，内边界曲线(也称为螺杆的端面曲线)依螺槽的形状而不同，可以是比较复杂的曲线。对于均匀变化的螺杆或螺杆中均匀变化的一段，其不同截面可看成是由一个截面沿同一轴心旋转不同角度所形成。因此，若将每一横截面通过适当的变换变为相同的圆环形区域，则单螺杆挤出机流体流动区域就变为圆环柱形区域，从而使得求解区域不仅与时间无关，而且成为较规则的三维区域，非常便于计算包含漏流层的单螺杆挤出机流体流动问题。为了使得变换光滑，且变换后的方程尽可能简洁，共形变换^[6-8]是一种较好的变换方法。

对于不可压缩流体流动问题，近年来兴起的投影算法^[9]是一个相对高效的算法，其核心步骤是求解流动区域上的泊松方程。为了说明将共形变换应用于投影算法计算单螺杆挤出机流体流动问题时可较大简化问题的计算，本文以端面曲线为椭圆的单螺杆挤出机为例，计算了其流动区域上的三维泊松方程定解问题^[4]，并与精确解作了比较。计算结果表明，此方法简单有效。

1 共形映射的计算

通过共形映射^[8]使得区域规则化的方法有很多种，本文借鉴文献[6]的方法计算单螺杆挤出机流体流动区域截面到圆环区域的共形映射。

1.1 理论公式

在复变函数理论中，有以下熟知的结论^[6]。

定理1  如果h在复平面区域Ω上是单值解析函数，且在Ω∪Γ连续，Γ为Ω的边界且Γ是光滑的，则有

$ \frac{1}{{2{\rm{\pi i}}}}\int_\mathit{\Gamma } {\frac{{h\left( w \right)}}{{w - z}}{\rm{d}}w = \frac{1}{2}h\left( z \right),z \in \mathit{\Gamma }} $

(1)

设Ω为双连通区域，Ω的边界为Γ，方程为z=z(t)。其中外边界为Γ₀，方程为z=z₀(t)；内边界为Γ₁，方程为z=z₁(t)。理论上，对任一双连通区域，存在单值解析函数将其映射为同心圆环区域。设映射函数为f(z)，且f(z₀(t))=e^iθ₀(t), f(z₁(t))=μ₁e^iθ₁(t)，即f将Γ₀映为单位圆周，将Γ₁映为半径μ₁的圆周(1/μ₁为Ω的模，是个待定的参数)，从而f将Ω映射为同心圆环。若能解出θ₀(t), θ₁(t)和μ₁，则利用Cauchy积分公式即可得到f(z)。

利用定理1, 加上适当的推导，对于双连通区域，可得到

$ \begin{array}{l} \;\;\;\;\;\;\;g\left( {{z_0}} \right) + \int_{{\mathit{\Gamma }_0}} {N\left( {{z_0},w} \right)g\left( w \right)\left| {{\rm{d}}w} \right|} - \int_{ - {\mathit{\Gamma }_1}} {{P_0}\left( {{z_0},} \right.} \\ \left. w \right)g\left( w \right)\left| {{\rm{d}}w} \right| = 0,{z_0} \in {\mathit{\Gamma }_0} \end{array} $

(2)

$ \begin{array}{l} \;\;\;\;\;\;\;g\left( {{z_1}} \right) + \int_{{\mathit{\Gamma }_0}} {{P_1}\left( {{z_1},w} \right)g\left( w \right)\left| {{\rm{d}}w} \right|} - \int_{ - {\mathit{\Gamma }_1}} {N\left( {{z_1},} \right.} \\ \left. w \right)g\left( w \right)\left| {{\rm{d}}w} \right| = 0,{z_1} \in {\mathit{\Gamma }_1} \end{array} $

(3)

其中T(z)=z(t)/|z′(t)|, g(z)=T(z)f′(z)，P₀(z, w)、P₁(z, w)、N(z, w)是与z、w、T(z)、μ₁相关的函数。

为了保证唯一性，可设f(z₀(0))=1, 即

$ \left\{ \begin{array}{l} {\mathop{\rm Re}\nolimits} \left[ {g\left( {{z_0}\left( 0 \right)} \right)} \right] = 0\\ {\mathop{\rm Im}\nolimits} \left[ {g\left( {{z_0}\left( 0 \right)} \right)/\left| {g\left( {{z_0}\left( 0 \right)} \right)} \right|} \right] = 1 \end{array} \right. $

(4)

最后基于边界对应函数的单调递增性，有

$ \left\{ \begin{array}{l} \int_0^{2{\rm{\pi }}} {\left| {g\left( {{z_0}\left( 0 \right)} \right){{z'}_0}\left( t \right)} \right|{\rm{d}}t} = \int_0^{2{\rm{\pi }}} {{{\theta '}_0}\left( t \right){\rm{d}}t} = 2{\rm{\pi }}\\ \int_0^{2{\rm{\pi }}} {\left| {g\left( {{z_1}\left( t \right)} \right){{z'}_1}\left( t \right)} \right|{\rm{d}}t} = 2{\rm{\pi }}{\mu _1} \end{array} \right. $

(5)

1.2 数值求解

设参数t的变化范围为[0, 2π]。在Γ₀上对t等距离散n个点，在Γ₁上对t等距离散m个点，用梯形法则的Nystrom方法将式(2)~(5) 离散化，得到

$ \left( {\begin{array}{*{20}{c}} {{\mathop{\rm Re}\nolimits} A}&{ - {\mathop{\rm Im}\nolimits} A}\\ {{\mathop{\rm Im}\nolimits} A}&{{\mathop{\rm Re}\nolimits} A} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{\mathop{\rm Re}\nolimits} x}\\ {{\mathop{\rm Im}\nolimits} x} \end{array}} \right) = \left( \begin{array}{l} 0\\ 0 \end{array} \right) $

(6)

$ \sum\limits_{j = 1}^m {\left( {{\mathop{\rm Re}\nolimits} {x_{1\tilde j}} + {\rm{i}}{\mathop{\rm Im}\nolimits} {x_{1\tilde j}}} \right)} = 0 $

(7)

$ \sum\limits_{j = 1}^n {\sqrt {{{\left( {{\mathop{\rm Re}\nolimits} {x_{0j}}} \right)}^2} + {{\left( {{\mathop{\rm Im}\nolimits} {x_{0j}}} \right)}^2}} } = n $

(8)

$ \sum\limits_{\tilde j = 1}^m {\sqrt {{{\left( {{\mathop{\rm Re}\nolimits} {x_{1\tilde j}}} \right)}^2} + {{\left( {{\mathop{\rm Im}\nolimits} {x_{1\tilde j}}} \right)}^2}} } = m{\mu _1} $

(9)

$ {\mathop{\rm Re}\nolimits} {x_{01}} = 0 $

(10)

$ {\mathop{\rm Im}\nolimits} \left[ {{x_{01}}/\sqrt {{{\left( {{\mathop{\rm Re}\nolimits} {x_{01}}} \right)}^2} + {{\left( {{\mathop{\rm Im}\nolimits} {x_{01}}} \right)}^2}} } \right] = 1 $

(11)

其中x_0j=|z′₀(t_j)|g(z₀(t_j))(j=1, 2, …, n), $ {x_{1{\rm{ }}\tilde j}}=|{z'_0}({t_{\tilde j}})|g({z_0}({t_{\tilde j}}))\left( {\tilde j=1, 2, \ldots, m} \right) $, $ \mathit{\boldsymbol{A}}=\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_n} + \mathit{\boldsymbol{B}}}&{ - \mathit{\boldsymbol{C}}}\\ \mathit{\boldsymbol{D}}&{{\mathit{\boldsymbol{I}}_m} - \mathit{\boldsymbol{E}}} \end{array}} \right) $, B=h₀(|z′₀(t_i)|N(z₀(t_i), z₀(t_j)))_n×n, C=h₁(|z′₀(t_i)|P₀(z₀(t_i), z₀(t_j)))_n×n, D=h₀(|z′₀($ {t_{\tilde i}} $)|N(z₀($ {t_{\tilde i}} $), z₀(t_j)))_n×n, C=h₁(|z′₀($ {t_{\tilde i}} $)|P₀(z₀($ {t_{\tilde i}} $), z₀($ {t_{\tilde j}} $)))_n×n, h₀=2π/n, h₁=2π/m, x=(x₀₁, x₀₂, …, x_0n, x₁₁, x₁₂, …, x_1m)^T。

x分量的实部和虚部加上参数μ₁，共有2(m+n)+1个未知量，式(6)~(11) 共有2(m+n)+5个方程，所以此方程组为超定方程组，可采用最小二乘法求解。

解出x后，即可求得θ₀(t)、θ₁(t)，从而得到边界上的点在映射f下的值。对于内部的点，利用柯西积分可计算出f的值，即

$ \begin{array}{l} \;\;\;\;\;\;\;f\left( {zz} \right) = \frac{1}{{2{\rm{\pi i}}}}\int_0^{2{\rm{\pi }}} {\frac{{f\left( {{z_0}\left( t \right)} \right)}}{{{z_0}\left( t \right) - zz}}{{z'}_0}\left( t \right){\rm{d}}t} - \frac{1}{{2{\rm{\pi i}}}}\int_0^{2{\rm{\pi }}} {} \\ \frac{{f\left( {{z_1}\left( t \right)} \right)}}{{{z_1}\left( t \right) - zz}}{{z'}_1}\left( t \right){\rm{d}}t,zz \in \mathit{\Omega } \end{array} $

以Ω的外边界为单位圆周，内边界为椭圆周为例进行了计算，用大范围信赖域方法求解式(6)~(11) 构成的超定方程组的非线性最小二乘问题。当内边界椭圆长半轴为0.9、短半轴为0.6时，取m=n=48，计算得到μ₁=0.783176；计算的最优结果下，最小二乘的残差平方和为4.136×10^-9。图 1为计算结果图形。可以看出，此方法能得到较高精度的共形映射。

2 三维流动区域的变换

对于均匀变化的螺杆或螺杆中均匀变化的一段，其不同截面可看成是由一个截面沿同一轴心旋转不同角度形成，如图 2所示。将每个截面通过第1章的共形变换变为相同的圆环形区域，则单螺杆挤出机流体流动区域就变为圆环柱形区域。

设机筒为圆柱面，螺杆的轴心与机筒的轴心一致，以此轴心为z轴，z=0为流体进口的横截面，并取其为xOy坐标面，建立直角坐标系{O, x, y, z}。以挤出机的螺槽是单线螺槽的情况为例，螺距为H₀。当z=0时，共形映射f将流体流动截面区域Ω₀映射为圆环形区域D。取b=H₀/(2π), 则螺杆表面可看成是Ω₀的内边界Γ₁随z的增大而绕z轴沿逆时针方向旋转角度z/b得到的，且当z=H₀时恰好转过2π弧度角；当z=h时, 流体流动截面区域Ω_h也可看成是Ω₀绕z轴旋转角度h/b得到。这样，流体流动区域G上的变换可以取为

$ \left\{ \begin{array}{l} \xi = {\mathop{\rm Re}\nolimits} \left( {{{\rm{e}}^{{\rm{i}}g\left( {x,y,z} \right)}}f\left( {\left( {x + {\rm{i}}y} \right){{\rm{e}}^{ - {\rm{i}}z/b}}} \right)} \right)\\ \eta = {\mathop{\rm Im}\nolimits} \left( {{{\rm{e}}^{{\rm{i}}g\left( {x,y,z} \right)}}f\left( {\left( {x + {\rm{i}}y} \right){{\rm{e}}^{ - {\rm{i}}z/b}}} \right)} \right)\\ h = z \end{array} \right. $

(12)

其中g(x, y, z)可取为任意的实值光滑函数。这时，区域G就变为{ξ, η, h}下的圆环柱体区域。若G={(x, y, z)|(x, y)∈Ω_z, 0≤z≤H}，对{ξ, η, h}再作柱坐标变换

$ \xi = r\cos s,\eta = r\sin s,h = h $

(13)

则G变为{r, s, h}坐标下的长方体区域。

3 求解泊松方程

为了说明第二章的三维变换对求解复杂的单螺杆挤出机三维流体流动区域G上的偏微分方程有效和简便^[10]，以端面曲线为椭圆的单螺杆挤出机为例，应用该变换求解G上的Dirichlet问题，并与精确解作比较。

3.1 问题变换与计算方法

设G={(x, y, z)|(x, y)∈Ω_z, 0≤z≤H}，在G区域上所求解的问题为：

$ \left\{ \begin{array}{l} {U_{xx}} + {U_{yy}} + {U_{zz}} = p\left( {x,y,z} \right)\\ U\left| {_{\partial G} = q\left( {x,y,z} \right)} \right. \end{array} \right. $

(14)

设H₀为螺距，b=H₀/2π为参数，Ω₀为进口(z=0) 处的流体流动截面区域，则Ω_z为Ω₀绕z轴逆时针旋转角度z/b得到的流体流动截面区域。设Ω₀的外边界为单位圆周，参数方程为x=cos t和y=sin t(t∈[0, 2π]); 内边界为椭圆周，其方程为x=ccos t和y=dsin t(t∈[0, 2π])，其中c, d > 0为椭圆的两个半轴长。

设共形映射f将Ω₀映成内半径为μ₁、外半径为1的圆环区域，且取

$ \left\{ \begin{array}{l} \xi = {\mathop{\rm Re}\nolimits} \left( {{{\rm{e}}^{{\rm{i}}z/b}}f\left( {\left( {x + {\rm{i}}y} \right){{\rm{e}}^{ - {\rm{i}}z/b}}} \right)} \right)\\ \eta = {\mathop{\rm Im}\nolimits} \left( {{{\rm{e}}^{{\rm{i}}z/b}}f\left( {\left( {x + {\rm{i}}y} \right){{\rm{e}}^{ - {\rm{i}}z/b}}} \right)} \right)\\ h = z \end{array} \right. $

(15)

则此变换将G映射成内半径为μ₁、外半径为1、高为H的圆环柱V。进一步作式(13) 所示的柱坐标变换，设变量复合后变换为

$ \left\{ \begin{array}{l} r = r\left( {x,y,z} \right)\\ s = s\left( {x,y,z} \right)\\ h = z \end{array} \right. $

(16)

则G变换为{r, s, h}坐标下的长方体区域B=[μ₁, 1]×[0, 2π]×[0, H], 问题(14) 变成

$ \left\{ \begin{array}{l} {a_{11}}{U_{rr}} + {a_{21}}{U_{ss}} + {a_{33}}{U_{hh}} + 2{a_{12}}{U_{rs}} + 2{a_{13}}{U_{rh}} + \\ \;\;\;\;\;\;\;2{a_{23}}{U_{sh}} + {b_1}{U_r} + {b_2}{U_s} + {b_3}{U_h} + cU = \\ \;\;\;\;\;\;\;\tilde p\left( {r,s,h} \right)\left( {\left( {r,s,h} \right) \in B} \right)\\ U\left| {_{r = {\mu _1}} = {g_{1a}}\left( {s,h} \right)} \right.,U\left| {_{r = 1} = {g_{1b}}\left( {s,h} \right)} \right.\\ U\left| {_{s = 0}} \right. = U\left| {_{s = 2{\rm{\pi }}}} \right.,{U_s}\left| {_{s = 0}} \right. = {U_s}\left| {_{s = 2{\rm{\pi }}}} \right.\\ U\left| {_{h = 0} = {g_{3a}}\left( {r,s} \right)} \right.,U\left| {_{h = H} = {g_{3b}}\left( {r,s} \right)} \right. \end{array} \right. $

(17)

其中a₁₁=r_x²+r_y²+r_z², a₂₂=s_x²+s_y²+s_z², a₃₃=h_x²+h_y²+h_z²=1, a₁₂=r_xs_x+r_ys_y+r_zs_z, a₁₃=r_xh_x+r_yh_y+r_zh_z=r_z, a₂₃=s_xh_x+s_yh_y+s_zh_z=s_z, b₁=r_xx+r_yy+r_zz, b₂=s_xx+s_yy+s_zz, b₃=h_xx+h_yy+h_zz=0, c=0

式(17) 的方程形式虽然比式(14) 复杂，但不会产生实质性困难，反而因为复杂的求解区域变成了简单的长方体区域，可方便地用差分方法求解，极大地简化了问题的计算。

计算步骤如下：

1) 计算共形映射f  得到θ₀(t), θ₁(t)和μ₁后，以样条函数形式给出w₀(t)=f(z₀(t))=e^iθ₀(t)和w₁(t)=f(z₁(t))=μ₁e^iθ₁(t)；

2) 计算x=x(r, s, h)和y=y(r, s, h)  因为f的逆映射f^-1也是解析函数，可利用Cauchy公式给出f^-1，进一步，由

$ x + {\rm{i}}y = {{\rm{e}}^{{\rm{i}}h/b}}{f^{ - 1}}\left( {\left( {r\cos s + {\rm{i}}r\sin s} \right){{\rm{e}}^{ - {\rm{i}}h/b}}} \right) $

计算出x=x(r, s, h)和y=y(r, s, h)，并以三元样条函数形式给出；

3) 对式(17) 离散计算  采用等步长中心差分格式，得到的线性方程组用块追赶法求解。

3.2 数值计算结果

设精确解为U_e=1+x²y+sin z，先计算出相应的p(x, y, z)=2y-sinz和q(x, y, z)，然后分别就c=0.5, d=0.8和c=0.6, d=0.9对式(17) 进行计算。令U_ij=max_k{|U_ijk-U_eijk|}，其中U_ijk为U(r_i, s_j, h_k)的计算值，U_eijk为U_e(r_i, s_j, h_k)的值。

计算误差的情况如下。

1) c=0.5, d=0.8

当(r, s, h)网格数取为5×10×100时，最大模误差为0.02793，均方根误差为0.00534；网格数取为10×20×200时，最大模误差为0.02748，均方根误差为0.00145。h方向的最大模误差分布如图 3所示。

2) c=0.6, d=0.9

当(r, s, h)网格数取为5×10×200时，最大模误差为0.10456，均方根误差为0.00862；网格数取为10×20×400时，最大模误差为0.08182，均方根误差为0.00497。h方向的最大模误差分布如图 4所示。

4 结论

通过本文所述的共形变换将流动区域的截面变为圆环形区域，再构造三维区域变换，可以将单螺杆挤出机流体流动区域变为圆环柱形区域，从而使得求解区域不仅与时间无关，而且变成了较规则的三维区域，大大便于问题的计算；对此这类区域上的三维泊松方程进行数值实验的结果表明，此方法简单有效。当使用投影算法求解不可压缩流体流动问题时，利用本文方法求解核心步骤中的泊松方程，可以使投影算法更为简便高效。