粘弹性材料的模拟方法

摘要

本发明涉及模拟以树脂或橡胶为基体，掺混了填料的粘弹性材料形变的方法。该方法包括的步骤有：将粘弹性材料分割成有限数量的单元以形成粘弹性材料模型的步骤；在设定条件下进行粘弹性材料模型的形变计算的步骤；以及从形变计算中获得所需的物理量的步骤。其中，将粘弹性材料分割成有限数量的单元的步骤包括：将至少一种填料分割成有限数量的单元以形成填料模型的步骤，以及将基体分割成有限数量的单元以形成基体模型的步骤。

说明书

                              技术领域

本发明涉及很精确地模拟用于分析的粘弹性材料的形变等的方法。

                              背景技术

以橡胶为代表的粘弹性材料被广泛用于，如轮胎以及工业用品(如体育用品)。粘弹性材料在施加载荷时，发生很大的形变，而当载荷被完全移走或卸载时，粘弹性材料又恢复到初始的状态。粘弹性材料在静态载荷下具有非线性弹性行为，而在周期性载荷下具备滞后的速率依赖性或粘弹性行为。为了减少实验室制备的困难和费用，人们用电脑来模拟，如粘弹性材料的形变过程。日本公开专利No.2002-365205披露了一种传统的粘弹性材料的模拟方法。

以上提到的专利聚焦于粘弹性材料对应于不同的应变速度表现出不同的纵向弹性模量这一事实。更具体的说，是在测试条件为预先假定粘弹性材料的实际使用状态的前提下，测定了粘弹性材料产生的应变、应变速度和应力。这样，就得到了纵向弹性模量和应变速度之间的对应关系。对于作为分析目标的粘弹性材料的模量，先给出预定的应变速度，再根据上述的对应关系，计算出纵向弹性模量来进行形变大小的运算。

粘弹性材料的模拟方法描述于，如，下面的文章中。

《橡胶弹性材料的高延展性行为的三维结构模型(AThree-dimensional ConstitutiveModel for the Large Stretch Behavior of Rubber Elastic Materials)》，Ellen M.Aruuda和MarryC.Boyce著，《固体力学与固体物理学杂志(Journal of the Mechanics and Physics ofSolids)》，第41卷第2期，第389-412页(1993年2月)。

下面简要的介绍该文章的内容。

上述文章基于分子链网络模型理论，该理论认为粘弹性材料具备微观的网络结构。如图35所示，粘弹性材料“a”的网络结构包括一系列连接于联结点b的分子链c。联结点b包括分子间的化学联结点，如，化学交联点。

每个分子链c都由一系列链段e构成，一个链段e就是最小的重复组成单元。而且，一个链段e由一系列单体f连接组成，单体中，碳原子之间以共价键连接。每个碳原子可绕其间的键轴相互自由转动，因此，链段e可作为一个整体而被弯曲成各种形状。

在上述文章中，相对于原子的波动周期，联结点b的平均位置长时期内不会改变。因此，联结点b的波动可忽略不计。而且，分子链c的终端矢量具备两个联结点b，在两个末端上的联结点b被假定为随着被嵌入在分子链中的粘弹性材料“a”这一连续体而发生形变。

Aruuda等人还提出一个八链段橡胶弹性体模型。如图6所示，这个模型将粘弹性材料的微观结构定义为立方体网状结构体h，其间，聚集了多个微观八链段橡胶弹性体模型g。在单个八链段橡胶弹性体模型g中，分子链c从位于立方体中心的一个联结点b1延伸至每一个位于立方体各个顶点的联结点b2，如图6右边放大部分所示。

模拟时，粘弹性材料被定义为超弹性体，其体积几乎不改变，且在载荷去除后能恢复至初始形状。超弹性体，如下面的方程式(1)所表示的，被定义为是具有应变能量函数W的一种物质。函数W被格林(Green)应变的一个分向量Eij微分，得到共轭基尔霍夫(Kirchhoff)应力S_ij。换言之，应变能量函数表明粘弹性材料形变时其存储的势能的存在。因此，超弹性体的应力和应变之间的关系可由应变能量函数W的微分斜率得到。

$$ >>>S>ij>>=>>>\partial >W>>>>\partial >E>>ij>>>·>·>·>·>·>·>>(>1>)>>>$$

基于非高斯统计理论，Aruuda等人认识到，当粘弹性材料的形变增加时，其熵变也显著增加(分子链伸展并取向)，此时，橡胶弹性体的应变能量函数W如方程式(2)所示。而且，通过将应变能量函数W代入上述方程式(1)后，就可得到弹性材料的应力和应变之间的关系。

$$ >>W>=>>C>R>>N>>(>>>>I>1>>>3>N>>>>·>\beta >+>ln>>\beta >>sinh>\beta >>>)>>->T>·>n>·>c>·>·>·>·>·>·>>(>2>)>>>$$

C^R＝n·k_B·T

(n：每个体积单元内所含分子链的数目；k_B：波尔兹曼常数；T：绝对温度)

I₁：应变初始不变量，I₁＝λ₁²+λ₂²+λ₃²(λ₁、λ₂、λ₃为伸长率)

N：每个分子链的平均链段数目

$$ >>\beta >=>>L>>->1>>>>(>>>r>chain>>>N>·>a>>>)>>>$$

(r_chain：单个分子链的末端间的距离，

$$ >>>>>N>·>>I>1>>>3>>>·>a>;>>$$a：链段的长度；L：朗之万(Langevin)函数，

$$ >>L>>(>x>)>>=>coth>x>->>1>x>>)>>$$

应用Aruuda等人给出的应力和应变的关系，对粘弹性材料进行单轴拉伸形变的模拟，可得到应力和应变间的非线性关系，例如，如图36所示。其结果与加载形变时的实际测量值非常相关。

在粘弹性材料被用作工业用品时，通常掺混填料(填充剂)，如碳黑和二氧化硅。为了能很精确地模拟掺混了填料的粘弹性材料的形变，就不能忽略填料的存在。然而，在传统的模拟方法中，填料都未被考虑在内。

各种试验的结果表明，在填料和基体的界面上，会产生各种现象，例如填料和基体的滑移和摩擦，而且在这些区域，会产生较大的能量损失。因此，为了精确地进行粘弹性材料的模拟，计算中将上述发生在界面上的现象加以考虑就显得非常重要。

                              发明内容

考虑到上述的各种缺点，本发明的目的在于提供一种粘弹性材料的模拟方法，该方法能很精确的模拟粘弹性材料模型的形变状态，该模拟则基于两个步骤，一个步骤是将粘弹性材料分割成有限数量的单元，另一个步骤是将至少一种填料分割成有限数量的单元以形成填料模型。

本发明提供了一种模拟基体为橡胶或树脂，掺混有填料的粘弹性材料形变的方法。该方法包括以下步骤：将粘弹性材料分割成有限数量的单元以形成粘弹性材料模型的步骤；在预先设定的条件下，进行粘弹性材料模型形变计算的步骤；以及从形变计算中获取必要物理量的步骤。其中，将粘弹性材料分割成有限数量的单元的步骤包括将至少一种填料分割成有限数量的单元以形成填料模型的步骤以及将基体分割成有限数量的单元以形成基体模型的步骤。

