CN103218471B - 用于模拟塑性材料和气体的多相流的方法 - Google Patents

Abstract

一种用于模拟腔室中的塑性材料的流动状态的计算机化方法，其包括下述步骤：生成作为腔室的有限体积模型的腔室模型；定义塑性材料的材料模型；定义气体的气相模型；以及对腔室模型中的材料模型和气相模型进行流动计算，其中，所述材料模型设定有塑性材料的粘度、并且所述气相模型设定有比气体实际粘度更大的粘度。

Description

用于模拟塑性材料和气体的多相流的方法

技术领域

本发明涉及用于模拟塑性材料和气体的多相流的计算机化方法，其有助于分析于腔室（Chamber）中混炼的例如未硫化橡胶、树脂等塑性材料的流动状态。

背景技术

近年来，已经提出了各种计算机模拟法。例如，下述非专利文献1提出了用于计算在班伯里密炼机中混炼时在交联之前的塑性材料例如未硫化橡胶或树脂的流动状态的方法。

[非专利文献1]“附聚物的分散混合的数值和实验研究”，V.Collin(1)^*、E.Peuvrel-Disdier(1)等人。

非专利文献1的方法的前提是腔室完全被塑性材料充满（即填充率=100%）。

然而，通过使用班伯里密炼机等难以使混炼在填充率为100%的条件下进行，并且已知随着填充率的改变，混炼的塑性材料的品质和状态也改变。

尽管存在这些事实，但非专利文献1的方法不能通过将塑性材料的填充率设定为除100%以外的不同值来实施。的确，要使混炼在填充率100%的条件下进行几乎是不可能的，因此，非专利文献1的模拟结果不能通过实验来验证。

另一方面，塑性材料的粘度很大程度上不同于空气的粘度，因此，如果通过流动计算（Flowcalculation）来模拟腔室中的塑性材料和空气两相，则计算倾向于变得不稳定。在通过利用空气的实际粘度进行流动计算的情况下，即使在不被打断的情况下完成计算，计算时间也变得非常长，并且需要高昂的计算成本。

发明内容

因此，研究出本发明来解决上述问题，并且本发明的目的在于提供一种用于模拟存在于腔室中的塑性材料和气体的多相流的计算机化方法，在所述方法中，可以在不大幅增加计算时间的条件下稳定地进行各种填充率下的多相流计算，并且由此可实现精确模拟。

根据本发明，用于模拟在腔室中的塑性材料的流动状态的计算机化方法，包括：

生成作为腔室的有限体积模型的腔室模型的步骤、

定义塑性材料的材料模型的步骤、

定义气体的气相模型的步骤、

对在腔室模型中的材料模型和气相模型进行流动计算的步骤，其中，

材料模型设定有塑性材料的粘度，并且

气相模型设定有大于气体实际粘度的粘度。

所述腔室可以是班伯里密炼机的混炼空间，在其中设置至少一个转子，并且塑性材料可以是待混炼的未硫化橡胶或树脂材料。

在流动计算中，优选进行材料模型的剪切生热的计算，但是不进行气相模型的剪切生热的计算。

在塑性材料模型和腔室模型之间的界面上，可对滑移边界条件进行定义。

附图说明

图1是用于混炼塑性材料的班伯里密炼机一部分的示意性横截面图。

图2是根据本发明的方法的一个例子的流程图。

图3是腔室模型的透视图。

图4是用于说明腔室模型的腔壁（Wall）和内部空间的啮合的横截面图。

图5是单独显示腔室模型的功能部分的横截面图。

图6是显示配置了材料模型和气相模型的腔室模型状态的横截面图。

图7是流动计算程序的流程图。

图8是显示通过根据本发明的模拟方法获得的塑性材料的混炼状态的图形。

具体实施方式

现在结合附图对本发明的实施方式加以详细说明。

