（固体力学专业论文）复合材料液体模塑成型工艺三维非等温数值模拟.pdf

中文摘要 l c m 工艺是一类树脂基复合材料成型过程的统称，主要包括树脂传递模塑 成型，结构反应注射成型，以及由两者衍，i 钳j 的t 岂如真空辅助树脂注射、树 脂膜熔渗透j _ = 艺等。该工艺技术含量较高，涉及的工艺参数较多，通过反复试 验来选定材料和确定工艺参数显然既费时又不经济。因此建立l c m 工艺成型过 程的数学模型，使用计算机对该工艺进行仿真模拟是实现低成本和优质l c m 工 艺构件的有效途径。本文的主要成果包括以下儿个方面： 针对l c m 工艺充模过程中树脂流动前沿的模拟方法，本文提出了新的模拟 算法一隐式有限元法。该算法采用填充因子追踪流动前沿，不需要使用控制体 技术。该算法计算的流动前沿分布与时间步长无关，时间步长越大计算时m 越 短。这方法特别适用于复杂模腔的充模模拟，与传统的控制体有限元法相比， 该算法能够有效地提高模拟计算速度。 本文建, 2 y 反映l c m 工艺非等温充模过程热化学和物理现象的耦合数学模 型，提出了控制方程的解耦算法。由于l c m 工艺的能量方程是典型的对流占优 的扩散方程，为避免使用传统有限元方法求解此类方程产生的数值振荡问题， 本文采用线流迎风有限元法来消除数值振荡。为了尽可能地降低截断误差，本 文还采用了时间子步方案，即将充模时间步分解为若干步长相等的温度时间子 步。 针对l c m 工艺固化过程，本文使用有限单元法进行了三维瞬态数值模拟。 编制的模拟程序可以计算得到固化过程中复合材料在任意时刻的温度分布及固 化度分布，为进一步研究在成型工艺中由于温度, r n 化度的分布引起的残余应 力奠定了基础。 作者采用f o r t r a n 高级程序语言，编制了能够模拟l c m 工艺非等温充模 过程和固化过程的三维模拟程序l c m s i m 。 关键词：l c m 工艺 数值模拟有限单元法对流扩散方程 武汉理 ：火学硕 ：学位论文 a b s t r a c t l i q u i dc o m p o s i t em o l d i n g ( l c m ) i so n eo ft h ep r o c e s s e sf o rr e s i nm a t r i x c o m p o s i t e s ，w h i c hc o n s i s t so fr e s i nt r a n s f e rm o l d i n g ，s t r u c t u r a lr e a c t i o ni n j e c t i o n m o l d i n g ，v a c u u m a s s i s t e d r e s i n i n j e c t i o n ，r e s i n f i l m i n f u s i o n ，e t c m m a y p r o c e s s i n gp a r a m e t e r sa r ei n v o l v e di nl c m t r a d i t i o n a le x p e r i m e n t a lm e t h o d sf o rt h e o p t i m i z a t i o n o fm o l da n d p r o c e s sd e s i g n c a nb et o ot i m e c o n s u m i n g a n d e c o n o m i c a l l yp r o h i b i t i v e c o m p u t e r b a s e ds i m u l a t i o no ft h e p r o c e s s c a n h e l p m a n u f a c t u r e r si m p r o v et h ep r o d u c t q u a l i t y a tl o wc o s t an e w a p p r o a c h - i m p l i c i tf i n i t e m e t h o d w a sp u tf o r w a r dt os i m u l a t et h ef l o w b e h a v i o ro ft h er e s i nd u r i n gm o l df i l l i n go fl c m t h em e t h o du s e dt h ef i l l i n gf a c t o r t od e s c r i b et h ef l o wf r o n t ，w h i c hw a sf r e eo fu s i n gc o n t r o lv o l u m e t e c h n i q u e s ， c o m p a r e dw i 恤t r a d i t i o n a lc o n t r o lv o l u m e f i n i t ee l e m e n tm e t h o d ，t h em e t h o dc a n e f f e c t i v e l y i n c r e a s et h e c o m p u t a t i o n a ls p e e d ，w h i c h w a s e s p e c i a l l y s u i t e df o r s i m u l a t i o no f c o m p l e xm o l d f i l l i n g a c o u p l e dm a t h e m a t i c a lm o d e ld e s c r i b i n gt h et h e r m a la n dp h y s i c a lp h e n o m e n a d u r i n gn o n i s o t h e r m a lm o l df i l l i n g o fl c m p r o c e s sw a se s t a b l i s h e d 1 1 1 ee n e r g y e q u a t i o no f l c mw a s t y p i c a l l yc o n v e c t i o n d o m i n a t e dd i f f u s i o ne q u a t i o n i no r d e r 幻 o v e r c o m et h e o s c i l l a t o r yn u m e r i c a l r e s u l t sc a u s e d b y t r a d i t i o n a lf i n i t ee l e m e n t m e t h o d t h es t r e a m l i n eu p w i n dp e t r o v - g a l e r k i nm e t h o dw a sa d o p t e d t 1 1 es u b t i m e s t e pm e t h o d w a su s e dt oa v o i dt h et r u n c a t i o ne l l o ra sm u c hu sp o s s i b l e af i n i t ee l e m e n tf o r m u l a t i o nf o rt h r e e d i m e n s i o n a lt r a n s i e n tc u r es i m u l a t i o no f l c mw a se s t a b l i s h e d n l e t e m p e r a t u r e a n d d e g r e e o fc u r ed i s t r i b u t i o nc a nb e o b t a i n e da ta n yt i m eb yt h es i m u l a t e dc o d e ，1 1 1 en u m e r i c a ls i m u l a t i o nm e t h o dc a nb e u s e dt o o p t i m i z et h ec u r i n gt e m p e r a t u r ec y c l e ，w h i c hc a ni m p r o v et h eq u a l i t yo f c o m p o s i t ep r o d u c t s t h em e t h o d l a i dt h eb a s ef o rs i m u l a t i o no fr e s i d u a ls t r e s s d i s t r i b u t i o n as i m u l a t i o nc o d e l c m s i m b a s e do nf o r t r a np r o g r a mw a sd e v e l o p e d ， i i 武汉理工人学硕士学位沦义 w h i c hc a l lb eu s e dt os i m u l a t et h en o n i s o t h e r m a lm o l df i l l i n ga n dc u r i n gp r o c e s so f l c m k e y w o r d s ：l i q u i dc o m p o s i t em o l d i n g n u m e f i c a js i m u l a t i o nf i n i t ee l e m e n t m e t h o d c o n v e c t i o n d i f f u s i o ne q u a t i o n n 武汉理，】：人学硕士学位论文 第1 章绪论 1 1 复合材料液体成型工艺介绍 航空航天、国防、汽车、船舶等工业领域对强度高、重量轻的构件的需 求促使了复合材料液体成型工艺( l i q u i dc o m p o s i t em o l d i n g ) 的出现，并推 动了它的发展。l c m 工艺是一类树脂基复合材料湿法成型过程的统称，它主 要包括树脂传递模塑成型( r e s i nt r a n s f e rm o l d i n g ) ，结构反应注射成型 ( s t r u c t u r a lr e a c t i o ni n j e e t i o nm o l d i n g ) ，以及由两者衍生出的工艺如真空辅助 树脂注射( v a c u u m a s s i s t e dr e s i ni n j e c t i o n ) 、树脂膜熔渗透工艺( r e s i nf i l m i n f u s i o n ) 等。 进入上世纪9 0 年代以来，随着先进复合材料的低成本化( c o s t - e f f e c t i v e c o m p o s i t e s ) 及买得起的先进复合材料( a f f o r d a b l e a d v a n c e d c o m p o s i t e s ) 等 概念的提出，l c m 工艺得到了迅猛的发展。美国国家宇航局( n a s a ) 的先 进复合材料技术( a d v a n c e d c o m p o s i t e s t e c h n o l o g y ) 研究项目就是通过工艺、 材料和设计的综合，最终实现复合材料的高强度、高减重和低成本。m c d o n n e l l d o u g l a s 航空公司在a c t 计划的支持下，采用r f ! 工艺制造商业运输机的机 翼结构j ：欧洲空中客车公司采用r t m 工艺批量a 3 2 0 发动机吊架尾部的整 流锥【3 】：在汽车工业领域，l c m 工艺更是得到了广泛的应用。生产诸如保险 杠、车盖、仪表盘、车扶手等部件h j 。 尽管l c m 工艺中的每种具体工艺有所差别。但其基本原理都是将预催化 的树脂用压力注入铺设有增强材料预成型体的闭合模腔内，固化后脱模成制 品。l c m 工艺主要包括四个方面：1 、预制件的制造和铺敷，即将增强纤维 按要求制成一定形状，然后放入模具中；2 、充模过程，即模具闭合锁紧后， 在一定条件下将树脂注入模具，树脂浸渍纤维增强体；3 、固化阶段，在模具 充满后，常温或通过加热使树脂发生交联反应；4 、脱模，当固化反应完毕后r 一一些堡堡：! 叁堂堡主堂堡堡壅 打开模具取出制品。整个工艺过程如图1 所示。 1 、预制件的制作和铺敷2 、充模阶段 3、固化阶段4、脱模 图1 1 l c m 工艺主要过程示意图 由于充模过程和纤维预成型体的设计是独立的，工装和增强材料的设计 自由度大，所以l c m2 1 2 艺是生产高性能网状结构复合材料的重要成型工艺。 与其它传统复合材料生产工艺相比，它具有许多优点，如能制造高质量、高 精度、低空隙率高纤维含量的复杂复合材料构件；模具制造容易、材料选择 范围广；成型过程中挥发少，有利于劳动保护和环境保护【5 j 。 l c m 工艺是一种技术含量较高的复合材料成型工艺，涉及的工艺参数较 多，如注射压力、模具温度、纤维预制件的渗透率、注入口和溢料口的位置 等。这些工艺参数都会对制品的性能产生直接或间接影响。通过反复试验来 选定材料和确定工艺参数显然既费时又不经济，因此，建立l c m 工艺成型过 程的数学模型，应用计算机对l c m 工艺进行仿真模拟是实现低成本和优质 l c m 工艺构件的有效途径。 1 2l c m 工艺数值模拟的研究现状 为科学地掌握和实施l c m 工艺，国内外学者积极地开展该工艺的理论和 实验研究，并且在计算机仿真分析方面取得了显著的进展，开发出了一些模 拟程序。这些数值模拟的研究工作主要是集中在l c m 工艺充模过程的仿真分 亟堡堡三盔堂堡主堂堡笙塞 析，少数提到了固化阶段的模拟方法。现在已知的充模过程的数值模拟方法， 大致可分为两类：等温充模模拟 9 - 2 5 1 和非等温充模模拟2 6 3 9 1 。所谓等温充模模 拟是指模拟树脂在预制件的浸溃过程中，忽略热化学现象；而非等温充模模 拟则在充模模拟时，考虑由树脂发生固化反应产生的热化学现象。 1 2 1 充模过程的等温模拟方法 等温充模模拟方法主要集中在两个方面：l 、如何描述树脂在纤维预制件 的流动行为；2 、如何确定树脂的流动前沿。对于前者，目前学者们达成了共 识，即采用d a r c y 定理来定量描述树脂流动行为l “。而对于后者，学者们则提 出了许多不同的解决方案。 在l c m 中，当树脂浸渍预制件后，饱和区树脂的形状及流动前沿都随时 间不停地变化，因此它属于移动界面问题( 自由边界问题) 。处理这种瞬态流 体的自由面或移动边界的算法主要有三类：欧拉法、拉格朗日法和欧拉一拉 格朗日混合法。欧拉法使用静坐标系来描述求解域形状的连续变化，因此欧 拉法的网格是固定不变的。拉格朗日法属于移动嘲格法，网格结点与流体一 起移动，网格边界与自由面吻合。欧拉一拉格朗日混合法吸取了欧拉法和拉 格朗日法的优点而避免了两者的缺点。 拉格朗日法可以对流体的运动边界进行准确地描述，但对于变形较大的 运动，为了避免网格的畸形，计算到若干步时，可能就需要重新进行网格划 分。因此对于大变形的流体计算，该方法计算成本大，应用并不广泛。对于 此类问题，欧拉法的固定网格具有很大的优势，更易于适用在复杂的模腔上。 它的主要缺点是当流体自由面不在网格边界上时，就不能准确确定界面的位 置。为了能用固定网格描述流体自由面，需要引入特定的方法，如m a c 法、 f a n 法、v o a 法【7 8 】。对于l c m 工艺充模阶段的数值模拟，大多数学者采用 的是欧拉法。 b r u s c h k e 和a d v a n i 【9 】将预制件中流动的树脂视为牛顿流体，采用d a r c y 定理对树脂的流动行为进行建模。他们采用有限元控制体方法( f i n i t ee l e m e n t ，c o n t r o lv o l u m e ) 对r t m 工艺和s r i m 工艺的充模阶段进行了数值模拟，并 且通过实验验证了模拟方法的可靠性。 一 亟塑里三! 三查兰| 兰堡主堂焦堡苎 u m 和l e e ( i 。j 采用边界元法对r t m 工艺的充模过程进行了模拟。在每一 时间步进行求解时，边界上的结点都被更新以反映流动前沿的扩展。为了验 证模拟方法，他们进行了矩形模腔在变注射压力条件下的充模实验。实验得 到的不同时刻的流动前沿及总的充模时间与数值模拟结果相吻合。 y o u n g 等 1 l - 1 2 采用有限元控制体法对s r i m 和r t m 在等温条件下的充模 过程进行了二维和三维模拟，并且编制了计算程序。他们给出了模拟的具体 方法，包括控制方程的建立、控制体的建立、边界条件、求解算法等。他们 进行了实验以验证模拟程序的准确性。 f r e d e r i c k t l 3 1 基于d a r c y 定理对r t m 工艺充模过程进行了模拟。他采用 f a n 法追踪流动前沿，这一方法类似控制体有限元法。作者指出当采用非等 参单元时，求解压力场的g a l e r k i n 有限元法和控制体有限元法完全等价，但 是当采用等参单元时，前者形成的系数矩阵是对称的，而后者形成的系数矩 阵则未必对称。 m a r k 等( 1 4 1 通过引入部分饱和度( p a r t i a ls a t u r a t i o n ) 概念，建立了一种新 的基于d a r c y 定理的有限元控制方程，实现了对r t m 等温充模过程的模拟。 他们采用隐式时间积分方案，因此可以使用固定的时间步长。 x i a o 等【1 5 】采用有限元结点控制体法( f i n i t ee l e m e n t n o d a lc o n t r 0 1 ) 对 l c m 工艺的充模阶段进行等温仿真分析。根据这个新算法，他们对现有的通 用型有限元软件进行二次开发，成功开发出了模拟程序，并通过数值算例验 证了程序的可靠性。 j o s h i a 等 1 6 1 讨论了使用传统的有限元控制体法可能会造成的质量不守恒 问题，提出了解决的方法，并将模拟的数值解与解析解作对比，证明了该方 法的有效性。 s h o j a e i a 等l j 7 】采用控制体有限元法对r t m 充模过程进行7 - - 维等温模 拟。他们采用了两种方法来处理流动前沿：传统的亚稳态方法和部分饱和度 法，并且编制了基于这两种不同方法的程序。通过一些数值算例说明基于部 分饱和度的模拟方法更适用于复杂模腔充模过程的仿真分析。 m o h a r t 等【1 8 - 1 9 1 采用d a r c y 定理描述树脂的流动，提出了追踪流动前沿的 4 一一 墓望望三盔堂堡主兰堡笙茎 显式方法和隐式方法，对r t m 工艺用于制造复杂几何形状的构件的等温充模 进行了数值模拟，并讨论了两种方法的优缺点。 国内一些学者从上世纪9 0 年代中期也开始了对l c m 工艺充模过程数值 模拟的研究 2 0 - 2 5 1 。但这些模拟工作主要集中在二维的模拟，并且研究的深度 也不够，与国外同行相比有一定的差距。 1 2 2 充模过程的非等温模拟方法 与等温充模模拟相比，非等温充模模拟要复杂得多。非等温充模模拟研 究工作的焦点除了在1 2 i 节中提到的两点外，更主要的是集中在建立反映橱 脂在充模过程中发生固化反应的热化学数学模型、控制方程的解耦算法、树 脂的流变行为等方面。 l e e 等【2 叫开发了用于模拟r t m 和s r i m 非等温充模过程的计算程序。他 们用d a r c y 定理来描述在预制件中树脂流动行为，分别建立了基于局部热平 衡和局部热非平衡理论的热传递方程，并采用控制体有限元法来处理流动前 沿。他们使用迎风方法对热传递方程进行修正，以避免由于方程中占优的对 流项所导致的数值解振荡。 y o u n g l 2 7 2 踟对l e e 的工作进行了扩展，编制了r t m 三维非等温充模模拟 程序。他使用的算法与l e e 基本一致，并通过数值算例，给出了不同工艺参 数对树脂流动行为及模具内温度场的影响。在文献 2 8 1 中，y o u n g 分析了不同 材质的模具对r t m 充模阶段和固化阶段树脂温度分布及固化过程的影响。 l i u 和a d v a n i t 2 9 1 对薄壁结构的r t m 工艺进行了非等温模拟。流动行为的 模拟是基于传统的控制体有限元法，而热传递和固化模拟则采用有限差分控 制体法，这样热传递和固化模拟就是三维的。 d e s s e n b e r g e r 和t u c k e r 3 0 - 3 1 采用局部体积均匀法建立了反映l c m 工艺物 理化学现象的质量守恒方程、动量方程( 采用d a r c y 定理) 和基于局部热平 衡的能量方程。他们通过实验讨论了热弥散现象对成型工艺的影响，并认为 在数值模拟中应该考虑热弥散。t u c k e r 在文献【3 1 】中提到了有关对流项占优的 抛物方程采用g a l e r k i n 有限元法求解所导致的振荡闯题。 武汉理工人学硕士学位论文 l o o s 等建立了r f i 非等温充模过程的三维有限元分析模型，并考虑了 预制件所受的压实力、模具的弹性变形、树脂问的热传递以及热压罐中模具 与周围环境间的热传递等对树脂在预制件中的流动和固化的影响。 m a l 等【3 3 1 对r t m 薄壳结构充模过程进行了非等温模拟。在他们的模拟中， 将弥散现象引入到了热传递方程和物质生成方程中。m a l 等通过对一些实例的 研究表明弥散现象对一般的热传递会产生影响。 l a m 等【3 4 】扩展了他们在等温模拟的研究工作，他们对通用有限元软件进 行了二次开发，实现了r t m 充模过程的非等涅模拟。文中详细地讨论了能量 方程和物质生成方程中对流项的处理，数值算例表明该方法求解稳定，精度 较高。 k a n 2 和l e e 等【3 5 - 3 7 1 采用传统的控制体对r t m 工艺充模过程进行了非等 温模拟，并进行了实验研究。在文献 3 2 】中，他们提出了一种新的捕捉流动前 沿的方法一浮动虚拟结点单元法( f l o a ti m a g i n a r y n o d e sa n de l e m e n t s ) ，该法 可以提高模拟的流动前沿的精度。 l i n 等 3 8 j 重点讨论如何解决在r t m 数值模拟中采用常规有限元法所导致的 数值振荡问题。他们采用一阶迎风来解决振荡问题，同时引入了所谓的内结点 方法( i n t e r n a ln o d em e t h o d ) ，来避免由迎风修正法产生的伪传播问题( a r t i f i c i a l d i f f u s i o n ) 。 s h o j a e i a 等【3 9 】开发了分别基于控制体有限元法和结点饱和度概念的r t m 非 等温充模模拟程序，并通过数值算例来验证算法的可靠性。 国内学者在l c m 充模阶段的非等温数值模拟上基本上是空白，作者并没有 找到相关的文献。 1 ，2 3 固化阶段的模拟方法 l c m 工艺固化阶段与其它热固性树脂基复合材料成型工艺固化阶段基本 一样。固化阶段的树脂并不发生流动，因此可以忽略因树脂流动而产生的热量 传递。固化阶段只涉及到热传导和化学反应两个方面，两者相互耦合- 许多学者针对特定成型工艺的固化阶段提出了不同的数值模拟方法。 6 武汉理上人学硕士学位论文 l o o s h 叫等人针对层合板固化过程建立了一维数学模型，并采用隐式差分法求 解。b o g e t t i 4 q 等人使用差分法对热固性厚复合材料的固化过程进行了二维数值 模拟。s u n 9 1 4 2 等人在不考虑树脂流动的情b r , t 采用非线性热传导有限元法模拟 了固化过程。c h o i l 4 3 1 等人对r t m z e 艺固化阶段进行了三维数值模拟，并通过实 验验证了模拟结果的准确性。p a r k 4 4 】等人采用有限元法对厚复合材料结构的固 化过程进行了二维数值模拟。q i 4 5 j 等人建立了热一化学一粘弹性分析模型来模 拟固化过程中复合材料的热传导、固化、残余应力等现象。 国内学者也开展了有关热固性复合材料固化阶段的数值模拟研究工作，并 取得了一定进展【4 “”。 1 3 本文研究工作的思路和内容 一般而言，只通过实验的方法来优化l c m 工艺参数，成本往往较高，而 且耗费时间，有时甚至是不可能的，因此数值模拟逐步被大家认为是一种可 行的替代方法。实际上，由于数值模拟能够很快地确定不同工艺参数产生的 效果，所以它能使得实验更加有效和经济。正是因为这些优点，对l c m 工艺 开展数值模拟研究工作具有重要意义。然而由于l c m 工艺过程是一个多场 耦合的移动边界问题，所以要实现对l c m 工艺的数值模拟，必须解决一些关 键性问题。 首先，如何准确捕捉充模过程树脂流体的流动前沿? 传统的控制体有限 元法无需对求解域进行网格重新划分因而获彳导了广泛地应用。然而这一方法 将连续充模的瞬态过程视为一系列的准稳态所组成，为了确保计算的稳定性 必须满足c o u r a n t 条件，即需将每一个准稳态的时间步长严格地限定为当前步 只有个控制体被充满，这无疑使计算时间冗长，特别是进行三维复杂模具 的充模模拟时，其计算时间可能无法忍受。因此建立准确、高效的模拟算法 尤为重要，这能使得基于此算法编制的模拟软件具有较高的工程实际应用价 值。 其次，在充模过程中，由于一些树脂会发生固化反应释放热量，因此为 武汉理工大学硕士学位论文 了进行准确模拟，不能忽略这些热化学现象。充模过程中树脂流动，热传递 和化学反应之间是相互耦合的。因此，建立非等温成型过程的耦合控制方程 及其解耦算法也是实现l c m 准确模拟的关键性问题之一。 最后，充模过程中的能量方程是典型的对流项占优的对流扩散方程，接 近双曲型，采用传统的有限元法( b u b n o v g a l e r k i n ) 进行数值求解，往往会 引起寄生振荡( s p u r i o u so s c i l l a t i o n s ) 最终导致求解失败。因此，选择合适的 求解方法，消除数值振荡、增加求解的稳定性，也就成为实现l c m 非等温数 值模拟迫切需要解决的问题之一。 针对以上所述的关键性问题，本文采用有限单元法作为基本数值模拟工 具，从以下几个方面开展研究工作。 1 通过引入填充系数概念，对质量守恒方程进行变形处理，建立追踪流动 前沿的新算法。与传统的控制体有限元法相比，这一方法能够有效的提 高模拟计算速度，特别适合于复杂三维充模模拟。 2 建立l c m 充模过程中树脂流动、热传递和化学反应过程之间相互耦合 的数学模型，并提出解耦算法。 3 采用线流迎风有限元法s u p g ( s t r e a m l i n e u p w i n d p e t r o v - g a l e r k i n ) 消 除求解能量方程时的数值振荡问题。 4 编制了模拟程序l c m s i m ，对l c m 工艺的充模过程和固化过程进行全 过程数值模拟。 武汉理l ：人学硕十学位论文 第2 章l c m 等温充模过程的隐式有限元算法 2 1 引言 l c m 工艺是非等温的成型过程，涉及到树脂流动、热化学反应等几个方 面的相互作用，然而在某些情况下，把l c m 充模过程近似看作等温过程来分 析，尽管无法获取一些参数，但可以大大简化分析过程，加快模拟速度，节 约计算资源，只要处理得当仍可以达到足够的工程近似。等温充模过程主要 是把树脂浸透预制件的过程看作是不可压缩的牛顿流体通过多孔介质的过 程，应用d a r c y 定律来描述，忽略充模过程中的热化学现象。 树脂流体在模腔中的流动是移动边界问题。在1 2 1 节中介绍了处理移动 边界问题的方法主要有欧拉法、拉格朗日法和欧拉一拉格朗目混合法。因为 欧拉法采用固定网格，能够有效节约计算成本著能更好地适应复杂几何形状 的模具，所以该方法获得了广泛地应用。为了能用欧拉法描述流体自由面， 需要引入特定的方法，如m a c 法、f a n 法、v o a 法 v - s 。大多数学者采用的 控制体有限元法( c o n t r o lv o 】u m ef i n i t ee l e m e n tm e t h o d ) 9 1 1 - 1 3 , 1 7 , 2 0 - 2 5 1 实际上 就是基于有限元的f a n 法。 传统的控制体有限元法进行等温充模模拟的基本思路是：将瞬态问题转 化为一系列拟稳态问题，即分解为若干个时涮步。在每一时间步，先计算已 填充区域的压力场，然后根据压力场计算位于流动前沿的控制体的流量，再 根据流量确定下一时间步的时间步长并更新控制体填充系数，最后根据新的 填充系数确定新的求解区域，重复上述步骤，直到树脂流体充模整个模腔， 该算法的最大缺点是为了确保计算的稳定性必须满足c o u r a n t 条件，即需将每 一个拟稳态的时间步长严格地限定为当前步只有一个控制体被充满，这无疑 使计算时间冗长，特别是进行三维复杂模具的充模模拟时，其计算时间可能 无法忍受。为了提高模拟效率，节约计算成本，本文提出了隐式有限元模拟 武汉理工大学硕士学位论文 算法。该算法采用欧拉方法描述l c m 工艺充模过程，并使用有限元隐式时间 积分方法对基本方程进行了数值求解。 2 2 控制方程的建立 2 2 1 基于欧拉法的连续方程 为了采用欧拉方法来描述树脂在模腔内的流动，引入填充系数妒，它表示 某时刻某区域内树脂的体积含量。在没有树脂的区域妒为0 ，在树脂完全充满 的区域为1 。用欧拉方法描述流体质量守恒定理，即欧拉坐标系下的连续方 程： 微分形式：粤+ v ( p ) = o ( 2 - 1 ) 积分腻肛+ p 蚺- o ( 2 - 2 ) 其中p 为流体密度，n 为边界法向矢量，u 为流体速度矢量。 为了用上述方程描述模具内树脂的流动，将填充系数代入式( 2 - 2 ) 中， 蓐们姗p 叫肌划 将式( 2 - 3 ) 展开得到： 房妒拼膨脚+ p 伽叫皿_ o 将树脂看作定常不可压缩流体时，式( 2 4 ) 可简化为 弘+ p 螂。0 对于式( 2 - 5 ) 的第二项，利用散度定理可得 l o ( 2 3 ) ( 2 4 ) ( 2 5 ) 武汉理。】：大学硕十学位论文 t + 脚= 。 式( 2 - 6 ) 是树脂流动时质量守恒的另一种表达形式。 2 2 2 树脂流动模型 ( 2 6 ) 根据l c m 充模过程的特点，为了实现数值模拟，需要进行以下简化或假 设队3 2 1 ： ( 1 ) 纤维预制件看作刚性的多孔介质，不考虑其在充模过程中的变形； ( 2 ) 将树脂看作不可压缩的牛顿流体； ( 3 ) 忽略惯性作用和毛细作用；表面张力与粘性力相比可忽略不寸； 根据上述假设，大多数学者采用达西定理6 1 来描述树脂在预制件中的流动 行为 7 - 4 8 ： h ：一竺v p ( 2 7 ) “ 其中：v 为梯度算子，“为树脂表面速度矢，p 为压力，为树脂粘度。k 为 纤维预制件的渗透率，它是二阶张量。 将式( 2 - 7 ) 代入式( 2 - 6 ) t 2 妒丢v p 阎 协8 对于没有树脂填充或填充不满的区域( 矿 1 ) ，忽略压力梯度，只考虑树 脂完全充满的区域( = 1 ) 的压力梯度，因此，式( 2 - 8 ) 可简化为 弘。丢v p 瑚 式( 2 9 ) 为积分方程，它等价于如下的微分方程： 丝：v 墨v j d a f“ ( 2 9 ) ( 2 1 0 ) 武汉理t 大学硕士学位论文 式( 2 - 1 0 ) 是本文采用有限兀法模拟l c m 工艺过程的基本控制方程。当 式( 2 - 1 0 ) 的左边项为零时，该方程退化为拉普拉斯方程。传统的有限元控制 体方法采用的基本方程实际上是这一退化方程 1 4 , l 7 。 2 2 ，3 边界条件和初始条件 为了能够求解控制方程，还必须给出边界条件和初始条件。在流动前沿 认为压力为零，在模腔壁面认为壁面的法向压力梯度为零，即没有树脂流过 模腔壁。边界条件和初始条件的数学表达式为： r 1 、馗h 奉辟竺：o ， ( 2 - 1 1 ) ( 1 ) 模腔壁= ， ( 1 1 a 一 ( 2 ) 沣射口 恒压注射 恒速注射 ( 3 ) 流动前沿： ( 4 ) 初始条件 注射口： 巩：o = l 注射1 ：3 以外区域： 乱：o = 0 2 3 有限元格式的构造 ( 2 1 2 ) ( 2 1 3 ) ( 2 1 4 ) ( 2 15 ) ( 2 1 6 ) 对式( 2 1 0 ) 的左右两边项采用有限元法进行离散处理，通过填充系数来 追踪树脂的流动前沿，从而不需要使用控制体技术a 采用标准的伽辽金法建 立有限元方程，将求解的整个区域离散为有限个单元，假定单元的求解区域 为q e ，取压力与填充系数的近似函数为 p = ，只 妒= 。 ( 2 1 6 ) q 一 = 凼pv 旧了却 一一 万 尸 p s j ! 坚堡士叁堂堡主堂堡笙壅 其中n ，为形函数，# 节点的压力值，矿，为节点的填充系数 在每个单元内，根据伽辽金法对式( 2 - 1 0 ) 进行处理，取权函数w ：n ， 并将式( 2 - 1 6 ) 代入，同时使用散度定理得 p 峨姐鼍= 卜擎n 焉n i 鳓l + 嚣 p 吣n t 耵( 2 - 1 7 ) 式( 2 - 17 ) 也可以用如下的矩阵形式表示 阻，t 矿 + 陋】p f 尸 ：= i s l 8 ( 2 】8 ) 吖r 是单元的质量矩阵，计算时将质量矩阵对角化。陋r 是单元的刚度矩阵。 沙户是单元的节点载荷列阵。其中 阻】8 = 儿m m ( 2 1 9 ) n ( 2 ，2 0 ) ( 2 2 1 ) 对于式( 2 1 8 ) 中的莎对时间的一阶导数项，采用隐式算法向后差分法处 ；= 等 ( 2 2 2 ) 将所有的单元质量矩阵、单元刚度矩阵、单元节点载荷列阵拼装成总体 矩阵，并利用式( 2 - 2 2 ) ，可以得到 阻】( 移一移) “) + f 陋i j d = f 驴) ( 2 2 3 ) 式( 2 - 2 3 ) 是通过隐式有限元方法构造的线性方程组。进一步研究式 ( 2 - 2 3 ) ，可以推导出时间步长r 的大小并不影响所计算的流动前沿和总的充 炻 r ， a 则 岍 石一肛 讲 m 陬 面一 ，盯，r = = r r 瞳 塑望堡! 查堂堡堂垡堡壅 模时间陋1 钔，本章的算例也说明了这点。因此当预成型体离散成的单元数目 很大时，可以采用较大的时间步长以节约计算成本。 2 4 求解算法 由于式( 2 - 2 3 ) 涉及到两种未知场函数，因而采用迭代法求解这个方程组。 v o l l e r 等人h 9 1 从理论上证明了对这类方程采用迭代法对任意给定的时间步长 是收敛的。求解的步骤如下： ( 1 ) 每一个时间步的开始，认为当前迭代步中的填充系数等于上一时间步的 填充系数 移，) 。”1 = 移，) “ ( 2 2 4 ) 上标n + l 表示当前计算的时间步数，下标m 表示当前时间步内的迭代次 数。 ( 2 ) 将式( 2 - 2 3 ) 进行适当变形来求解压力场，同时引入边界条件 【j r 尸 = ( s + 蛐警 ( 2 - 2 5 ) ( 3 ) 求得压力场后，利用下式计算下一迭代步中的填充系数 移乙。”1 = 移 1 + 阻】- 1 ( f j f 一a t 江 r p ) ) ( 2 2 6 ) 注意填充系数o ( 】，否则进行修正。 ( 4 ) 检查当前时间步内前后两次迭代步中的填充系数之差是否小于容 差 占，如果小于占，则进行第5 步的计算，否则令 劬) 。”1 = 彩h ”1 ( 2 2 7 ) 从第二步开始进行新的一次迭代计算。 ( 5 ) 检查如果所有节点的填充系数为1 时退出计算，否则进行下一个时间步 计算直到所有的节点填充完毕。 武汉理】二火学硕士学位论文 整个算法的流程图如下图2 1 所示 图2 - 1 迭代算法流程图 墨堡望土堂堡圭堂笪堡壅 2 5 数值算例 为了比较本文所提出的隐式有限元算法与传统显式算法( 有限元控制体 法) 的优缺点，编制出了基于这两种算法的程序。通过下面的算例进行对比。 算例中预制件的材料参数如表2 一l 所示。 表2 - 1 预制件的材料参数 预成型体的渗透率m 2 k ：1 1 0 1 0 预成型体的空隙率 9 ；0 5 树脂的粘度，p as “= 0 1 9 2 5 1 圆板型预制件 半径矗= 0 1 阮厚度日= 0 3 c 所的圆板型模具采用恒压方式注入树脂，注 射压力p = 1 1 0 5 p a ，注入半径兄：o 1 5 c m 。在这种充模情况下树脂流动前沿半 径r 与时间，关系的解析表达式为 ，：丝l 壁l 。f 旦 一墅+ 壁f ( 2 2 8 ) k 尸i 2 l 月oj 44 考虑对称性，建立四分之一圆板的有限元模型。采用八节点六面体单元 划分网格。离散的单元总数为2 5 6 ，节点总数为5 1 2 。有限元模型如图2 - 2 所 示。采用传统的显示算法计算的流动前沿如图2 - 3 所示，采用隐式算法计算的 流动前沿如图2 - 4 所示( 时间步长a t = 2 s ) 。理论充模时间为1 7 5 7 s ，采用显 式算法和隐式算法计算的充模时间分别为1 6 7 6 s 和1 6 8 s ，可以看出两者与解 析解求得的总的充模时间均吻合较好( 误差不超过5 ) 。显式计算时间为 8 9 0 s ，隐式计算时间为6 8 0 s ( 时间步长a t = 2 s ) ，从计算成本角度考虑，隐式 算法明显优于显式算法。隐式算法的另一好处是流动前沿的分布及总的充模 时间与时间步长的大小无关。分别取时间步长r 为2 s ，6 s ，1 2 s 进行计算， 计算时间分别为6 8 0 s ，4 7 0 s ，2 8 0 s 。分别绘出不同时间步长模拟出的流动半 径对应时间的关系图，并与解析解对比如图2 5 所示。从图中可以看出，不同 的时间步长计算出的流动前沿结果一样，并与解析解吻合很好。因此只关心 总充模时间时，可以选用较大的时间步长以节约计算成本。 图2 - 2 有限元网格 图2 - 3 显式算法( c v f e m ) 得到的流动前沿 图2 4 隐式算法得到的流动前沿( a t = 2 s ) 武汉理工人学硕士学位论文 量 重 里 芒 皇 量 图2 5 不同时间步长计算的流动前沿与解析解的比较 2 5 2 半球体型预制件 空心半球体模具采用恒流注入方式注入树脂。模具的半径r = o 1 5 m ，注 射半径为r 。= 0 0 1 5 m ，注射流量q = 1 l o m 3 s 。由于采用恒流方式注入树脂， 易知流动半径r 与时间，满足以下解析表达式 尺( f ) = ( 鬈十r 0 3 ) ( 2 2 9 ) 当树脂的流动行为为理想的辐射流时，注入口压力p 与时间，的关系如下 2 7 1 p = 篆【去一( 旦+ r 0 3 r2 7 r 0 广，3 】 ( 2 _ 3 0 ) 2 斌、 7。 采用八节点六面体单元划分网格。离散的单元总数为8 6 4 ，节点总数为 1 0 8 9 ，有限元模型如图2 - 6 所示。采用隐式算法，分别取时间步长r 为l j ， 2 s ，4 s ，6 s ，1 2 s 进行计算，其中a t = i s 时，模拟的流动前沿如图2 7 所示。 绘出不同时间步长计算的流动半径对应时间的关系图，并与解析解对比如图 2 - 8 所示。从图2 8 中可以看出，不同时间步长计算的流动前沿结果一致，并 与解析解吻合很好。再次表明隐式算法计算的流动前沿与时间步长的大小无 关。不同时间步长的隐式计算时间及显式方法的计算耗时如表2 2 所示。从表 2 2 中可以看出，隐式算法的计算效率较高，时间步长越大计算耗时越少。由 于采用恒流注入方式进行充模，所以注射压力随时间变化。采用不同时间步 武汉理i 人学硕十学位论文 长计算得到的注射压力与时问的关系与解析解的结果对比如图2 - 9 所示。从图 2 1 9 中可以看出，隐式算法中时间步长的大小影响压力值的精确程度，步长越 小精确度越高。这表明当需要确定模具内的压力分布时，需适当地选择较小 的时间步长= 以获取较精确的数值解。另外，图2 - 9 还表明，不同时间步长计算 出的注射压力差异随充模时间递增而逐渐减小直至为零。因此，只需确定充 模后期模具的压力分布时，依然可以采用较大的时间步长，既可满足精度要 求，也可节约计算成本。 表2 2 不同时阳j 步长隐式及显式算法的计算时问 图2 - 6 有限元网格 图2 7 隐式算法得到的流动前沿( f i l s ) 图2 - g 不同时间步长计算的流动前沿与解析解的比较 武汉理工大学硕士学位论文 图2 - 9 不同时间步长计算的注射压力与解析解的比较 2 5 3 工字梁型预制件 工字梁型预成型体的几何尺寸及有限元网格分别如图2 一1 0 ( a ) ，2 一1 0 ( b ) 所示。采用恒流注入方式注入树脂，在腹板上下两侧中央各有一个注入点， 注射流量均为o ：1 1 0 一s m 。采用八节点六面体单元划分网格。离散的单元 总数为1 2 0 0 ，节点总数为2 2 0 1 。采用隐式算法计算的流动前沿分布及最终时 刻压力场分布如图2 - 1 1 ，2 - 1 2 所示( a t = 5 s ) 。理论上总的充模时间为l o o s ， 模拟计算的结果也为1 0 0 s 。隐式算法计算的时间为3 7 1 8 s ，而传统显式算法为 1 4 7 8 0 s 。这再次表明隐式算法在计算时间上具有显著的优越性。因此隐