本发明包括一个将填料分割形成填料模型的步骤，因此，填料的影响被考虑在粘弹性材料模型的形变计算结果里面了。

将粘弹性材料分割成有限数量的单元的步骤还宜进一步包括在基体模型和填料模型之间配置具有与基体模型不同粘弹性能的界面模型的步骤。此界面模型可以具有比基体模型更软的粘弹性能，或者是比基体模型更硬的粘弹性能，以及滞后损失比基体模型大的粘弹性能。在基体模型中，应力和应变之间的关系宜被限定于模拟基体橡胶的能量损失。

填料模型宜包括至少两个填料模型：第一填料模型和第二填料模型，两个模型设置成彼此之间有一段距离。在第一模型和第二模型之间，设置有填料间模型，其填料间吸引力随距离而变。在粘弹性材料模型的形变计算中，将填料间相互作用考虑在内是十分有用的。

所述填料间引力宜设定设定在抛物线函数上，其包括随第一填料模型和第二填料模型间距的增大，抛物线有一逐渐上升的区域，曲线平滑上升，到达一个顶点；随间距的进一步增大，曲线有一逐渐下降的区域，曲线平滑下降，在一个预先设定的特征长度上，到达零点。而且，在曲线逐渐下降区域的卸载变形中，填料间引力从卸载初始值线性递减至零值这一步骤也包括在内，以将第一填料模型和第二填料模型间产生的滞后损失在模拟时考虑进去。

                              附图说明

本发明，及其目的和优点可在参照下面优选的实施方式的描述，并结合附图得以最好理解。其中：

图1为本发明所用的电脑设备的一个具体实例的透视图；

图2为本发明优选实施方式的程序流程图；

图3为粘弹性材料模型(微观结构)的优选实施方式的示意图；

图4为建成粘弹性材料模型(微观结构)的步骤的流程图；

图5为分子链联结点断裂示意图；

图6为粘弹性材料网状结构的一个实例及八链段橡胶弹性体模型的透视图；

图7为参数λ_c和每一分子链所含平均链段数目N的关系图；

图8为碳黑形状示意图；

图9为用均化法表示的完整结构和微观结构之间关系的视图；

图10为形变模拟程序的流程图；

图11为基体模型、填料模型以及界面模型的各自的应力-应变曲线图；

图12为粘弹性材料模型(完整结构)的模拟出的实际的应力-应变关系曲线图；

图13为另一种实例的基体模型和界面模型的应力-应变关系曲线图；

图14(A)至图14(J)为示意图，直观地显示了粘弹性材料模型的形变过程；

图15(A)为未将界面层考虑在内的模型的拉伸应变分布示意图；

图15(B)为将界面层考虑在内的模型的拉伸应变分布示意图；

图16(A)为未将界面层考虑在内的模型的能量损失分布示意图；

图16(B)为将界面层考虑在内的模型的能量损失分布示意图；

图17为掺混有不同含量的碳黑的橡胶材料的拉伸测试结果示意图

图18为粘弹性材料模型(微观结构)的另一实施方式的示意图；

图19为图18主要部分的放大图；

图20(A)为界面的截面视图；

图20(B)为引力T_n和界面间距d的关系曲线图；

图21为界面间的放大视图；

图22为填料模型间，引力T_n和间距d的关系曲线图；

图23为填料模型间，加载形变和未加载形变的示意图；

图24为完整的粘弹性材料模型模拟出的实际的应力-应变关系曲线图；

图25为各填料模型间产生的引力和间距d的关系曲线图；

图26为另一模拟实施方式的流程图；

图27为一系列碳黑二次粒子的放大示意图；

图28为碳黑初级粒子建模示意图；

图29为填料模型的放大示意图；

图30(A)为填料模型基本结构的透视图；

图30(B)为组成基本结构的网状物排列的平面图；

图31(A)和图31(B)为设定设定填料模型实例的示意图；

图31(C)为填料模型横截面的示意图；

图32为计算步骤的过程的一个具体实例的流程图；

图33为填料模型间的间距和势能总和的关系曲线图；

图34为填料模型间的间距和填料模型作用力的关系曲线图；

图35为粘弹性材料的放大示意图，包括了一个分子链结构，一个放大的分子链，以及一个放大的链段；以及

图36为传统模拟方法模拟出的应力-应变曲线图；

                            具体实施方式

下面，将参照附图，对本发明的优选实施方式进行描述。图1是本发明用以进行模拟的计算机设备1的示意图。计算机设备1包括主体单元1a、键盘1b、用作输入设备的鼠标1c、以及用作输出设备的显示器1d。尽管图上没有显示，但主体单元1a应包括中央处理器(简称CPU)、只读存储器(ROM)、工作存储器、一个大容量存储设备如硬盘，以及CD-ROM或软盘的驱动器1a1和1a2。大容量存储设备用于存储执行模拟时的过程步骤(即程序)，这在后面描述。计算机设备1优选工程工作站之类。

图2为本发明的模拟程序的一个实例。在此实例中，首先设定粘弹性材料的模型(步骤S1)。图3中，则直观地显示以微观结构呈现的粘弹性材料模型2的一个实例。

为设定粘弹性材料模型2，将用于分析的粘弹性材料的微观区域(无论其是否存在)分割成有限的小单元2a、2b、2c……。每个单元2a、2b、2c……都赋予用统计分析法进行形变计算所必需用到的参数。统计分析法包括，如，有限数量的单元法、有限体积法、有限微分法、边界元法等。参数包括，如，联结点座标值、单元形状、和/或单元2a、2b、2c……的材料性能。而且，单元2a、2b、2c……中，用三角形单元或四边形单元作为二维平面，用四面体单元或六边形单元作为三维单元。这样，粘弹性材料模型2就成为电脑设备1可用的数字数据。在此例中，所展示的是粘弹性材料的二维模型。

图4则显示了设定粘弹性材料模型2的详细步骤。在本实施方式中，粘弹性材料模型2由如下几个步骤设定：将基体分割成有限数量的单元以形成基体模型3的步骤S11；将至少一种填料分割成有限数量的单元以形成填料模型4的步骤S12；以及形成与基体模型3有不同粘弹性能的界面模型5的步骤S13。

因此，图示的粘弹性材料模型2由基体模型3、填料模型4(白色部分)、以及界面模型5(浅黑色部分)所组成。基体模型3中，基体由橡胶或树脂组成并被模型化；填料模型4则分布并掺混于基体模型3中，其中，填料被模型化了。界面模型5则位于基体模型3和填料模型4之间，以在模型间形成界面。

图中，基体模型3为黑色部分。本例中，基体模型3构成了粘弹性材料模型2的主体，并被分割成三角形元或四边形元。如下列关于材料性质的方程式(3)所示的应力和应变的关系被用于基体模型3的每个单元。

$$ >>ver>>S>▿>>ij>>=>-ver>>p>·>>>\delta >ij>>+>>R>ijkl>>ver>>ϵ>·>>kl>>·>·>·>·>·>·>>(>3>)>>>$$

其中，

$$ >>>R>ijkl>>=>>1>3>>>C>R>>>N>>\{>>(>>>\xi >c>>>N>>>->>>\beta >c>>>\lambda >c>>>)>>>>>B>ij>>>B>kl>>>>B>mm>>>+>>(>>>\beta >c>>>\lambda >c>>>>(>>\delta >ik>>>B>jl>>+>>B>ik>>>\delta >jl>>)>>)>>\}>>$$