本发明的目的在于提供一种用于通过使用计算机（未图示）来评估设置于腔室中的塑性材料的流动状态的模拟方法。

塑性材料，例如交联之前的未硫化橡胶、树脂或弹性体，具有粘度。然而，其并不限于这些材料。只要其具有稳定的流动状态，所述塑性材料可以是任意材料。

在为交联之前未硫化橡胶的情况下，可将在80℃左右充分混炼的混合物的状态作为稳定的流动状态。

所述腔室为基本上闭合的空间，在其中设置有塑性材料，并且所述塑性材在其中流动。所述腔室可以具有任意的构造（Configuration）。所述腔室的一个典型例子是班伯里密炼机的混炼空间。如图1所示，本实施方式中的腔室4形成于班伯里密炼机1的外壳2和设置于外壳2中的一对转子3之间，并且所述腔室4具有像数字8的截面形状。然而，所述腔室并不限于这样的构造。

图2是根据本发明的模拟法的一个例子的流程图。

*生成腔室模型的步骤s1。

在该实施方式中，首先，计算机生成腔室模型5，其是腔壁和腔室4内部三维空间的有限体积模型。

图3是腔室模型5的透视图。

图4是其横截面图。

腔室模型5具有由下述部分封闭而成的三维空间：由外壳2的内表面限定的外周面5o、由两个旋转的转子3的外周面限定的内周面5i、以及转子3沿轴向两侧的两个端面5s。

该三维闭合空间被分（被离散）成三维单元（e）例如四面体单元、六面体单元和多面体单元。

对于各单元（Element），计算塑性材料（材料模型）的物理量例如压力、温度和/或速度。

所述外周面5o和两个端面5s不被移动（Move）。然而，内周面5i根据转子3的旋转而移动，并且相应地，腔室模型5的构造被改变。

在图5中显示的例子中，腔室模型5由四个功能部分组成：一对旋转部分5A及5B；夹在它们之间的中间部分（In-between part）5C；以及包绕这些部分（5A、5B和5C）的外框架部分5D。

各旋转部分5A/5B皆为筒状并且具有圆筒形外周面5Ao/5Bo以及对应于一个转子3外周面的内周面5i。

将旋转部分5A和5B设置在外框架部分5D中，并且将它们限定为可环绕各中轴线Oa和Ob旋转，来显示由转子3的旋转引起的腔室体积构造的变化。

与此相反，所述中间部分5C保持固定在旋转部分5A和5B之间，并且具有紧靠于各旋转部分5A和5B的两个凹面j。在所述凹面j和各圆筒形外周面5Ao和5Bo上，定义了作为滑移面的边界条件。

这使得发生在旋转部分5A和5B之中的物理作用（力、热等）经过凹面j被传递至存在于中间部分5C之中的材料模型。

外框架部分5D为筒状并且包绕旋转部分5A、5B和中间部分5C。两个轴末端均被两个端面5s封闭。

在外框架部分5D和旋转部分5A与5B之间的界面上，以及在外框架部分5D和中间部分5C之间的界面上，定义了作为滑移面的边界条件。

这使得产生于旋转部分5A和5B中的物理作用（力、热等）穿过它们之间的界面被传递至外框架部件5D。

外框架部分5D因转子的运行而受到相对较大的剪切力。因此，为了更详细地计算材料的速度等，将组成外框架部分5D的元件的尺寸制得比旋转部分5A和5B以及中间部分5C的尺寸更小。由此，可更详细地计算靠近腔室模型5内表面的材料模型的速度分布图等。

定义材料模型的步骤s2^*

接着，通过计算机定义并且储存材料模型。材料模型是在腔室4中流动或移动的塑性材料的模型。

在材料模型上，对塑性材料的物理特性例如剪切粘度、比热、热导率和比重进行定义并且储存于计算机中。

如上所述，本实施方式中的塑性材料是处于塑化状态的橡胶混合物，因此，预先输入并且储存这样的橡胶混合物的物理特性。

