基于有效介质理论的物理性能计算模型的软件实现*

孙楠楠 施展† 丁琪 许伟伟 沈洋 南策文

1)(厦门大学材料学院,厦门 361005)

2)(厦门大学航空航天学院,厦门 361102)

3)(清华大学材料学院,北京 100084)

1 引 言

近年来,随着计算机技术的发展以及理论方法的完善,科研人员对一些材料的结构形成、结构与性能之间的定量关系进行了计算模拟[1−5].传统的逐点实验方法不仅耗费大量的人力物力,而且实验周期漫长,其结果也存在偶然性.预先对实验进行理论模拟计算不仅可以更好地指导实验工作的开展,也可以对实验结果进行验证,还可以缩短实验周期、节约人力物力[6,7].

复合材料具有从原子尺度到显微结构尺度再到宏观尺度的多尺度、多层次的结构特征.广义上讲所有的非均质材料都可以称为复合材料,但是不同尺度上的处理方法不同.显微结构尺度是指材料结构为纳米级以上甚至达微米级的尺寸范围,该尺度上主要研究材料显微结构的形成与演变及表征、显微结构与性能之间的定量关系以及进一步的材料显微结构设计[2].狭义上讲复合材料指的是显微结构尺度下的多相复合材料,讨论夹杂物、析出物、增强相与性能之间的定量关系.这类材料具有可调节性和可设计性,其性能也并非其组成材料性能的简单加和平均[8].因此,可以通过改变组成材料的种类、组合方式来改变复合材料的有关性能,从而设计出具有某种优异性能的新材料.

对于确定材料结构与物理性能之间定量关系的理论方法已有诸多的研究,并且形成和发展了各自的理论体系,但每种理论都有其前提条件和适用范围[2].主要的理论方法有:1)第一性原理[9,10]及均化方法[11,12]是通过先得到局部尺度(非均匀性尺度)水平上的精确解,然后利用几何周期性得到宏观尺度上的线性性能.该方法主要适应于周期结构.2)非均质材料的细观力学[13]方法是以Eshelby[14]的等效夹杂原理为基础,最先应用于解决复合材料线弹性问题,后来逐步推广到解决非线弹性问题等方面.3)有效介质理论[15]是用于确定材料显微结构与物理性能定量关系的理论方法,它假设第二相随机分布于基体相当中.有效介质理论从有效介质近似[16]起,经历了多阶段的发展过程,现已成为一种较为成熟的理论方法.近年来,国内外学者将多重散射理论[17,18]应用于复合材料物理性能计算当中,如有效质量密度[19−22]、有效弹性模量[2,23,24]、有效介电常数[2,25,26]等,推动了有效介质理论的进展.改进的有效介质理论方法[2]将有效介质理论与细观力学等其他方法结合,构造了一个系统地描述和预示复合材料显微结构与性能定量关系的理论框架.在改进的有效介质理论中,通过较少的假设,尽可能多地考虑各项显微结构因素(包括增强体长径比、体积含量、取向分布角、宏观取向角、界面性质与厚度),涵盖力学、电学、磁学等物理场的单一场性能,以及多场之间的交叉耦合性能(如压电、磁致伸缩、磁电耦合),它可以对常见的多种结构的复合材料进行计算.

目前复合材料物理性能计算理论较多,但这些理论大多模型复杂,计算步骤繁琐,缺少操作方便的模拟计算软件.本文基于改进的有效介质理论,针对复合材料弹性模量和介电常数两项物理性能,推导了理论公式,并设计开发了Composite Studio物理性能计算软件.该软件除了实现基本的内核功能之外,采用了求解和分析独立运行的工作方法.首先构建大数据量的计算组合进行求解,然后再设定不同的分析提取方式,分析显示计算结果.该工作方式可以有效地提高计算效率,降低使用者的使用门槛,利于软件的推广使用.

2 基本原理

2.1 复合材料物理性能与平均场理论

材料的物理性能包含了热、力、磁、电场下的物理性能.这些物理性能是材料前期设计时需要考虑的重要数据.

(1)式为材料物理性能的本构方程.其中,K为物理性能,J为响应场,F为内源场.性能K可以理解为单位源场F下产生的响应场J.单一场性能指的是源场F和响应场J都属于同一类型的物理场,比如弹性模量、介电常数、磁导率.

通常,在平衡(稳态)、无内源场情况下,材料对外场的响应J是一个无散量,即

复合材料的有效性能是复合材料的整体表现,通常用平均场的方式来定义,即平均场产生的平均响应的大小,如(3)式所示:

从(3)式可以看出,复合材料的有效性能求解,主要是平均场和平均响应的求解.在一些理想结构中,场的分布容易求解.对于一般的显微结构,精确地求解场分布,是一个复杂的物理问题.

采用格林函数[2,15]可以有效地对多种显微结构进行平均场、平均响应的求解,核心问题是采用格林函数求解平衡方程,即(2)式.对于非均匀介质,可以看作是在一个均匀介质的微区中引入异质颗粒,在这个局部位置上,均匀介质的均匀性遭到破坏,根据微扰观点,局部的K(x)可以表示为

(4)式中,Ko为均匀介质的相应性能参数,它不依赖于空间位置;K'(x)是一个微扰项,它包含了K(x)中的所有随机变化.进一步得到非齐次平衡方程

其中,∇2表示拉普拉斯算子,y是对应场F的势函数.该方程的解可直接由对应的格林函数求出

其中,ψo不依赖于变量K′,它是该均匀介质中的均匀势,通过分部积分可以得到局部场的解为

引入一个多重散射“t矩阵”张量

其中I为单位阵,进一步可把局部场写成

通过对局部场变量取平均,可最终得到复合材料的有效性能

(11)式是对复合材料的有效性质的普适解,但是实际上复合材料的颗粒夹杂多数是多颗粒问题,难以得到格林函数和t矩阵的精确解,因此,根据参考均匀介质的不同,通常采取两种近似来计算有效性质.当忽略颗粒间的相互作用时,可将Ko近似为基体相应的物理性能,该种近似适合于低浓度的颗粒弥散结构,对于各向同性球形颗粒的输运问题,(11)式可变为 Maxwell-Garnett[2,15,27]方程.若考虑颗粒的相互作用,可将Ko近似为复合材料的有效性能K*,称为耦合势近似,或自洽近似,此时参考介质的性质即为复合材料的性质本身,颗粒嵌入参考介质引起的微扰影响最小,对于各向同性球形颗粒的输运问题,(11)式可变为Bruggeman[2,15,28]方程.因此,Maxwell-Garnett方程与 Bruggeman方程均为有效介质理论普适解的两个特例.

2.2 复合材料几何模型与显微结构因素

复合材料主要由基体、增强体和界面组成,其真实显微结构十分复杂,为了简化,假设基体为无穷大的连续介质,而增强体用简单的旋转椭球体来模拟[2].在这种几何模型的假设下,显微结构因素包括增强体的体积分数、长径比、截止取向分布角、宏观取向角.复合材料整体基于增强体随机分布的假设,采用统计平均的方式求解等效性能,适用于大部分常规意义上的复合材料.当涉及到某些性能取决于精密的周期性结构,如光子晶体的特征散射行为,本模型并不适用.

在体积分数和长径比方面,增强体用旋转椭球体等效后,旋转椭球体的体积分数和长径比是结构的两个重要因素,这也是复合材料调整性能的两个重要手段.

在截止取向分布角方面,增强体在基体中通常是无规分布的,如图1所示.为了描述这种分布的混乱情况,引入纤维的取向分布角的概念,即每根纤维的局部坐标系的轴和宏观坐标系X3轴的夹角.而截止取向分布角θcutoff,就是所有纤维的取向分布角的最大值,这个值可以衡量取向分布的混乱程度.在计算中,假设纤维取向在截止分布角构成的分布锥内均匀分布.在这样的定义下,θcutoff为 180°表示完全无规分布,θcutoff为 0°表示完全有序,纤维全部平行排列.在实际材料中,影响取向分布的通常是工艺,例如复合材料的注塑过程的流体流动经常使增强体产生显著的择优取向.

图1 qcutoff规分布纤维的截止取向分布角Fig.1.Cut-off orientation distribution angle of randomly distributed fibers.

在宏观取向角方面,当复合材料呈现各向异性后,弹性模量和介电常数随测量方向的变化也是一个需要关注的问题.宏观取向可以由3个欧拉角来严格描述,由于复合材料通常具有∞ mm的对称性[2],为了简化软件操作,本文的计算中宏观取向角只考虑了影响最大的章动角,即两个坐标系x3坐标轴的夹角.其中x3轴为∞ mm的旋转对称轴.

3 软件设计与界面

3.1 软件设计

Composite Studio物理性能计算软件由C++语言/Qt编程完成.采用C++语言编写了计算内核,采用Qt设计了友好的人机界面.C++语言可以方便日后的升级和功能扩展,以及引入其他计算模块.Qt是一种图形化程序设计框架,方便获得可视化界面.