C^R＝n·k_B·T

(n：每单位体积内所含分子链的数目；k_B：波尔兹曼常数；T：绝对温度)N：每个分子链所含的平均链段数目

$$ >>>\xi >c>>=>sup>>\beta >c>2sup>>>1>-sup>>\beta >c>2sup>>csc>>h>2>>>\beta >c>>>>,ver>>p>·>>=>->K>ver>>ϵ>·>>mm>>,>K>=>\psi >max>>(>>R>ijkl>>)>>>$$

$$ >>>\beta >c>>=>>L>>->1>>>>(>>>\lambda >c>>>N>>>)>>,sup>>\lambda >c>2sup>>=>>1>3>>>(sup>>\lambda >1>2sup>>+sup>>\lambda >2>2sup>>+sup>>\lambda >3>2sup>>)>>>$$

λ₁²+λ₂²+λ₃²＝I₁，I₁：应变初始不变量

B_ij：剩余柯西-格林(Cauchy-Green)形变张量

现在简要介绍方程式(3)的推导过程。粘弹性材料如橡胶在形变过程中体积的变化非常小，在计算时可以忽略不计。因此，在计算时粘弹性材料的密度被作为一个常数。这样，基尔霍夫应力S_ij可由方程式(4)来表示。方程式(4)中，E_ij是格林应变的一个分向量，p是液静压，X_i为应力应变均为零的C0状态下，任意目标点P的位置，X_j为形变状态C下，目标点P的位置。

$$ >>>S>ij>>=>>>\partial >W>>>>\partial >E>>ij>>>->P>>(>>>>\partial >X>>i>>>>\partial >x>>m>>>)>>>(>>>>\partial >X>>j>>>>\partial >x>>m>>>)>>·>·>·>·>·>·>>(>4>)>>>$$

下面方程式(5)中表达了柯西应力分向量σ_ij和基尔霍夫应力分向量S_mn之间的关系，因此，可从方程式(4)推得方程式(6)。

$$ >>>\sigma >ij>>=>>1>J>>>>>\partial >x>>i>>>>\partial >X>>m>>>>S>mn>>>>>\partial >x>>j>>>>\partial >X>>n>>>·>·>·>·>·>·>>(>5>)>>>$$

$$ >>>\sigma >ij>>=>>>\partial >W>>>>\partial >E>>mn>>>>>>\partial >x>>i>>>>\partial >X>>m>>>>>>\partial >x>>j>>>>\partial >X>>n>>>->>p\delta >ij>>·>·>·>·>·>·>>(>6>)>>>$$

在方程(5)中，如果体积为常数，则代表目的物体积变化速率的“J”就被定为“1”。而且，由于方程式(2)的应变能量函数仅仅是剩余柯西-格林形变张量Aij的应变初始不变量I₁的函数，因此，又可从方程(6)推导得方程式(7)

$$ >>>\sigma >ij>>=>2>>>\partial >W>>>>\partial >I>>1>>>>A>ij>>->>p\delta >ij>>·>·>·>·>·>·>>(>7>)>>>$$

而且，由下列关系式(8)、(9)(10)，方程式(7)所显示的粘度形式又可由方程式(11)来表达。

$$ >ver>>\sigma >▿>>=ver>>\sigma >·>>->W\sigma >+>\sigma W>>$$

$$ >>ver>>I>·>>1>>=>>(>trver>>A>·>>)>>=>2>A>·>D>·>·>·>·>·>·>>(>8>)>>,>>(>9>)>>>(>10>)>>>$$

$$ >ver>>A>·>>=>>(>D>+>W>)>>A>+>A>>(>D>->W>)>>>$$

$$ >>ver>>\sigma >▿>>ij>>=>-ver>>p>·>>>\delta >ij>>+>>1>3>>>C>R>>>N>>\{>>(>>>\xi >c>>>N>>>->>>\beta >c>>>\lambda >c>>>)>>>>>B>ij>>>B>kl>>>>B>mm>>>+>>>\beta >c>>>\lambda >c>>>>(>>\delta >ik>>>B>jl>>+>>B>ik>>>\delta >jl>>)>>\}>·>·>·>·>·>·>>(>11>)>>>$$

在体积为常数的条件下，方程式(11)所表示的柯西应力的Jaumann粘度可以被基尔霍夫应力的Jaumann粘度替换。而且，通过用应变粘度张量替换形变粘度张量D，即得粘弹性不可压缩材料粘度的表达形式的基本方程式，即方程式(3)。

然而，如果仅仅使用方程式(3)，则与Aruuda等人的文章一样，不能将能量损失这一粘弹性材料的特征在模拟时考虑进去。就是说，当模拟从图36中点Q移走负载时，应变恢复曲线与加载时的曲线大致相同。这样，就有滞后回线没有体现这一缺点。

发明者尝试对Aruuda等人提出的模型进行各种改进。如上所述，粘弹性材料可以承受非常大的应变，通过彼此杂乱纠结的大分子链C的拉伸，其应变可以达到百分之几百。发明者假设：加载形变过程中，粘弹性材料中相互纠结的分子链C会有部分变得松弛并消失(联结点b的数量减少)，或者在移走负载时，分子链间纠结可能会再次产生(联结点b的数量增加)。

如图5所示，当拉伸应力Y施加到联结在联结点b的分子链c1、c2、c3、c4上时，c1到c4的每个分子链都被拉长，联结点b受到很大的拉力，可能会发生断裂(消失)。如图5右边所示，分子链c1和c2合并成了一个长分子链c5，c3和c4也是如此。这种现象在橡胶材料的加载形变过程中相继发生，并可能引起大量的能量损失。

为了模拟基体橡胶的能量损失，将上述现象优选应用到方程(2)中。为达到此目的，每个分子链c中所含的平均链段数N被定义为可变参数，在加载形变和卸载形变中是各不相同的。在这里，加载形变指的是在极短的时间内，基体模型3的应变的增加，卸载形变指应变的减少。

图6中的微观三维网络结构体h再一次被用以解释上述观点。网络结构体是这样的体系：在其轴向上、高度方向上、深度方向上各结合有k个八链段橡胶弹性体模型。需要指出的是，k是个足够大的数字。包括在相关的网络结构体h中的联结点b的数目被称为“联结数”。假设联结数为“m”，网络结构体h中的分子链c的数目(即基体模型3的单位体积内的分子链数目)为“n”，m和n可分别表达为

m＝(k+1)³+k³    ……(12)

n＝8k³          ……(13)

由于k是个足够大的数字，当k的三次方以外的项被省略后，上面的两个方程式可分别由下面的方程式(14)、(15)来表达：

m＝2k³            ……(14)

n＝8k³            ……(15)

由方程式(14)和(15)的关系，可知联结数m可由n来表达，如方程式(16)。

m＝n/4            ……(16)

而且，由于甚至基体模型3发生变形时，模型3中所含的分子链所含的链段总数N_A也不改变，因此，可得方程式(17)和(18)。

N_A＝n·N           ……(17)

N＝N_A/n＝N_A/4m    ……(18)

由方程式(18)可知，在加载形变和卸载形变(应变的恢复过程)时，橡胶分子链的联结数m应改变以使平均链段数目N为其可变参数。这样，以平均链段数目N为在加载形变和卸载形变中可取不同值的可变参数，就可以模拟基体模型3的能量损失。