[剪切粘度]

为了得到剪切粘度，首先，在多个温度条件下测定塑性材料的粘弹性（G′和G″）。然后，根据Cox-Merz规则，将测得的粘弹性换算成剪切粘度。

根据幂次法则，如上得到的剪切粘度可通过下述式（1）近似得到。

式（1）：

η=mγ′^n-1

其中，

η：剪切粘度、

m：绝对温度的函数的系数、

γ′：剪切速度、以及

n：系数。

分析对象（塑性材料）的比热可以通过绝热连续加热法（以25℃）得到，并且预先将测得的比热值输入并储存在电脑中。

分析对象（塑性材料）的热导率可以通过热线法（以25℃）得到，并且预先将测得的值输入并储存在电脑中。

定义气相模型的步骤s3^*

接着，通过计算机定义并且储存气相模型。所述气相模型是存在于腔室中的气体的模型。

在本发明中，在腔室中塑性模型的填充率小于100%，因此，为了实现流动计算，未充满塑性模型的部分用气相模型充满。

本实施例中的气相模型设定了粘度和比重。

将气相模型的比重设定为气体比重的实际值。

然而，通过计算机将气相模式的粘度设定为大于气体粘度实际值的值。

当腔室中的塑性材料的填充率小于100%时，有必要对塑性材料和气体的多相流进行流动计算。

通常，气体和塑性材料之间的粘度差别非常大。因此，如果将因为速度引起的剪切生热纳入考虑来进行计算，则在气相模型和材料模型之间的界面上剪切生热会超常增加，并且结果是计算容易变得不稳定。

根据本发明，因此，为了使计算稳定，在不会给计算结果带来负面影响的范围内，尽可能地提高气相模型的粘度。结果是，有可能模拟中增加时间间隔（时间增量）来加快计算速度。

本发明的发明人仅改变气相模型的粘度，进行模拟来比较腔室模型5中的压力场。

结果，发现如果气相模型的粘度超过实际粘度的10倍，则腔室模型5中的压力开始增加，并且给压力场带来负面影响。

然而，如果气相模型的粘度低于实际粘度的5倍，则变得难以使计算稳定。

因此，优选将气相模型的粘度设定在气体实际粘度的5～10倍的范围内。

定义边界条件的步骤s4^*

接着，对流动计算所需的包括边界条件在内的各种条件进行定义。

在腔室模型5的表面上，边界条件可能包括流速边界条件和温度边界条件。

关于流速边界条件，可根据模拟的目的、所需的精度等对（a）壁表面无滑移条件，或（b）壁表面滑移条件进行定义。

在壁表面非滑移条件的情况下，在腔室模型5的表面上的材料模型的流速始终是零。

在壁表面滑移条件的情况下，在腔室模型5的表面上的塑性模型的流速可具有非零值、也可具有零值。

在这种情况下，可将下述式（2）（Navier′sLaw）应用于在材料模型和腔室模型5之间的界面上的滑移现象。

式（2）：

f_s＝-F_slip(v_wall-v_s)|v_s-v_wall|^eslip-1

式（2），：

τ_ω＝F_slip[v_slip-v_wall]

其中

f_s：在腔室模型5的表面上的材料模型的剪应力、

v_slip：在腔室模型5的表面上的材料模型的滑移速度(切向速度)、

v_wall：在腔室模型5的表面上的材料模型的切向速度、

F_slip、e_slip：材料特性(如果e_slip=1，则为线性；如果0<e_slip<1，则为幂次法则)。

为了定义在腔室模型5的表面上的材料模型的滑移速度v_slip，首先，得到初始阶段的速度场，然后由该速度场（或修正后的速度场）得到在腔室模型5的表面上的材料模型的剪切力，然后通过上述式（2）得到在腔室模型5的表面上的材料模型的滑移速度。