计算内核包含了格林函数、角度取向平均、弹性模量计算、介电常数计算、T矩阵计算以及基本的数学函数集.功能界面主要包含材料参数库、材料种类选择、功能参数选择和计算结果做图分析等.

3.2 软件界面

3.2.1 软件主界面—材料库

图2为材料参数库界面,在此界面可以添加或者删除材料,也可以对材料参数进行编辑.对于各向异性材料,勾选“anisotropic”可以输入介电常数张量、弹性常数张量的各个分量.这些材料可以作为后续计算中的基体、增强体等.

图2 材料库输入界面Fig.2.Input interface of the material library.

3.2.2 软件计算流程

图3为Composite Studio物理性能计算软件计算流程图,首先设定复合材料的显微结构因素,然后设定计算模型的类型,并求解计算,最后对结果做图分析.这些步骤都由软件的可视化界面完成,由人机交互完成计算模型信息的输入.软件预设了颗粒复合、短切纤维复合、长纤维复合三种复合材料类型.其中颗粒复合和长纤维复合分别对应了长径比为1、长径比为无穷大的特例.选择基体与增强体性质的步骤如图4(a)所示,即从材料库中直接导入基体材料和增强体材料的性质参数.图4(b)为关键的复合材料结构参数设置界面,包括了体积分数、长径比、截止分布角、宏观取向角四个结构参数.设定参数范围和间距,可以构建相应数量的组合.例如,体积分数设定10个点、长径比设定10个点、截止分布角10个点、宏观取向角10 个点,则构建的计算组合数为 10×10×10×10=104个.通过减小间距,可以构建数目非常大的计算组合数.这些组合进入求解内核,选择物理性质以及自洽/非自洽方式进行复合材料有效性质的计算.其中,自洽模型认为(4)式中的参考介质就是复合材料的性能本身(Ko=K∗),通过迭代求解.而非自洽模型通常以固定的参考介质进行求解,忽略颗粒间的相互作用.这两种模型是常用的复合材料近似方式.计算得到所有组合的有效性质后,对计算结果进行作图分析,得到有效性质随各个显微结构因素的变化关系曲线.

图3 计算流程图Fig.3.Flow chart of the calculation.

图4 显微结构参数输入界面Fig.4.Input interface of microstructure parameters.

3.2.3 计算结果分析

目前Composite Studio软件开发了复合材料有效弹性模量与有效介电常数的计算模块.图4的参数设置界面构建了一个数量很大的参数组合集合,如果把每个可变参数当作一个维度,计算结果实际上是一个多维数组.为便于分析,软件设计了对每一个参数的单独分析模式,从结果的多维数组中提取出二维的性质-参数变化曲线.

图5和图6分别为玻璃纤维/环氧树脂复合材料的有效弹性模量,以及CaCu3Ti4O12/聚氨酯复合材料的有效介电常数的计算结果.可以看出,计算结果包含了4个结构参数(包括体积分数、长径比、截止分布角、宏观取向角)下的性质-参数变化曲线,可以满足分析需要.为了便于使用,设计了参数联动同步变化的功能,即如果改变4个显微结构参数中的某一个参数,其他3个图形相应的结构参数也会进行同步地改变,省去了对参数的反复设置.最后,软件还设计了单点分析功能,可以对任意一个组合的性质进行提取.如图6所示,在“单点分析”区域,当选定确定的结构参数后,其物理性能计算结果直接呈现在下面的物理性能矩阵列表中,可供研究者直接获得计算结果.

图5 Composite Studio 物理性能计算软件计算结果—弹性模量Fig.5.Calculation results of physical performance calculation software Composite Studio—elastic modulus.

图6 Composite Studio 物理性能计算软件计算结果—介电常数Fig.6.Calculation results of physical performance calculation software Composite Studio—dielectric constant.

4 结 论

本文基于改进的有效介质理论,采用C++/Qt混合编程,开发出了一款可跨平台应用的复合材料物理性能模拟计算软件—Composite Studio,包含弹性模量和介电常数两个计算模块.软件采用了计算和分析独立工作的运行方式,创建显微结构参数叠加组合的104量级以上的结构组合数进行高通量计算,后续分析筛选不再重新计算,提高了计算效率,降低了使用者的门槛.Composite Studio物理性能计算软件可以作为一种通用的计算工具,未来将嵌入大型服务器平台,开放用于复合材料的材料设计.