平均链段数目N可以有多种方法确定。例如，平均链段数目N可在加载形变过程中应变的相关参数的基础上增加。与应变相关的参数没有特别的限制，可以是，比如，应变、应变速度、或者应变初始不变量I₁。在本发明的实施方式中，平均链段数目N由下面方程式(19)给出定义。这个方程式表明，平均链段数目N是相应基体模型3单元中应变初始不变量I₁(更具体的说，是参数λ_c平方根)的函数。

N(λc)＝A+B·λc+C·λc²+D·λc³+E·λc⁴    ……(19)

其中，$$ >>>\lambda >c>>=>>>>I>1>>3>>>>$$

方程式(19)是经多次试验设定起来的，根据如，粘弹性材料的单轴拉伸测试结果，可以很容易的确定常数A～E的具体值。例如，可以首先得到粘弹性材料或是被分析的物体的应力-应变曲线，然后，在负载移去时，曲线中引入n和N。这样，就可以确定分子链中链段的总数N_A(＝n·N)。在这里，由于加载形变过程和卸载形变过程中，分子链c中链段总数N_A是相同的，根据加载时的曲线推出平均链段数N。而根据加载时的平均链段数N，可以确定方程式(19)中的参数A～E。本实施方式中，N＝6.6，以上常数依次设定为：

A＝+2.9493

B＝-5.8029

C＝+5.5220

D＝-1.3582

E＝+0.1325

图7中，展示了基体模型3的单元加载形变过程中平均链段数N和参数λ_c之间的关系。当与应变相关的参数λ_c增大时，平均链段数N也逐渐增大。在此实例中，参数λ_c的上限是2.5。在下文中要描述的粘弹性材料模型2的形变模拟中，在发生加载形变时，基于稳定的基础上计算基体模型3的每个单元的参数λ_c。将λ_c的计算值代入方程式(19)，即可计算出相关单元的相关应变状态下链段平均数N。这里，基体模型3的卸载形变过程中的链段平均数N的值是常数。

在本实施例中，填料模型4中，碳黑也被模型化了，需要指出的是，填料不限于碳黑，它还可以是，如，二氧化硅等。本实施例中，填料模型4的物理形状是设定在用电子显微镜观察到的填充在实际的橡胶中的碳黑的形状的基础上。图8为碳黑6的二次粒子的示意图。更具体的说，二次粒子具有这样的结构：每个由碳原子组成的、直径大约在10nm左右的一系列球形初级粒子7无规则的键合成三维的结构。

碳黑6的硬度(纵向弹性模量)比基体橡胶高数百倍。因此，本实施方式中，填料模型4被定义为弹性体而不是粘弹性体。因此，填料模型4就以纵向弹性模量作为材料性能，且在形变计算中，其应力和应变也是成比例的。填料模型4的数量应取决于所分析物体的粘弹性材料中掺入的填料数量。

界面模型5位于基体模型3和填料模型4之间。界面模型5无需限定于连续的围绕在填料模型4周围，但最好是在整个范围内围绕在填料模型4周围。在此例中，界面模型5的厚度非常小。界面模型5的厚度t在1～20nm之间，优选在5～10nm之间。

填料和橡胶基体的界面上，具备如此厚度的界面层的物理结构实际上是观察不到的。然而，填料和橡胶基体的界面上发生的各种导致如滑移和摩擦等能量损失的现象是可以观察到的。从而，为了在模拟时，将这些现象考虑在内，模拟计算时就将界面模型定义为粘弹性材料。因此，与基体模型3类似的，方程式(3)所表达的应力和应变之间的关系也可应用于界面模型5。

需要指出的是，界面模型5与基体模型3的粘弹性能并不相同。本实施例中的界面模型5比基体模型3的粘弹性更软。因此，当在界面模型5和基体模型3上施加同样的应力时，界面模型5的应变比基体模型3的应变要大很多。而且，界面模型5的粘弹性的滞后损失(每一个应变循环产生的能量损失)比基体模型3的要大。这可用于，如，调整方程式(19)的参数以计算方程式(3)所将用的链段平均数N。

下面，设定粘弹性材料模型2的形变条件。(步骤S3)

本实施例中，粘弹性材料模型2在单轴向拉伸的条件(平面应变状态)下发生形变。因此，粘弹性材料模型2在如图3所示Z轴方向上不发生应变。形变条件包括粘弹性材料模型2形变时的应变速度和最大应变。而且，形变计算应在如图2所示的粘弹性材料模型2的微观片段上进行，但优选在如图9所示的、粘弹性材料模型2(如图3所示)的微观结构在垂直方向和水平方向上周期性重复而组成的粘弹性材料完整模型M上进行。

当进行如上所述的粘弹性材料完整模型M的形变计算时，比较合适的方法是均化法，在粘弹性材料完整模型M中，微观结构的重复程度非常接近。因此，仅用有限数量的单元法很难直接将完整的模型M分割。而在均化法中，用到了如图9所示的代表粘弹性材料完整模型M微观标尺X₁，X₂和代表上述微观结构的微观标尺Y1、Y2的两个独立变量。在均化法中，不同的微观标尺Y1、Y2和微观标尺X₁，X₂中的各独立变量为渐进展升，因此，可获得具备一定尺寸的、周期性含有如图3所示的微观结构的模型的粘弹性材料完整模型M的平均动力学反应特性。

渐进展开均化法已经在数值计算法中被确定，这种方法被详细的描述于，如，下面的文章中。

Higa，Y.和Tomita，Y.，Computational Prediction of Mechanical properties of Nickel-basedsuperalloy with gamma Prime Phase Precipitates，Proceedings of ICM8(Victoria，B.C.，Canada)，Advance Materials and Modeling of Mechanical Behavior，(Ellyin.F.，和Prove.J.W主编)，III，1999年，第1061-1066页，Fleming印刷出版社。

粘弹性材料模型2的一个面的尺寸为300nm×300nm，而图9所示的粘弹性材料完整模型M的面为2mm×2mm的长方形。将恒定的应变速度施加于图3所示的微观结构，以得到均匀的单轴拉伸应变E₂(其应变速度为1.0×10^-5/s)。当图9中X₂方向上的应变E₂达到0.65，即设定起在相同的应变速度下将应变逐渐降为零的条件。

形变计算(模拟)用设定的模型和条件来进行(步骤S3)。图10中，展示了形变计算的一个具体程序的实例。形变计算时，数据首先输入电脑设备1(步骤S31)。输入的数据包括，如，构成粘弹性材料模型2的数值数据，以及预设的一系列边界条件。

随后，形成各单元的刚性基体(步骤S32)，接着，组合整个结构中的刚性基体(步骤S33)。再对完整结构的刚性基体导入已知的节点位移和节点力(步骤S34)，并对刚性方程式进行分析。接着确定未知的节点位移(步骤S35)，并计算和输出如应变、应力、以及各单元的主应力等物理量(步骤S36，S37)。步骤S38中，确定计算是否完成，如果计算没有完成，则重复步骤S32以后的步骤。在有效工作的原则基础上，方程式(20)用作为单元的方程式。

$$ >\underset{}{}$$