在此情况下，如果需要，为了控制滑移速度的突然增加，优选通过将滑移速度乘以小于1的松弛系数，在常规的SIMPLE算法的循环之中使腔室模型5各表面的壁表面滑移速度逐渐收敛于真实值。（见后述）

关于温度边界条件，可定义为：（a）其中热量不会经过腔室模型5的表面从腔室模型5逃逸至外界的绝热条件，或（b）其中整个腔室模型5的表面皆具有恒温（例如50℃）的条件。

进一步地，所述条件可能包括塑性模型的初始温度（在本实施方式中为20℃（293K））、由旋转部分5A和5B的转数代表的转子的转数、腔室模型5的表面的滑移比、塑性模型相对于腔室模型5的体积的填充率等。

此外，所述条件可以包括流动计算的初始状态、用于计算的时间间隔、内部处理中的迭代数、计算的最大周期等。

关于初始状态，例如参见图6中所示，可定义为

气相模型的区域A，所述区域A在水平界面s的上侧，所述水平界面被定义为延伸穿过腔室模型5，以及

液体模型的区域M，所述区域M在水平界面s的下侧。

因此，通过改变界面s的水平面，可调节塑性材料（材料模型）的填充率。可根据模拟的目的任意地定义这些条件。

进行流动计算的步骤s5^*

接着，在腔室模型5中设置材料模型和气相模型，并且基于上述条件进行流动计算。

在流动计算中，计算限定材料模型运动状态的沿三个坐标轴方向（x、y、z）的分速度、以及限定材料模型内部状态的温度T和压力p。因此，待测的未知量为这5个物理量。

本实施方式中的流动计算是通过使用用于不可压缩性流动的纳维尔斯托克斯方程来进行的，将各个材料模型和气相模型的密度设定为恒定值。

在本发明中，在整个可能的温度范围内，将材料模型视作液体。

因此，纳维尔斯托克斯方程，即质量守恒和能量方程的公式，必须作为流体的方程组来求解。

在本实施方式中，为了处理气体和塑性材料的多相流，采用用于计算自由界面的流动的VOF（流体体积）法。VOF法不直接计算两种流体之间的界面的运动。在VOF法中，自由界面通过对在各单元的体积范围内的各流体的填充率（体积分数）进行定义而表示。

支配方程如下。

[运动方程]

在本实施方式中，将在腔室模型5中气相模型和材料模型汇流的两相流视为单相流。

在此情况下，待求解的运动方程为沿三个坐标轴方向x、y和z的下述式（3）。通过VOF法对所述双相进行平均，并将其视为单相，结果是上述运动方程的求解变得可能。

式(3)：

其中，

u：多相流模型的速度、

p：多相流模型的压力、

ρ：多相流模型的密度、

g：重力加速度、

T：多相流模型的绝对温度、

F：外力。

在此处使用的在材料模型和气相模型之间界面附近的密度ρ和粘度系数μ的值通过各相（即材料模型和气相模型）所占的体积来确定其权重，并且然后通过显示在下述式（4）中的公式进行平均。

式(4)：

ρ＝∑α_qρ_qμ＝∑α_qμ_q

其中

α_q：各单元中各相的体积分数、

ρ_q：各单元中各相的密度、

μ_q：各单元中各相的粘度。

[质量守恒方程]

关于质量守恒方程（连续性方程）和压力方程，仅求解沿三个坐标轴方向的一组方程就足够。因此，根据本实施方式中的模拟法，即使是多相，计算机也能将其视作单相来计算流场（flow field）。换句话说，待求解的是其中材料特性随着位置（体积分数）改变的流动。

各相的位置可以根据由计算结果得到的体积分数的分布来估计。

[能量方程]

可以通过下述能量方程(5)得到材料模型的温度。

方程(5)：

其中

E：焓、

k：热导率、

S：源项（sourceterm）。

[体积分数的输送方程]

体积分数的分布决定了两相（或气相模型和材料模型）之间的界面的位置。

可以通过精确求解下述式（6）来得到体积分数αq。

式(6)：

如果在任意单元（e）中的体积分数αq=0，则意味着q相（q相：材料模型）在该单元（e）中不存在。

如果体积分数aq=1，则意味着单元（e）的整个体积都被q相充满。

如果0<αq<1，则意味着单元(e)被q相和其他相（气相模型）充满，即，所述单元为多相，其具有界面。

此方程可以通过修正的HRIC(隐式)求解，在“ANSYS Fluent User′s Manual，26.2.9 Modified HRIC Scheme”给出了其详细的描述。

在本实施方式中，上述的各方程通过基于压力的分离型解法（Partitionedmethod）求解。

为了将压力方程和运动方程进行联立，优选使用SIMPLE（Semi-Implicit Methodfor Pressure-LinkedEquations（用于压力耦合方程组的半隐式法））算法。

图7显示了通过计算机进行的流动计算的一个例子的流程图。

**步骤s51

在本实施例中，首先，设定压力梯度和速度的上下限，并且计算速度梯度和压力梯度。

**步骤s52

接着，通过当前的压力场的离散化，对运动方程进行定义，并且通过迭代法求解运动方程。即，计算沿三个坐标轴方向的材料模型（或多相）的速度。关于迭代法，可以使用高斯-赛德尔法。

**步骤s53

接着，为了检验上述速度是否满足质量守恒方程，首先计算在腔室模型5的单元表面处的材料模型的未修正的质量流量。此处，所述“未修正的质量流量”是暂时用于SIMPLE算法的循环开始时的质量流量。这样的质量流量可能会有大的误差，因此称其为“未修正的质量流量”。

所述质量流量是确定是否满足质量守恒定律所必需的。

**步骤s54

接着，使用SIMPLE算法，将速度场和压力场联立，并且产生用于修正压力场的下述压力修正方程。

▽[k▽φ]＝src

**步骤s55

接着，通过迭代法例如AMG求解器、共轭梯度法（CG）或双共轭梯度法（Bi-CG）等来求解压力修正方程，并且计算压力修正量p′。

**步骤s56

接着，基于得到的解，通过下述公式对压力场进行修正：

pⁿ⁺¹=pⁿ+ωp′

其中，

p：压力、

n：当前的时间步进数、

ω：松弛系数，在本实施例中其为0.3。

**步骤s57

接着，对界面的边界条件进行修正（或更新）。具体地，由修正的压力场得到压力梯度。将得到的压力梯度设置为边界条件。

**步骤s58

接着，通过下述公式对单元表面处的质量流量进行修正：

mfⁿ⁺¹=mf^*+m′f

其中，

mfⁿ⁺¹：修正后的单元表面处的质量流量、

mf^*：未修正的单元表面处的质量流量、

m′f：质量流量的修正值。

**步骤s59

接着，通过下述式对速度场进行修正：

其中，

V：单元的体积、

v^*：从运动方程得到的修正前的中间速度场、

▽p′：压力梯度的修正量、

运动方程的矩阵的对角分量。

**步骤s60

接着，通过求解上述能量方程（5），计算材料模型的温度和粘度。

**步骤s61

进一步，通过下述式（7）计算腔室模型5的壁表面的滑移速度v_slip。

式(7)：

其中,

τ_w：壁表面剪应力、

V_slip：在腔室模型5的表面处的材料模型的滑移速度、

v_wall：在腔室模型5的表面处的材料模型的切向速度、

F_slip：表示滑移倾向的材料特性（系数）。如果其值增大，则变成壁表面非滑移条件。

**步骤s62

接着，判断计算的解是否收敛。预先定义了收敛的判断标准。

**步骤s63

如果判断为收敛（在步骤s62中的Y），则将每个时间间隔调前一个小时，增加一个时间步进。