其中，为基尔霍夫应力速度，σ_mj为柯西应力，v_i，m为位移速度曲线斜率，为表面作用力的速度，Φ为单位表面积内的不可逆功。

在均化法的计算中，要用到下面的微观平衡方程式(21)和性能位移函数(Y-周期)方程式(22)

$$ >sup>>L>ijkl>Hsup>>=>>1>>|>Y>|>>>>\int >Y>>[>>L>ijkl>>->>L>ijkl>>>1>2>>>(>sup>>>\partial >\chi >>p>kjsup>>>>\partial >y>>q>>>+>sup>>>\partial >\chi >>q>kjsup>>>>\partial >y>>p>>>)>>]>dY>>$$

$$ >sup>>\tau >ijkl>Hsup>>=>>1>>|>Y>|>>>>\int >Y>>>(>>\sigma >ij>>>\delta >kj>>->>\sigma >mj>>sup>>>\partial >x>>j>klsup>>>>\partial >y>>m>>>)>>dY>·>·>·>·>·>·>>(>21>)>>>$$

$$ >>>\int >Y>>[>>L>ijpq>>>1>2>>>(>sup>>>\partial >\chi >>p>klsup>>>>\partial >y>>q>>>+>sup>>>\partial >\chi >>q>klsup>>>>\partial >y>>p>>>)>>+>>\sigma >qj>>>\delta >pi>>sup>>>\partial >\chi >>p>klsup>>>>\partial >y>>q>>>]>>>>\partial >\delta v>>i>>>>\partial >y>>j>>>dY>=>>\int >Y>>>(>>L>ijkl>>+>>\sigma >ij>>>\delta >ki>>)>>>>>\partial >\delta v>>i>>>>\partial >\delta y>>j>>>dY>>$$

x^kl……Y-周期

$$ >>>\int >Y>>[>>L>ijkl>>>1>2>>>(>>>>\partial >\phi >>k>>>>\partial >y>>l>>>+>>>>\partial >\phi >>l>>>>\partial >y>>k>>>)>>+>>\sigma >mj>>>>>\partial >\phi >>i>>>>\partial >y>>m>>>]>>>>\partial >\delta v>>i>>>>\partial >y>>j>>>dY>=>>\int >Y>sup>>P>ij>\prime sup>>>>>\partial >\delta v>>i>>>>\partial >\delta y>>j>>>dY>>$$

φ……Y-周期                                    ……(22)

形变的计算，可通过有限数量的单元法，用工程系统分析应用软件(如，LS-DYNA等，LivermoreSoftware Technology Corporation(US)研发并改进)来进行。

本实施方式中，方程式(3)的常数设定如下：

C_R＝0.268

N＝6.6

T＝296

K_B＝1.38066×10^-29

n＝6.558×10²⁵

N_A＝4.328×10²⁶

填料模型的体积容量：30％

填料模型的纵向弹性模量E：100Mpa

填料模型的泊松比v：0.3

形变计算时，如上所述，计算出每个应变状态的平均链段数，将其值代入方程式(3)随后进行的计算。在粘弹性材料模型2和粘弹性材料完整模型M中，都用到Aruuda等人提出的三维八链段橡胶弹性体模型，其厚度方向(图3所示的Z轴方向)也未加改变。而且，基体模型3和界面模型5的平均链段数N设定如下：

<基体模型>

加载形变过程中平均链段数N

N＝-3.2368+20.6175λ_c-21.8168λ_c²+10.8227λ_c³-1.9003λ_c⁴

卸载形变过程中平均链段数N(常数)

N＝6.6

分子链的链段总数N_A(常数)

N_A＝4.3281×10²⁶

<界面模型>

加载形变过程中平均链段数N

N＝-5.9286+20.6175λ_c-21.8168λ_c²+10.8227λ_c³-1.9003λ_c⁴

卸载形变过程中平均链段数N(常数)

N＝3.91

分子链的链段总数N_A(常数)

N_A＝3.203×10²⁵

形变计算时，在每个恒定的时间增量里面，都计算应变、应力以及单元的其它此类变量。随后记录其数值(步骤S37)。这样，就得到了必需的物理量，并将其用于分析(图2所示的步骤S4)。

图11中，展示了基体模型3、填料模型4、以及界面模型5的独立的形变计算结果。填料模型4显示出最高弹性，且没有能量损失。与基体模型3相比，界面模型5显示出较软的粘弹性能，且有更大的能量损失(曲线形成的环状表面积)。

图12中，展示了粘弹性材料完整模型M的应变和实际应力之间的关系。曲线La表明了界面模型5比基体模型3更软时的关系结果。曲线La中，第一条非线性曲线La1是在加载形变过程中得到的，第二条曲线La2则不同于第一条曲线La1，它是在卸载形变过程中得到的。第二条曲线La2比第一条曲线La1更软(弹性更低)，且产生了滞后回线。通过计算此回线内的封闭表面积，可得到一个拉伸形变循环内的能量损失。另一方面，在每个分子链的平均链段数为常数的前提下，进行与Aruuda等人相同的形变计算时，加载形变和卸载形变的应力-应变曲线形状是相同的，如图36所示，没有滞后回线产生。

通过适当的改变方程式(19)中系数A～E的数值，可以得到各种不同形状的第一曲线La1和第二曲线La2。因此，通过根据被分析的粘弹性材料，设定系数A～E，可以精确的研究出各材料的能量损失及此类物理量。这对于以粘弹性材料为主体的工业品如轮胎、高尔夫球等的性能的改进是非常有用的。

而且，图12中，给出了将与基体模型3相同的材料性能赋予界面模型5(也就是说，该模型未将界面考虑在内)后的模型计算结果，如由长短破折号交替组成的曲线Lb所示。在这个模型中，由于不存在软质的应变剧烈的界面区域，故与曲线La相比，其应力-应变曲线的斜率上升。而且，由于界面上不发生如界面模型5那样的较大的能量损失，所以整个模型的能量损失较小。

而且，图12中，给出了具有与曲线La相反粘弹性能的模型的计算结果，如由虚线组成的曲线Lc所示。在此模型中，如图13所示，界面模型5具有比基体模型3更硬的粘弹性能。正如预期的那样，从图中也可看出曲线Lc具有高弹性。在此实施方式中，各模型的平均链段数N设定如下：

<基体模型>

加载形变过程中平均链段数N

N＝-3.2368+20.6175λ_c-21.8168λ_c²+10.8227λ_c³-1.9003λ_c⁴

卸载形变过程中平均链段数N(常数)

N＝6.6

分子链的链段总数N_A(常数)

N_A＝4.3281×10²⁶

<界面模型>

加载形变过程中平均链段数N

N＝-5.4800+20.6175λ_c-21.8168λ_c²+10.8227λ_c³-1.9003λ_c⁴

卸载形变过程中平均链段数N(常数)

N＝4.36

分子链的链段总数N_A(常数)

N_A＝8.569×10²⁶

在图14中，直观的展示了粘弹性材料模型2在拉伸形变模拟(曲线La)时，单个微观结构(单位晶格)从形变状态过渡到恢复状态的顺序过程。需要指出的是，在形变计算时此过程不是连续的，它具有一定的时间间隔。图14(A)至图14(E)表明了加载形变过程，图14(F)至图14(J)表明了卸载形变过程。在图14中，应力的级别以颜色的变化来表示。颜色变白的区域表明此区域有较大的应变。此结果正确的表明，较大的应变集中在填料模型4的界面之间或者填料模型4，4之间。尤其是，可以看出，在填料模型4，4间距很小的区域有较大的应变发生。