可以基于修正后的质量流量的总数是否在预定的误差范围内来判断收敛。如果时间没有到达预定时间（步骤s64中的N），那么计算机再次执行步骤s51和后续步骤。

另一方面，在步骤s62中，如果判断为未收敛，那么计算机再次执行步骤s51和后续步骤。

在本实施方式中的流动计算中，计算了材料模型的剪切生热，但是没有计算气相模型的剪切生热。材料模型的剪切速度相对较大，因此，需要计算剪切生热。剪切速度常常在材料模型和气相模型之间的接触面处增加。

剪切生热正比于剪切粘度乘以剪切速度的平方，因此，有可能异常大的剪切速度值会不利地影响温度的计算（能量方程）。

因此，在本实施方式中，剪切生热的计算仅限于材料模型。

更详细地说，仅对腔室模型5的一些单元进行计算，其中，材料模型的体积分数αq的值不小于恒定值（本实施方式中为0.90）。

在本实施方式中，在材料模型的表面与腔室模型5的表面之间的界面上定义壁表面滑移条件。

在本实施方式中，按如下所述一点一点地将滑移速度给予材料模型。

通过计算机对基于当前的压力场而定义的运动方程进行求解，得到剪应力。

由得到的剪应力，通过使用式（7）来计算待设定的壁表面滑移速度（上述步骤s61）。

在下一步中，如果就这样将计算出的滑移速度v_slip给予接触面，则有可能由式（2）给出的壁表面剪应力τw变成异常值，并且由此在壁表面上定义出异常的滑移速度。

因此，如下述式(8)中所示，优选设置松弛系数β，并且通过改变松弛系数β将滑移速度逐渐增加至目标值。

式(8)：

其中，

0<β<1、

n：计算的重复数。

然后，通过使用SIMPLE算法等，计算机进行压力场和速度场联立，给出计算的滑移速度v_slip，并且产生用于修正压力场的压力修正方程（上述步骤s52）。

比较试验

为了确认本发明的效果，基于图3中显示的腔室模型和表1中显示的条件，进行塑性材料的流动模拟。

关于分布、混炼材料的分散和温度以及转子的转矩，将计算精度评价为三个等级，其中等级1表示良好、并且等级3表示差。结果显示于表1。

表1

根据结果，证实了根据本发明，计算精度得以显著改善，同时保持计算时间不变。

图8显示了根据本发明的方法，通过设计的混炼模拟得到的混炼塑性材料的状态的暂时改变，所述状态的暂时改变从显示于表6的初始状态开始。在图8中，实心黑色部分表示塑性材料。

如上所述，根据本实施方式，可精确模拟塑性材料的混炼状态，并且可输出视觉化的模拟结果。

Claims (3)

1.一种用于模拟腔室中塑性材料的流动状态的计算机化方法，包括：

生成作为腔室的有限体积模型的腔室模型的步骤，

定义塑性材料的材料模型的步骤、

定义气相模型的步骤，所述气相模型是存在于腔室中的气体的模型、

对处于腔室模型中的材料模型和气相模型进行流动计算的步骤，

其中，

所述材料模型设定有塑性材料的粘度，并且

所述气相模型的粘度设定在气体实际粘度的5-10倍的范围内，

在进行流动计算的步骤中，使用流体体积法即VOF法将在腔室模型中气相模型和材料模型汇流的两相流视为单相流，所述单相流的粘度和密度通过各相即材料模型和气相模型所占的体积来确定权重，并且进行平均，以及

在流动计算中，仅对腔室模型的一些单元进行剪切生热计算，其中，材料模型的体积分数不小于0.90。

2.如权利要求1的计算机化方法，其特征在于，

腔室是设置有至少一个转子的班伯里密炼机的混炼空间，并且

塑性材料是待混炼的未硫化橡胶或树脂材料。

3.如权利要求1或2的计算机化方法，其特征在于，

在塑性材料模型和腔室模型之间的界面上定义滑移边界条件。