图15中，直观的展示了在各单元的拉伸方向上，最大应变下的应变分布。其中，图15(A)与图12的未将界面考虑在内的曲线Lb相对应；图15(B)与图12的将界面考虑在内的曲线La相对应。图上显示，较大的应变在颜色较浅的区域发生。未将界面考虑在内的图15(A)中，发生的应变较小，且发生在整个很广泛范围；而将界面考虑在内的图15(B)中，较大的应变集中在填料模型4的周边(或者是界面附近)。

图16中，用颜色显示了各单元一个应变循环的能量损失及其大小。图中表明一个循环的能量损失更多的发生在颜色较浅的区域。模型的形状显示为图15的最大形变形状。图16(A)与图12的曲线Lb相对应；图16(B)与图12的曲线La相对应。图16(B)将界面考虑在内，从图中可以看出，较大的能量损失集中在填料模型4的界面上。

现在将解释本发明的另一实施方式。本实施方式用于分析的物体也是由橡胶掺混碳黑作为粘弹性材料的。图17中，显示了掺混不同含量CB的碳黑的橡胶材料的应力-应变曲线。碳黑掺混量CB大的橡胶材料，具有较大的滞后回线的表面积。因此，填料掺混量大的橡胶材料趋向于具有大的能量损失。其原因被认为可能是因为填料颗粒之间的相互作用。也就是说，在橡胶基体中，与范德华力类似的填料间引力被认为在间距为几个纳米的两种填料颗粒之间作用。这个理论目前比较权威。在本实施方式中，形变计算时，做了将填料间引力考虑在内的尝试。

图18中，直观的显示了以微观结构呈现的粘弹性材料模型2的一个实例。本例中的粘弹性材料模型2包括：基体橡胶被分割的基体模型3；填料模型4A和4B(统称“填料模型4”)，其中，至少有两种其间有空隙的填料颗粒被分割；一个填料间模型10，它位于填料模型4A和4B之间，并具有随填料模型4A和4B间距而变的填充间引力。这种粘弹性材料模型不含界面模型。为了提高计算效率，设定粘弹性材料模型2的微观结构(单元晶格)，以使其在X轴和Y轴的中心线上对称。

计算时，基体模型3和填料模型4按与前面所述的实施方式相同的方法定义。在本例中，两个填料模型4均为椭圆形，且大小相同。填料模型4A、4B以最小距离d(d≠0)隔开，彼此之间互不接触。填料模型4A和填料模型4B之间的初始最小距离d宜选择的范围为，如，约1～3nm。需要指出的是，距离值可按填料模型4的位置、模拟条件等而改变。

如图18，以及放大图19所示，填料间模型10位于填料模型4A和填料模型4B之间。因此，填料间模型10和基体模型3在Z轴方向上重叠。在此例中，填料间模型10均为四方元。在一个填料间模型10中，有两个顶点与填料模型4A的外周边上的顶点共享，剩余的两个顶点与另一个填料模型4B的外周边上的顶点共享。因此，在形变计算时，如果填料模型4A和4B之间的间距d改变了，则填料间模型10各个顶点的座标也就随之改变，所以其形状(表面积)也就改变了。

填料间引力由填料间模型5确定。例如，填料间引力可根据下面的文章来确定，该文章中，明确的描述了一种当相II粒子从金属基体层键解时作用于界面上的粘附力。

《通过掺杂物键解进行中空成核的连续模型》(A Continuum Model for Void Nucleationby Inclusion Debongding)，A.Needleman著，《应用力学杂志》(Joumal of AppliedMechanics)1987年9月，第54卷，第525-531页。

图20(A)为Needleman等人提出的，界面A、B开始键解的示意图。在本实施方式中，界面A、B分别对应于填料模型4A和4B的外周边表面。图20(B)中，显示了A、B界面间产生的引力T_n和间距d的关系。图20(B)的水平轴上的δ是当界面A、B的间距d增加至引力变为零时的极限长度。这个长度在后面称为特征长度。界面A、B之间的引力Tn由一个抛物线函数决定，该抛物线包括：一个逐渐上升的区域，此区域内，随着间距d的增加，引力平稳的增大，并到达一个顶峰(σ_max)；一个逐渐下降的区域，此区域内，随着间距d的进一步增加，引力平稳下降，在预设的特征长度上到达零值。函数被表达为二维的方程式(23)。

$$ >>>T>n>>=>>27>4>>>\sigma >max>>\{>>>u>n>>>\delta >n>>>[>1>->2>>>>(>>>u>n>>>\delta >n>>>)>>2>>+>>>(>>>u>s>>>\delta >s>>>)>>2>>>+>>>(>>>u>n>>>\delta >n>>>)>>2>>+>>>(>>>u>s>>>\delta >s>>>)>>2>>]>\}>·>·>·>·>·>·>>(>23>)>>>$$

这里，u_n，u_s分别为界面A、B中法线和切线的位移，δn，δs为法线和切线的特征长度。法线n和切线s如图21所示。

本实施方式中，填料间引力由上述方程式(23)确定。如果对于粘弹性材料模型2，假设填料模型4A和4B之间，初始状态的间距d不等于零，而填料间引力为零，则对方程式(23)作如下的校正。

首先，为了调整填料模型4A和4B的初始间距d，图20(B)中的曲线在水平轴上平移-d(填料模型的初始间距)，在垂直轴上平移-Td(引力)。新曲线由方程式(24)表达。这样，在本实施方式的模拟中，与填料模型4A和4B的间距d对应的填料间引力由方程式(24)得出。

$$ >>>T>n>>=>>27>4>>>\sigma >max>>\{>>>>u>n>>+>>d>n>>>>\sigma >n>>>[>1>->2>>>>(>>>>u>n>>+>>d>n>>>>\sigma >n>>>)>>2>>+>>>(>>>>u>s>>+>>d>s>>>>\sigma >s>>>)>>2>>>+>>>(>>>>u>n>>+>>d>n>>>>\sigma >n>>>)>>2>>+>>>(>>>>u>s>>+>>d>s>>>>\sigma >s>>>)>>2>>]>>$$

$$ >>->>>d>n>>>\delta >n>>>[>1>->2>>>>(>>>d>n>>>\delta >n>>>)>>2>>+>>>(>>>d>s>>>\delta >s>>>)>>2>>>+>>>(>>>d>n>>>\delta >n>>>)>>2>>+>>>(>>>d>s>>>\delta >s>>>)>>2>>]>\}>>$$

$$ >>>T>n>>=>>27>4>>>\gamma \sigma >max>>\{>>>>u>s>>+>>d>s>>>>\sigma >s>>>[>1>->2>>>>(>>>>u>n>>+>>d>n>>>>\sigma >n>>>)>>2>>+>>>(>>>>u>s>>+>>d>s>>>>\sigma >s>>>)>>2>>>+>>>(>>>>u>n>>+>>d>n>>>>\sigma >n>>>)>>2>>+>>>(>>>>u>s>>+>>d>s>>>>\sigma >s>>>)>>2>>]>>$$

$$ >>->>>d>s>>>\delta >s>>>[>1>->2>>>>(>>>d>n>>>\delta >n>>>)>>2>>+>>>(>>>d>s>>>\delta >s>>>)>>2>>>+>>>(>>>d>n>>>\delta >n>>>)>>2>>+>>>(>>>d>s>>>\delta >s>>>)>>2>>]>\}>·>·>·>·>·>·>>(>24>)>>>$$

图22中，显示了当填料模型4A和4B的初始间距d₀改变成各种值时，用以确定引力Tn的函数的曲线实例。特征长度δ为，例如，在初始间距d₀＝0时大约是11nm，而当初始间距d₀＝2.0nm时大约是4.8nm。因此，初始间距d₀和特征长度δ宜根据具体条件而定。

图23中，给出了初始间距d₀＝1.0nm时，填料间引力Tn的曲线实例。在填料模型4A和4B的加载形变中，表达填料间引力Tn的方程式(24)为抛物线函数，该抛物线包括：一个逐渐上升的区域YLa，在此区域内，随着间距d的增大，引力从d＝1.0nm的点P1开始逐渐增大，并到达一个顶峰σmax；一个逐渐下降的区域YLb，在此区域内，随着间距d的进一步增加，引力逐渐下降，并在预设的特征长度δ的点P2上到达零值。当间距d大于特征长度δ时，填料模型4A和4B之间不再有填料间引力作用。

在逐渐下降的区域YLb的卸载形变过程中，填料间引力Tn被设为从卸载初始值线性递减至零。例如，如果载荷从逐渐下降的区域YLb中任意一点P3开始卸去，如虚线箭头所指示的那样，随间距d的变化，填料间引力值Tn线性递减，将相关的点P3和P1连成直线。这样的卸载过程中，引力Tn由下面的方程式(25)得到。

$$ >>>T>n>>=>>27>4>>>\sigma >max>>\{>>>>u>n>>+>>d>n>>>>\delta >n>>>>(>1>->>>2>v>>max>>+sup>>v>max>2sup>>)>>->>>d>n>>>\delta >n>>>[>1>->2>>>>(>>>d>n>>>\delta >n>>>)>>2>>+>>>(>>>d>s>>>\delta >s>>>)>>2>>>+>>>(>>>d>n>>>\delta >n>>>)>>2>>+>>>(>>>d>s>>>\delta >s>>>)>>2>>]>\}>>$$

$$ >>>T>n>>=>>27>4>>>\gamma \sigma >max>>\{>>>>u>s>>+>>d>s>>>>\delta >s>>>>(>1>->>>2>v>>max>>+sup>>v>max>2sup>>)>>->>>d>s>>>\delta >s>>>[>1>->2>>>>(>>>d>n>>>\delta >n>>>)>>2>>+>>>(>>>d>s>>>\delta >s>>>)>>2>>>+>>>(>>>d>n>>>\delta >n>>>)>>2>>+>>>(>>>d>s>>>\delta >s>>>)>>2>>]>\}>>$$其

中，$$ >>v>=>>>>(>>>>u>n>>+>>d>n>>>>\delta >n>>>)>>max>>+>>>(>>>>u>s>>+>>d>s>>>>\delta >s>>>)>>max>>>>$$(u_nu_s位于卸载初始处)

……(25)

因此，当填料间引力Tn超过顶峰σmax，就在加载和卸载这一循环中形成封闭回线。由所述封闭回线的函数确定的填料间引力表明在第一个填料模型和第二个填料模型之间产生了滞后损失。如果填料间引力Tn没有超过顶峰σmax，或者是在递增区域的卸载形变中。引力也由与加载形变相同的线路确定。

填料间模型10不像基体模型3具有自身刚性。换言之，提供填料间模型10以用方程式(24)和(25)计算填料间引力。填充间模型10不在经由填料模型4A和4B的方向上被分割。因此，填充间模型10的形状仅仅取决于填料模型4A和4B各自发生形变时的相对位置，如上所述。

由方程式(24)和(25)计算填料模型4A和4B之间单位表面积的填料间引力(相互作用力)。此引力可(如图21所示)被作用于四个顶点1～4的作用力取代。因此，填料间引力作用于被夹在填料模型4A和4B之间的基体模型3上，并压缩基体模型3。

形变的计算则以上述的粘弹性材料模型2，在与上述相同的条件下进行。图24中，显示了周期性包含图18的微观结构的粘弹性材料完整模型中，实际应力和应变的关系。图24(原文图18有误)中，实线表明填料间引力被确定在填料模型4A和4B之间的粘弹性材料模型2的结果，虚线表明引力不确定的模型2的结果。与忽略填料间引力相比，本实施方式考虑了填料间引力，其模拟要困难的多。如果考虑填料间引力，表明能量损失的封闭回线的表面积要比忽略填料间引力时高出0.929％。这对应于包含了填料间引力所带来的能量损失的数量。这个数值可以通过设定参数而改变。而且，在本实施方式中，在微观结构(晶格单元)内安排了两个填料模型4A和4B，但是，通过根据，例如，被分析的橡胶材料中的填料掺混比例来改变填料模型4的密度，有可能可以得到更为精确的能量损失计算值或者是对比值。

图25中，显示了每个填料间模型10间产生的引力和间距d的关系。图中所示的参照符号r1，r2，r3，r4和r5分别对应于图19的填料间模型的元素r1，r2，r3，r4和r5。元素r1，r2和r3不能超过顶峰σmax，因为间距d的改变较小。因此，在加载形变和卸载形变过程中，引力Tn被认为是相同的。至于元素r5，由于其初始间距d原本就很大，所以不产生填料间引力。而至于元素r4，在引力Tn超过顶峰后，负载被卸去，所以被认为加载形变和卸载形变所经过的路径不同，并形成滞后回线。

下面解释将填料间相互作用量化的模拟方法。在本实施例中，也同样使用掺混了如碳黑的橡胶作为粘弹性材料，它是被分析物体。本实施方式中，势能被算为代表填料相互作用的变量。在模拟时，优选使用如图1所示的计算机设备1。

在图26中，显示了本实施方式的模拟方法的一个程序的实例。在本实施方式中，首先执行的是设定第一填料模型和第二填料模型的模型设定步骤(步骤S100)。

图27中，显示了由电子显微镜观察到的分散于橡胶聚合物中的许多碳黑6的二次粒子。如上所述，碳黑6的二次粒子由许多球形初级粒子7通过三维无规键合而组成。本实施方式的模拟集中于两个彼此间隔的二次粒子6A和6B。二次粒子6A和6B最接近的部分(即，如图27所示，二次粒子6A的一个初级粒子7a和二次粒子6B的一个初级粒子7b)被模型化，以分析如此接近的部分的相互作用。更具体的说，如图28直观所示，初级粒子7a被模型化为第一填料模型Fm1，而初级粒子7b被模型化为第二填料模型Fm2。

碳黑初级粒子7是碳原子的集合体。因此，在计算时，初级粒子7被模型化的填料模型Fm(同时指第一填料模型和第二填料模型时，简称为填料模型Fm)优选表示许多个碳原子。本实施方式中的填料模型Fm，如图29直观所示，实质上是六边形基本结构11的组合体。

图30(A)为形象化的基本结构11的透视图。基本结构11由许多个网络状排列13叠加而成。在此例中，有3～5层的网络状排列13叠加。网络状排列13包括许多个碳原子模型12(填料原子模型)和许多个将碳原子模型12在同一个平面内键合成六边形的臂状模型14。如图30(B)所示，网络状排列13是通过臂状模型14，将大约90个碳原子模型12排列并键合在规则六边形的每个顶点上而构成的。由于有臂状模型14，各碳原子模型12的相对位置不会发生变化。

填料模型Fm可以是空心的，但最好有碳原子模型12排列在其内部，这基于目前碳颗粒的常识。图31展示了设定填料模型Fm的方法的一个实例。首先，如图31(A)所示，在三维座标上给出了一个构成填料模型Fm中心的碳原子模型12。接着，将基本结构11排列在碳原子模型12的周围，以包围相应的碳原子模型。这里，基本结构11不需要形成非常规则的多面体的表面，可适当的留出空白带15。按顺序重复此步骤，就得到内部具多个碳原子模型12的填料模型Fm，如图31(C)的横截面图所示。按实际形状，填料模型Fm的最大直径宜在，如，10～200nm之间，且一个填料模型Fm中优选含有一万至十亿个碳原子模型12。

将设定的填料模型Fm置于任意的X-Y-Z座标体系，可确定包含于相应填料模型Fm中的所有碳原子模型12的重力座标中心。在填料模型Fm1和Fm2中，从1开始的固有编号被按顺序的分配给各碳原子模型12。这些编号作为相应的碳原子模型12的座标等存储于电脑设备1中。

如图28所示，第一填料模型Fm1和第二填料模型Fm2被间距r所隔开。间距r宜分布于，如，约0.1～1.5nm的范围内。在本实施例中，初始值设为0.3nm。

本实施方式中，用电脑设备，执行了计算在第一填料模型Fm1和第二填料模型Fm2间作用的填料粒子间的相互作用的计算步骤。(步骤S200)。图32中，展示了计算步骤具体过程的一个实例。

在此实例中，首先将第一填料模型Fm1和第二填料模型Fm2之间的距离r设为预设初始值(本例中为0.3nm)(步骤S201)，并将变量i，j初始化为1(步骤S202，S203)。变量i，j分别对应于分配给碳原子模型12的固有编号。这包括值i_max，j_max，它们与各碳原子模型12的数值的最大值相等。

编号为i^th碳原子颗粒模型选自第一填料模型Fm1，编号为j^th的碳原子颗粒模型选自第二填料模型Fm2(步骤S204，S205)，并计算了两个模型间的势能(步骤S206)。作用于彼此分离的原子间的势能可由多个理论公式算得。本例中，用到的是如下方程式(26)表示的Leonard-Jones势能计算公式。

$$ >>\phi >=>4>\times >0.004783>[>>>(>>0.3345>R>>)>>12>>->>>(>>0.3345>R>>)>>6>>]>·>·>·>·>·>·>>(>26>)>>>$$

Leonard-Jones势能计算公式所表达的势能中，除了范德华引力产生的引力之外，还增加了一个近距离排斥力，因此公式具有较高的精确度，并被广泛使用。然而，其它的一些计算公式也是可用的。在公式(26)中，势能φ源自碳原子模型12之间的分隔距离R的一个函数。因此，在计算势能前，应先计算第一填料模型Fm1的编号为i^th的碳原子模型12和第二填料模型Fm2的编号为j^th的碳原子模型12的之间的距离R。

通过编号为i^th的碳原子模型12的座标和编号为j^th的碳原子模型12的座标，可以很容易的计算距离R。在计算出碳原子模型的间距R后，在编号为i^th的第一填料模型Fm1中的碳原子颗粒模型和编号为j^th的第二填料模型Fm2中的碳原子颗粒模型之间产生的势能就可以通过方程式(26)计算出来。

接下来，执行将计算出来的势能值加到势能总和存储器上的操作(步骤S207)。势能总和存储器属于如工作存储器的一部分，计算出的势能的值按顺序加到它上面。通过这些值，可以得到先前单独算出的势能的总值。

计算机设备1判断目前的变量j是否为最大值j_max(步骤S208)，如果判断结果为假，则将变量j加1(步骤S209)，并重复步骤S204～207。更具体的说，编号为i^th的第一填料模型Fm1中的碳原子颗粒模型和编号为j^th的第二填料模型Fm2中的碳原子颗粒模型之间的势能是按顺序计算的，此计算值也按顺序添加到势能总和存储器上。

如果步骤S208的判断结果为真，则计算第一填料模型Fm1中的所有碳原子模型和第二填料模型Fm2中的所有碳原子模型间的势能。在这种情况下，判断变量i是否为最大值i_max(步骤S210)，如果判断结果为假，则将变量i加1(步骤S211)，变量j则初始化为1(步骤S203)，并重复步骤S204～209。

如果步骤S210的判断结果为真，则计算第一填料模型Fm1中的所有碳原子模型和第二填料模型Fm2中的所有碳原子模型间的势能，并得到势能总和。

改变第一填料模型Fm1和第二填料模型Fm2之间的间距r，对其相互作用进行求值，求得并对比各个状态下的相互作用是非常有效的。因此，本实施方式中，间距r的初始值设为0.3nm，并通过0.001nm的增量α逐渐增大到1.2nm的最大值r_max。更具体的说，是判断第一填料模型Fm1和第二填料模型Fm2之间的间距r是否为最大值r_max(步骤S212)，如果结果为假，则间距r加上增量α(本例中为0.001nm)(步骤S213)。

势能总和存储器上的值作为当前间距r下的势能总和写下并存入磁盘等设备后，清除势能总和存储器上的值(步骤S214，步骤S215)。接下来，再计算间距r增加增量α后的势能总和。另一方面，如果当前间距r被判断为是步骤S12中的最大值r_max，则过程终止。程序转至图26所示的步骤S300。

通过上述操作，得到了第一填料模型Fm1中的所有碳原子模型12和第二填料模型Fm2中的所有碳原子模型12在各种间距下的各个势能的总和。因此，可以研究出具有最小势能的填料模型间的稳定间距。

图33中，展示了计算步骤得出的结果，其中，水平轴表示第一填料模型Fm1和第二填料模型Fm2之间的间距r，垂直轴表示势能总和。图34中，水平轴表示第一填料模型Fm1和第二填料模型Fm2之间的间距r，垂直轴表示填料模型间的作用力。图34可以由图33中的曲线微分求得。如图所示，当第一填料模型Fm1和第二填料模型Fm2之间的间距r大约为0.39nm时，第一填料模型Fm1和第二填料模型Fm2之间的作用力基本上为零。

本实施方式对于明确理解填充了如填料的橡胶组合物中相应填料间的相互作用是非常有用的。通过填料和橡胶键合的控制，可以调整填料间的距离(分散性)。因此，基于填料间相互作用的计算结果，可调整填料和橡胶的分散性能，并可用于提供能量稳定的填料填充橡胶。而且，由于填料颗粒将填料原子考虑在内，可以推导出更合适的计算结果。

在不背离本发明范围或宗旨的前提下，本发明可以其它多种具体形式来体现，这对于本领域的熟练技术人员应该是显而易见的。因此，本发明并不受这里给出的细节所限制，在所附的权利要求的范围或等效性以内可作修改。