（材料加工工程专业论文）注塑成型数值模拟的快速求解方法研究.pdf

郑州大学硕上学位论文摘要 摘要 注塑成型数值模拟过程中，利用有限元法对流动、保压和翘曲变形进行分析 时最终都归结为对稀疏线性方程组的求解，而且线性方程组的求解计算量在整个 数值模拟计算中占很大比重。随着单元数目的增加所要求解的稀疏线性方程组的 规模也越大，导致计算速度大大降低。尤其是在三维问题中，这个问题就更加突 出。因此，研究稀疏线性方程组的高效求解对于注塑成型数值模拟的实际应用具 有重要的意义和价值。 本文基于有限元方程组系数矩阵的稀疏性、对称性等特性，采用了链表式存 储结构，最大限度地节约了存储空间，提高了求解效率。在充分研究了大型线性 方程组的一般求解原理的基础上，采用不完全c h o l e s k y 分解预处理共轭梯度法， 保证了求解过程的高效性。由于该算法的收敛性受到矩阵条件数的影响，条件数 越大收敛速度越慢甚至不收敛。因此，本文从完全c h o l e s k y 分解出发，构造了 一种不完全分解的预处理方法，以此来降低矩阵的条件数，加快求解速度。 本文的主要工作：对有限元方程组的系数矩阵实现了链表式存储；用l d l t 直接法求得了方程组的精确解，利用g a u s s s e i d e l 迭代法求解方程组的近似解； 采用了不完全c h o l e s k y 分解预处理方法来构造预处理矩阵；使用预处理共轭梯 度法完成了对大型稀疏线性方程组的求解；通过算例比较了两种迭代法在同样的 求解精度的情况下的收敛速度，以此验证了预处理共轭梯度法与链表式存储结构 相结合求解大型稀疏线性方程组是一种高效的方法。 关键词：大型稀疏线性方程组预条件共轭梯度法链表式存储数值模拟 郑州大学硕士学位论文a b s t r a c t a b s t r a c t i nt h en u m e r i c a ls i m u l a t i o no fi n j e c t i o nm o l d i n g f e mi su s e dt oa n a l y z et h e p r o c e s so ff l o w , p a c k i n ga n dw a r p a g e ，w h i c hf i n a l l yc o m ed o w nt ot h es o l u t i o no f s p a r s el i n e a re q u a t i o n a n dt h ec a l c u l a t i o na m o u n ti sv e r yh u g ei nt h ew h o l e n u m e r i c a ls i m u l a t i o n w i t ht h ee l e m e n tn u m b e ri n c r e a s i n g ，t h es c a l eo fl i n e a r e q u a t i o nw i l lb eb i g g e r , a n df i n a l l yt h i sw i l lb r i n go nt h ef a l lo fc a l c u l a t i o ne f f i c i e n c y e s p e c i a l l yi nt h e3 - dp r o b l e m s ，t h i sw i l lb em o r ep r o m i n e n c e t h e r e f o r e ，i ti sv e r y v a l u a b l ef o rt h e o r e t i c a la n de n g i n e e r i n ga p p l i c a t i o nt os t u d yc o s t - e f f e c t i v ea l g o r i t h m a n di t so p t i m i z a t i o np r o c e d u r e st os o l v et h ee q u a t i o ng r o u p b a s e do nt h es p a r s ea n ds y m m e t r i c a lm a t r i xo ff e me q u a t i o n sa n da n a l y z i n g s o m ec u r r e n ts t o r a g e ，t h i st h e s i sc h o o s e st h ec h a i nl i s ts t r u c t u r e ，a n dw h i c hm a k e st h e l o w e s ts t o r a g er e q u i r e m e n tt oc o m p u t e ra n di m p r o v e st h es o l u t i o ne f f i c i e n c y o nt h e b a s i so ft h ea d e q u a t es t u d yo fg e n e r a lm e t h o dt os o l v el i n e a re q u a t i o ng r o u p , t h i s t h e s i se m p l o y si n c o m p l e t ec h o l e s k yc o n j u g a t eg r a d i e n t ( i c c g ) m e t h o df o rt h e e f f i c i e n c yo fs o l u t i o np r o c e s s a st h i sa l g o r i t h mc o n v e r g e n c ei sr e l a t i v ew i t ht h e c o n d i t i o nn u m b e ro f m a t r i x t h en u m b e ri sb i g g e r , t h ec o n v e r g e n c ei ss l o w e r s ot h i s t h e s i sp r o p o s e sa l le f f e c t i v ei n c o m p l e t ec h o l e s k yd e c o m p o s i t i o np r e c o n d i t i o n e d m e t h o dt of a l lt h ec o n d i t i o nn u m b e ra n di m p r o v et h es o l u t i o ne f f i c i e n c y t h em a j o rr e s e a r c hw o r ki n t h i st h e s i si n c l u d e sb u i l d i n g 叩t h el i s ts t o r a g e s t r u c t n r ef o rg l o b a ls t i f f n e s sm a n i xf r o mf e m ，u s i n gl d l tt og e tm e p r e c i s es o l u t i o n o fl i n e a re q u a t i o n ，u s i n gg a u s s s e i d e lt og e tt h ea p p r o x i m a t es o l u t i o n , c o n s t r u c t i n g i n c o m p l e t ec h o l e s k yp r e c o n d i t i o n i n gm a t r i x 。u s i n gp c gt of i n i s ht h ef a s ts o l u t i o no f t h el a r g e - s c a l es p a r s el i n e a re q u a t i o ng r o u p n u m e r i c a le x a m p l e ss h o wt h a tt h e c o m b i n a t i o no f p c gm e t h o da n dl i s ts t o r a g es t r u c t u r ei sh i g he f f e c t i v e k e yw o r d s ：l a r g e - s c a l es p a r s el i n e a re q u a t i o ng r o u p ，c h a i nl i s ts u u c t u r e ，n u m e r i c a l s i m u l a t i o n ，p c g i i 郑州丈学硕十学位论文绪论 第一章绪论 1 1注塑成型的原理及研究意义 塑料工业是国民经济中的一个非常重要的行业。由于塑料的机械性能和加工 性能优良，而且具有质量轻、耐腐蚀。电绝缘性能好、强度高等优点，引起了人 们的关注，获得了迅速的发展，在汽车、家电、仪器仪表、建筑装饰等领域得到 了广泛的应用。 注射成型是根据金属压铸成型原理发展而来的塑料制品的主要加工方法，使 用注塑机和注射模具把塑料原材料转变成塑料制品。基本原理是利用塑料的可挤 压性与可塑性，首先将松散的粒状或粉状物料从注射机的料斗送入高温的机筒内 加热熔融塑化，使之成为粘流态熔体，然后在柱塞或螺杆的推动下，以较大的流 速通过机筒前端的喷嘴注射进入温度较低的闭合模具中，经过保压、冷却后，开 启模具便可从模腔中脱出具有一定形状和尺寸的塑料制品。 注塑成型在整个塑料制品生产行业占有非常重要的地位。早期的注塑成型方 法主要用于生产热塑性塑料制品。随着塑料工业的迅速发展以及塑料制品应用范 围的不断扩大，注射成型方法己经推广应用到热固性塑料制品和一些塑料复合材 料制品的生产中。目前，除了少数几种塑料外，几乎所有的塑料都可以采用注塑 成型。据统计，注塑制品约占整个塑料制品总产量的3 0 ，全世界每年生产的注 射模数量约占所有塑料成型模具数量的5 0 。 模具在注塑成型过程中处于核心地位，作为聚合物成型的重要工艺装备， 其设计与制造水平直接关系到产品的质量、品种及更新速度。模具工业是国民经 济的基础工业，塑料模具是塑料工业发展的瓶颈之一。在橡塑制品加工中，材料 本身的特性以及复杂的加工条件使材料成型过程经历了相当复杂的变形历史，如 固体输送、熔融、熔体输送、流动、压实、相变、结晶、分子取向、纤维取向、 翘曲变形等，复杂的制品形状、边界条件及材料参数的不确定性使得问题更加复 杂化。由于问题的复杂性，长期以来，高聚物成型加工过程的控制和模具设计与 制造主要依赖于工艺人员和设计人员的经验和技巧，设计的合理性只能通过试模 才知道，制造的缺陷主要靠修模来纠正，致使模具及高聚物产品的设计与制造周 期长、成本高、档次低。随着新材料和新成型方法的不断出现，问题更加突出， 郑州大学硕上学位论文绪论 致使我国模具技术水平大大落后于国外。现在模具生产只能满足需要的6 0 左 右，每年迸日精密复杂模具高达数亿美元，这还未包括随设备和生产线带进来的 模具。随着汽车、电子、家电、通讯等的进步和其产品更新换代不断加快，对模 具产品的需求量正日益增大。模具生产技术水平的高低，已成为制约我国制造业 现代化发展的关键问题。特别是随着国防和国民经济的飞速发展，单靠跟踪与模 仿永远不能实现由制造大国向制造强国的转变，应该探索和创新成型和加工的工 艺过程，加强模具优化设计及其反问题研究对于提高模具的设计水平至关重要。 随着计算机技术的发展和人们对计算力学、流体力学、聚合物加工流变学。 传热学等学科研究的深入，成型模拟( c a e ) 技术与注塑成型技术相结合，为改 变这种状况提供了新的手段。c a e 技术通过建立高聚物成型过程的物理和数学模 型，构造有效的数值计算方法，借助于计算机仿真模拟确定加工条件的变化规律。 预测制品的结构和性能、确定高聚物制品和模具设计参数及工艺条件的最佳方 案，使高聚物成型加工及模具设计建立在科学分析的基础上，为优化模具设计和 控制产品成型过程以获得理想的最终“定构”提供科学依据和设计分析手段，指导 高分子成型，提高高分子材料使用水平。成型模拟技术使高聚物成型加工过程在 流场、力场、热场等作用下所出现的各种物理现象和化学变化的描述更加数学化 和定量化，从而使加工成型从一项实用技术变为- - n 应用科学。塑料成型c a e 技术对缩短产品的开发周期、提高塑料制品的质量、降低生产成本具有很重要的 意义【”。 1 2 注塑成型模拟的研究现状1 2 1 和发展趋势 由于塑料注射成型加工的重要性和塑料成型模拟、优化与控制技术对实际 生产过程的重要指导价值，许多国家的科研机构和高等院校都在这方面进行了研 究。 美国早在2 0 世纪5 0 年代就开始了对塑料成型过程的研究。m c k e l v e y ，g o r e 和s q u i r e s 等人在这方面都做了出色的工作。b e m h a r d 【3 】对2 0 世纪5 0 年代聚合物加 工模拟方面的成果进行了详细的总结。后来，m c k e l v e y 【4 1 基于质量守恒、动量 守恒和能量守恒对聚合物加工进行系统分析，t a d m o r 和k l e i n l 5 1 对注射机中的塑 化过程进行了系统的研究，并于1 9 7 0 年首次提出了一个包括固体输运、熔融、熔 体输送的完整的模型。k l e i n 和m a r s h e l l 于1 9 6 8 年写了第一本关于聚合物加工计算 2 郑州大学硕上学位论文绪论 机模拟的专著 g i l m o r e 和s p e n c e r 在2 0 世纪5 0 年代初期对注塑成型进行了研究，这可能是对 注塑成型进行的最早研究。此后，人们开始对注塑成型的许多方面进行研究。 t o o r 、b a l l m a n 及c o p p e r 于1 9 6 0 年采用数值方法对塑料熔体的充模过程进行了分 析。从2 0 世纪6 0 年代中期到7 0 年代中期，s t e v e n s o n 等人对塑料熔体在矩形和中心 浇口园盘薄壁型腔中的流动行为进行了详尽的研究，建立了熔体一维流动的数学 模型，为系统的开展注塑模c a d 的研究打下了坚实的基础。在2 0 世纪7 0 年代与8 0 年代初，不少研究者对熔体在任意形状薄壁型腔中的流动进行了研究。k a m a l 和 k e n i g 对半圆形模腔的注塑过程进行了模拟，他们的工作包括问题的提出、模型 的选择、方程的求解和结果的验证。w i l l i a m s 和l o r d 对塑料熔体的圆管一维流动 进行了详尽的研究，数值分析与实验结果吻合较好。对熔体充填过程进行二维分 析的研究始于2 0 世纪7 0 年代中期。b r o y e r 等人采用流动分析网络法( f l o wa n a l y s i s n e t w o r k ) 对塑料熔体在模腔中的二维等温流动过程进行了分析。k r u e g c r 对熔体 在含有嵌件的变厚度矩形型腔中的流动行为进行了分析。a u s t i n 首次对真实塑件 的成型过程进行了研究，并采用l a y f l a t 技术去处理复杂的几何形状。由于注塑件 的结构复杂，熔体在模具型腔中的流动是非等温非稳态过程，而且熔体表现出很 强的非牛顿特性，无法对该过程进行真实分析，必须进行适当的简化。h i e b e r 和 s h e n 作出了重要贡献，他们将h e l e s h a w 模型推广到非牛顿流体的非等温流动过 程而提出了广义h e l e s h a w 模型，并采用有限元、有限差分耦合法求解控制方程。 后来，w a n g 和h i e b e r 沿用h i e b e r 和s h c n 提出的数学模型并采用流动分析网络的基 本思想，采用控制体积法确定熔体流动前沿位置，对三维薄壁件的成型过程进行 了模拟。目前，h i e b e r 和s h e n 提出的有限元有限差分控制体积方法被推广应 用于气体辅助注射成型及纤维增强材料的注射成型中，在大多数情况下部可以较 好预测型腔内速度场、压力场、温度场、气穴、熔接线以及应力集中区域等。 利用广义h e l e s h a w 模型对充填过程进行分析可以得到近似的前沿位置及充 填时间、压力分布等。但该模型忽略了熔体在厚度方向上的流动，也没有考虑熔 体的弹性，因此存在以下问题：( 1 ) 对于厚壁件，流体在厚度方向上的流动不 能忽略，若仍按该模型计算，则结果存在较大的误差。但在对塑件的质量进行分 析时，常常需要对分子取向及由于流动引起的残余应力进行准确预测，这就需要 对流动| ； 沿附近的速度场进行分析，也就是说对喷泉效应进行分析，h e l e s h a w 郑州大学硕士学位论文绪论 模型无法解决这类问题。( 3 ) h e l e s h a w 模型无法对固壁附近、两股或多股流体 汇合处以及在浇口、筋、厚度突变的地方( 在这些地方流体的剪切应变和拉伸应 变对应力场的影响很大) 的流体流动行为进行很好的描述。 ( 4 ) 对于某些聚合 物，弹性行为对塑件的质量有较大的影响，对这种行为处理不好，往往会使塑件 产生表面缺陷，如，表面斑痕，光洁度不均匀等。除了在薄壁型腔中高速流动的 聚合物呈现出较弱的弹性行为外，在其他情况下，聚合物均表现出相当的粘弹行 为，h e l e s h a w s l 型忽略弹性的简化在一些情况下己显得不太合理。基于这些原 因，人们研究新一代的模型对塑料成型过程进行更合理的描述，从而将塑件的微 观结构和物理性能与流场和温度场结合起来。新一代的模型应该具有以下三个特 点：( 1 ) 考虑熔体的粘弹行为。( 2 ) 考虑材料的物性对温度、压力的依赖性。 ( 3 ) 对流场采用更精确的描述，如采用三维模拟，真实边界条件等。尽管对注 射成型过程的理解日益深入，但对型腔充填时喷泉流的处理及利用粘弹性模型对 成型过程进行分析( 包括收缩、取向) 的研究仍有待完善。 目前，注塑成型充模过程的数值模拟技术已经有了长足的发展，但该领域仍 面临许多亟待解诀的重要课题，主要表现在1 6 1 ： 1 对高聚物成型过程的认识还有待进一步加深。例如对于注塑制品的翘曲 和收缩机理还缺乏深入的物理解释；高聚物的粘弹性、结晶、分于取向以及各向 异性对收缩和翘曲的影响还没有一个合理的数学描述；温度、压力对熔体交汇处 制品微观结构及力学性能的影响还缺乏深刻的认识。 2 考虑温度、压力和冷却速率等因素对材料行为的影响，采用更精确的模 型对材料的热物理性能、流变性能等进行表征，提高数值分析结果精度。 3 对聚合物在复杂区域中，如浇口、筋和厚度突变等处的真三维粘弹流动 传热数值计算理论与方法的研究仍不成熟，在应变率很大的情况下，对粘弹性问 题的求解仍存在困难。对于小尺度或复杂形状制品的成型过程模拟，必须建立完 全的三维粘弹流动计算模型，同时构造高效和健壮的求解器和求解过程。 4 现在的流动分析大多是在已有模具设计或设计方案情况下，进行流动分 析的。然而借助流动模拟虽可对设计进行评价，提出修改意见，但因流动分析所 需时间长，反复运行代价昂贵，故大多只分析一次或几次，是很难找到最优设计 的。因此，如何改进流动模拟算法，提高计算速度，缩短分析时间，为模具优化 设计创造条件是有待进一步研究的问题。 4 郑州大学硕十学位论文绪论 1 3 数值模拟与大型方程组的求解 在聚合物加工中，由于制品形状复杂、材料本身变形的特性以及复杂的边界 条件和不确定因素使材料成型过程经历了相当复杂的变形历史和相态变化，如流 动、压实、固化、相交、分子取向、纤维取向、翘曲变形等，成型工艺条件对成 型过程中材料所经历的受力状态和热状态有着重要的影响，直接关系到制品的力 学性能和成型效率等。面对复杂的变形历史和相态变化，在对这些问题进行分析 时，传统的解析解法受到很大的制约。以有限元法为代表的数值模拟方法为改变 这种状况开辟了崭新的应用前景。 数值模拟方法的一个重要问题是线性方程组的求解问题。任何一种数值方 法，最终都归结为一个线性方程组的求解。从理论上讲，无论是有限元，还是有 限差分等数值方法对于任何频率、任何形状和大小的物休，只要建立了算子方程 模型，都能求出指定精度的解以满足工程的需要。但实际上，求解过程却受到计 算机内存、计算速度和数据长度等条件的限制，特别是计算速度的影响。对于小 规模的问题，由于所得到的线性方程组的阶数较小，无论使用任何一种方法求解， 都能快速而准确的求出方程的解【7 1 。但是现在由于注塑成型过程的复杂性，最后 由有限元、有限差分离散所得到的往往是一个超大规模的线性稀疏方程组。由于 该线性方程组的规模很大，用一些经典的解法，如高斯消元法，直接分解法等， 往往需要很大的计算内存以及计算机求解运算时间。甚至根本不能求出最终的 解。特别是当用有限元法离散网格划分的很细时，该线性方程组的阶数就会急剧 的上升，从而导致求解变得很困难，甚至由于不能求解或求解很慢而使得有限元 技术的应用范围受到很大的限制 引。因此如何解决计算精度与计算速度的矛盾一 直是人们特别关注的问题，吸引了众多的学者去研究和开发。采用各种加速手段 力争减少求解问题所需的计算量和存储量，以提高计算精度和计算效率。而共轭 梯度方法及其与其他迭代方法及各种加速手段的联合运用等，正是适应了这种要 求且在实际应用中都大大提高了计算效率，解决了很多实际问题，因此在快速求 解方面受到了众多研究者的青睐。 超大规模的线性方程组，由于直接求解很困难，一般都采用迭代方法进行求 解。迭代方法是求解注塑成型过程中复杂问题的最有效的方法之一，其具有通用 性好、存贮量小和精度易于控制等优点。再配合适当的预处理技术，迭代解法能 5 郑州大学硕士学位论文绪论 够在较少的步骤内收敛到问题的解，并达到精度要求，由于该方法充分利用了矩 阵的稀疏性、对称性等特点，所以迭代解法无论从时间还是空间上都是最优的， 在近些年得到了很高的重视，并广泛的应用于工程技术领域。但是，当所研究的 精度要求较高时，就要对所计算的问题进行更为精细的离散，从而导致所求未知 数个数的急剧增加，经典迭代方法的收敛速度将随着未知数个数和迭代次数的增 加而减小，这势必带来所用机时的猛增，给我们的数值模拟计算带来很大的约束， 从而使得其在工程应用中受到限制。因此，如何选择适当的迭代解法和相应的预 处理技术就成为高效求解有限元大型稀疏线性方程组的关键。目前对大型方程组 求解的研究主要有两方面的内容： 一是算法的加速求解问题 9 1 。主要有并行运算技术，预处理技术等。并行算 法就是对于同一个线性方程组，用多台计算机来同时并行求解。并行算法需要占 用较多的计算机资源。而预处理技术就是在线性方程组的迭代计算过程中，将所 求的线性方程组的系数矩阵通过一定的变换，来改变矩阵的某些特性，从而显著 提高算法的迭代求解速度。尤其是在对大型方程组进行快速求解研究时，共轭梯 度法的预处理技术的研究占有很重要的一方面。 另一种就是研究快速有效的迭代算澍1 0 1 。由于超大规模的线性方程组的求 解，一般都是运用迭代法求解，而不同的迭代法的收敛速度会有明显的不同，甚 至有可能会发散，因此研究一种迭代方法的性能以及研究一种新的快速算法变得 很重要。但是这与在现有的算法的基础上对算法进行适当的优化相比，实不是件 容易的事。 代数方程组的数值方法及其理论，是数值分析与科学计算的基础与核心。高 性能的线性与非线性代数求解器，是获得高计算效率和准确计算结果的需要和保 证。因此，随着计算机技术的迅速发展，并行算法、共轭梯度迭代算法及其预处 理技术，己成为当前数值代数的国际前沿研究方向。 1 4 本文的研究背景 注塑成型数值模拟过程中，在对流动过程、保压过程和制品的翘曲程度进行 分析时均用到有限元法。对所要研究问题的一般方程采用有限元法的基础理 论里兹变分法和加权残余法进行离散，得出每个单元的方程，然后组装成大 6 郑州大学硕士学位论文绪论 型方程组。对工程问题求数值解，要把整个连续体离散成单元的形式，然后对每 个单元进行求解，得出单元刚度矩阵，再和施加在该单元上的附加条件就形成了 单元的方程组，最后把单元刚度矩阵加入总体刚度矩阵中，也同样得到所要求解 问题的大型方程组。 理论上，随着单元网格划分得越密，单元的数目越多，最终所求得的解就越 精确。但是随着单元数目的增加，所要求解的线性方程组的规模也就越大，所形 成的系数矩阵的阶数就越高，从而在求解时就会造成计算速度的大大降低。尤其 是在求解三维问题时，方程组的规模更大，求解起来就更加的困难，所需要的求 解运算时间就更长，给数值模拟计算带来诸多不便。 目前在流动和保压分析中，主要采用g a u s s - - s e i d e l 迭代算法；而在翘曲分 析中一般采用直接法求解。但是直接法受到问题规模( 即内存量) 的限制，即在 存储数据时需保存系数矩阵中夹杂于非零元素之间( 即高度轮廓线以下) 的零元 素，因为在计算过程中它们要成为非零元素。因此，不仅增加了对计算机存储的 要求，而且影响了计算效率。在大型、超大型问题中这种零元素在系数矩阵中所 占比例很大，常常超过9 0 ，因此直接解法用于此类问题时求解效率是不高的。 另外，在计算过程中直接法不能对解的误差进行检查和控制。解的误差主要是由 计算机的有效字长，方程组的阶数，特别是系数矩阵的性态决定的。在一定的字 长条件下，方程组的阶数越高，计算过程中累积的误差越大。另一方面，系数矩 阵的性态越差，一定有效字长造成的截断误差对解的影响也很大。而系数矩阵的 性态主要体现在矩阵的条件数( 即系数矩阵的最大特征值和最小特征值之比) 。 计算规模越大，以及网格内不同单元间尺寸的差别、材料性质的差异、各单元自 身比越大，则系数矩阵的条件数越高，亦即矩阵性态越差，从而导致解的误差越 大，甚至导致求解失败。因此，这种方法一般不用在求解大规模的问题上。对于 大型、超大型的问题，采用迭代法进行求解可以避免上述的两个缺点，即不用存 储夹杂在非零元素之间的零元素和可以对解的误差进行检查和控制。但是，在注 塑成型数值分析中，g a u s s - - s e i d e l 迭代法在求解过程中会出现迭代收敛很慢甚 至发散的情况，导致计算效率低下。因此，为了快速求解大型线性方程组，就要 寻求一种收敛较快的迭代算法，而这也正是本文所要研究的内容。 7 郑州大学硕士学位论文绪论 1 5 本文主要工作 本文针对有限元大型稀疏线性方程组的特性，为了达到快速求解的目的，对 稀疏矩阵进行了适当的预处理，然后采用了不完全c h o l e s k y 分解预处理共轭梯 度法高效求解了线性稀疏方程组。基于有限元总刚矩阵的大规模、稀疏性、对称 性等特性，采用了链表式存储结构，不仅节省了存储空间，而且便于存取运算。 而且把这种链表式存储结构与预处理共轭梯度法相结合，经过算例的验证，大大 提高了有限元大规模稀疏线性化方程组的求解效率。 本文的主要工作概括如下： 1 、利用总刚矩阵稀疏性、对称性等特性，采用链表式存储策略，仅存储总 刚矩阵下三角部分非零元素，不仅存贮规模小，而且便于存取运算，有 效节省了计算时间。 2 、采用l d l l 直接法对线性方程组求精确解，为与使用迭代算法求出的解 作比照；在g a u s s s e i d e l 迭代法中引入了松弛因子，此松弛因子的引入 对加速方程组的收敛起了非常重要的作用，在与l d l o 直接法同样精度 的情况下求方程组的解。 3 、由共轭梯度法的基本原理出发，分析了影响该算法的迭代收敛速度的因 素，其中最主要的影响因素是矩阵条件数。条件数大的迭代收敛较慢甚 至发散，条件数小的则收敛的较快。 4 、对常用的预处理方法进行了总结，着重研究了不完全c h o l e s k y 直接分解 预处理方法。 5 、用v i s u a lc + + 编写了求解器，通过具体的算例验证了在同样的精度情况 下，预处理共轭梯度法比g a u s s s e i d e l 迭代法更加高效。 郑州丈学硕学位论文总刚矩阵的特点和存储方式 第二章总刚矩阵的特点和存储方式 2 1 有限元法概述 有限元法1 1 1 1 又称有限单元法，其作为一种强有力的工程分析方法被广泛应用 于各种研究领域( 机械、电子、流体力学、热分析等) ，是求解复杂的工程问题 的一种近似的数值分析方法。 t t 有限元，这个名词是c l o u 曲于1 9 6 0 年第一次使用的【1 2 1 ，从此以后，有限 元应用的著作逐渐增多，现在有许多杂志都主要地致力于有限元法理论的发展和 应用。5 0 年代中期至6 0 年代末，有限元法出现并迅猛发展，由于当时理论尚处 于初级阶段，计算机的硬件及软件也无法满足需求，有限元法和有限元程序无法 在工程上得到普及；到6 0 年代末7 0 年代初出现了大型通用有限元程序，它们以 功能强、用户使用方便、计算结果可靠和效率高而逐渐形成新的技术商品，成为 结构工程强有力的分析工具。 电子计算机的发展对有限元法的发展和应用有着决定性的影响。有限元法要 求解大规模的联立代数方程，未知数多达几万甚至几十万个，没有高速度大容量 的计算机是不可能实现的。电子计算机的发展和应用对有限元法的发展和推广应 用创造了良好的物质基础。 有限元法是基于里兹变分方法 1 3 1 与伽辽金加权残余法的一种数值模拟方法， 是在变分原理【1 4 1 的基础上结合有限差分方法中的区域离散思想而发展起来的。其 基本思想是将空间上的连续的求解区域离散为一组有限个、且按一定方式相互联 结在一起的单元的组合体。由于单元能按不同的联结方式进行组合，且单无本身 又可以有不同的形状，因此可以模型化几何形状复杂的求解域。有限单元法作为 数值分析的另一个重要特点是利用在每一个单元内假设的近似函数来分片地表 示全求解域上待求的未知场函数。单元内的近似函数通常由未知场函数或及其导 数在单元的各个结点的数值和其插值函数来表达。这样一来，一个问题的有限元 分析中，未知场函数及其导数在各个结点上的数值就成为新的未知量( 即自由 度) ，从而使一个连续的无限自由度问题变成离散的有限自由度问题。一经求解 出这些未知量，就可以通过插值函数计算出各个单元内场函数的近似值，从而得 到整个求解域上的近似解。 9 郑州大学硕+ 学位论文总刚矩阵的特点和存储方式 在工程实践中，有限元分析软件与c a d 系统的集成应用使设计水平发生了 质的飞跃，主要表现在以下几个方面： 1 ) 增加设计功能，减少设计成本； 2 ) 缩短设计和分析的循环周期： 3 ) 增加产品和工程的可靠性； 4 ) 采用优化设计，降低材料的消耗或成本： 5 ) 在产品制造或工程施工前预先发现潜在的问题； 6 ) 模拟各种试验方案，减少试验时间和经费； 7 ) 进行机械事故分析，查找事故原因。 实际上，有限元法是对古典近似计算方法的归纳和总结。它吸收了有限差分 方法中的离散处理的内核，又继承了变分计算中选择插值函数( 试探函数) 并对 区域进行积分的合理方法。插值函数是定义在单元内而不是整个区域，这样就克 服了古典变分法由于不作离散而不能求解复杂问题的缺点，在有限元方法中由于 积分运算在离散的单元中进行，就充分考虑了不同单元对节点的贡献，从而克服 了有限差分法中不考虑单元本身特性的缺点。 有限元法的最大优点是：对不规则几何区域的适应性好。其网格剖分十分灵 活，可以根据实际物体的复杂程度和各部分的重要系数调整网格的疏密和大小， 从而使得单元能够逼近物体的实际几何形状。同时随着单元数目的增加，单元尺 寸的缩小，或者单元自由度的增加及插值函数精度的提高，有限元解的近似程度 将不断增加。如果单元满足收敛要求，近似解将最终收敛于真实解。 有限元法求解的基本步骤： l 、结构的离散化。有限元的第一步，也是最重要的一步，所选择的离散方 式将影响到计算机内存需求、计算时间和数值结果的精确度。因此在把结构或连 续体分割成许多单元进行分析时，必须用适当的有限元素把结构模型化，并确定 单元的数量、类型、大小和布置。 2 、从区域或结构中取出其中一个单元来研究。选择适当的插值模式或位移 模式近似地描述单元的位移场。由于在任意给定的载荷作用下，一般情况，取插 值模式为线性插值模式，对于复杂结构的位移解不可能预先准确地知道，此时， 可以把插值模式取为多项式形式。 1 0 郑州大学硕上学位论文总刚矩阵的特点和存储方式 3 、建立方程组公式。导出方程组的公式的方法一般用里兹变分法和伽辽金 方法。这里我们首先考虑伽辽金方法。 在建立方程组公式过程中，需要三个步骤。首先，应用两种方法中的任一种， 写出单元方程式；其次，将单元方程对所有单元求和，得到方程组，这个过程叫 做组合；最后再应用边界条件来得到方程组的最终形式。在用计算机实现此计算 的过程中，这三个子步骤通常不是分立的，相反，而是相互交织在一起的。单元 矩阵的生成和边界条件的强加通常发生在组合过程中。 4 、方程组的求解 方程组的求解是有限元分析的最后一步，最终的方程组通常是下列形式： k 8 = p( 2 1 ) 该方程是确定型的，它通常是从非齐次微分方程或非齐次边界条件或从它们 两者兼有的问题中导出的。在注塑成型过程中，确定的方程组通常是所施加的载 荷和边界条件有关的。由有限元法导出的线性方程组通常为大型的线性方程组， 且为高度的稀疏矩阵，求解起来相当的困难。通常该步骤占总步数的9 5 左右。 因此，如何快速有效求解有限元方程组关系系到有限元法的实用性，具有非常重 要的意义【1 5 1 。 2 2 网格划分 利用有限元法对连续体进行离散，剖分网格【1 6 】时要遵循一定的原则或要求。 有限元网格应满足以下的要求： 1 ) 、单元之间不能相互重叠，要与原物体的占有空间相容，即单元既不能落 在原区域之外，也不能使原区域边界内出现空洞； 2 ) 、单元应精确逼近原物体。所有原区域的顶点都应取成单元的顶点；所有 网格的表面顶点都应落在原区域表面上；所有原域的边和面都被单元的边和面所 逼近； 3 ) 、单元的形状合理。每个单元应尽量趋近于正多边形或正多面体，不能出 现面积很小的二维尖角元或体积很小的三维薄元； 4 ) 、网格的密度分布合理。分析值变化梯度大的区域需要细化网格； 5 ) 、相邻单元的边界相容，不能从一个元素的边或面的内部产生另一元素的 郑州大学硕士学位论文总刚矩阵的特点和存储方式 顶点， 6 ) 、划分单元的个数，视计算要求的精度和计算机的容量而定。单元分得越 多，块越小精度越高，但需要的计算机容量越大。因此，需根据实际情况而定； 7 ) 、尽量利用对称性，以减少计算量。 2 3 总刚矩阵的合成 从单元刚度矩阵的意义可以看出；单刚中处于第i 行、第_ ，列的元素的值的 大小体现了第f 号结点对于第，号结点的影响或者说是“贡献”。所以单元刚度矩 阵在合成总刚度矩阵( 简称总刚) 的时候也应该体现该影响，由此我们可以得到 单刚合成总刚的方法。具体上讲就是：先求出各个单元的刚度矩阵k ，然后将 其各个元素按照结点的编号置于总刚度矩阵相应的位置上，对于在某个单刚中没 有出现的结点编号的位置用0 补上，然后将总刚矩阵相同位置上的元素进行相 加，就可以形成总刚度矩阵足，如图2 1 示， 1 ( 2 ( 1 0 ) 图2 1 单元的局部编码与整体编码 f i 9 2 1l o c a lc o d ea n dg l o b a lc o d eo f e l e m e n t 单元刚度矩阵是与单元的节点顺序相关的，为了使其在总刚中的位置近可能 集中在对角线附近，一般采用逆时针的顺序。上图中，1 、2 、3 代表其局部结点 号，( 3 ) 、( 7 ) 、( 1 0 ) 代表单元所对应的整体结点号。由整体结点号将所有单刚 相加，就构成了总刚： 足= k 8 ( 2 2 ) 1 2 郑州大学硕上学位论文总刚矩阵的特点和存储方式 其中， k l l k e ：i ， l 墨。 3 行 7 行 1 0 行 很显然，如果单元划的越密，那么节点数目就越大，总刚度矩阵的阶数也 将越大，所需要的存储空间就越大，计算所需要的时间就越长。 2 4 总刚矩阵的特点 由于总刚中非零元素体现的是不同节点之间的相互作用关系，而每个节点 只与该节点相邻的节点和单元之间存在相互影响关系，所以总纲中只有靠近对角 线的少数元素为非零元素，而其他位置上的元素为零。因此在求解有限元方程组 时就要考虑到刚度矩阵的存储方法。而且有限元的求解效率很大程度上不但取决 于线性方程组的解法的选取而且刚度矩阵的存储方案也同样影响着求解速度。所 以，为了提高求解效率，使求解过程达到最优，我们有必要研究一下总体刚阵的 一些特性： l 、大规模，在实际建模中，最后产生的矩阵一般为大规模的线性矩阵。当 然矩阵的大小及规模还与我们所要研究的注塑成型数值模拟过程中的物理模型 有关。所要模拟的区域越大，所划分的网格越细，刚度矩阵的规模就会越大。随 着有限元网格的划分越细，物理模型的尺寸越大，线性矩阵的维数就会急剧的增 列 2 2 2 mk玛局 ； 0 3 3 3 万一蚝一k憋 列一蜀翰捌一翰 riii叫纠il r - ，一耻 郑州大学硕十学位论文总剐矩阵的特点和存储方式 加。但是，随着现代计算机技术的不断发展，运算精度和速度的不断提高，求解 一些超大规模的方程组己成为现实； 2 、对称性，k 1 = k ，即在总刚矩阵中非零元素是关于主对角线对称的； 3 、稀疏性，连续体离散为有限个单元体时，每个结的相关单元只是围绕在 该结点周围为数甚少的几个，一个结点通过相关单元与之发生关系的相关结点也 只是它周围的少数几个，因此虽然总体单元数和结点数很多，刚度矩阵的阶数很 高，但刚度系数中非零元素却很少，一般的一行不会超过2 0 个非零元素。这样， 矩阵的维数越大，则非零元素的比例就越小，稀疏性就会更加的明显。 4 、带状性，刚度矩阵中所有的非零元素都集中在主对角线附近，呈现带状分 布，采用一定的单元编号规则后，使得邻近单元的结点编号相差不大时，带状特 性就更加的明显了，如图2 2 。 5 、奇异性，必须加入边界条件，才能消除奇异性，否则所要求解的方程组 无解。 图2 2 总刚矩阵示意图 f i g2 2s c h e m a t i cg l o b a lm a t r i x 因为有限元求解过程中，大部分时间都是用来求解线性代数方程组，所以我 们了解总体刚阵的上述特性对于提高求解效率是非常重要的。 1 4 郑州大学硕士学位论文总刚矩阵的特点和存储方式 2 5 总纲矩阵的存储方法 目前对总纲矩阵的存储主要是考虑到其稀疏、对称、呈带状分布等特点，已 经发展了多种压缩存储方法。利用矩阵非零元分布的带状性来存储带宽范围以内 的非零元，这种存储方法结构比较简单，因此现在有大量的有限元分析程序采用 此类方法对总纲矩阵进行存储。但是其受到单元剖分网格数目、节点编号、单元 编号以及节点间的连接等的限制。 采用直接法求解时，三角分解总纲矩阵所产生的非零元总是落在带宽范围 之内。因此对于结构简单、规模不大的，求解有限元方程时使用这种存储结构是 合适的。 而对于迭代解法，主要运算是矩阵向量相乘，存储带宽范围内的零元素是 没有作用的。且其一方面增加存储空间，另一方面零元素参与了与非零元的相乘 运算，即增加了程序的运算量。而注塑成型数值模拟的对象，其规模、复杂程度 以及缩减矩阵带宽的效果是不同的，因此，带宽范围内的零元素总量也是不同的 1 7 1 。对于矩阵阶数小、带宽比较窄、带宽内零元素比例较少的，采用这些传统的 存储结构进行程序的迭代求解是可取的。对于矩阵阶数较大的、零元素比较多的、 带宽比较宽的就要采用专门针对大型稀疏矩阵而设计的存储方式【墙】。 2 5 1 二维等带宽存储 对于n 阶的系数矩阵，求出其每一行的半带宽，然后求得所有带宽中最大 的半带宽d 为带宽，则上三角阵中的全部非零元素都将包括在这条以主对角元 素为一边的一条等宽带中。二维等带宽存储就是将这样一条带中的元素，以二维 数组的形式存储在计算机中，二维数组的界是n x d 。若将元素原来的行、列码记 为i 、_ ，它在二维数组中新的行、列码记为i 、_ ，+ ，则有 。= 。 ( 2 3 ) _ ，= j i + 1 采用二维等带宽存储，消除了最大带宽以外的全部零元素，较之于全部下三角阵 大大节省了内存。但是由于取最大带宽为存储范围，因此它不能排除在带宽范围 内的零元素。当系数矩阵大带宽变化不大时，采用二维等带宽存储是合适的，求 郑州丈学硕士学位论文总刚矩阵的特点和存储方式 解也是方便的。但当出现局部带宽特别大的情况时，采用二维等带宽存储时将由 于局部带宽过大而使整体系数矩阵的存储大大增加。此时采用这种存储方法就失 去了意义。如果对带状稀疏矩阵进行直接分解，不会在原带状结构之外引入新的 非零元，因此利用等带宽存储将非常合适。 2 5 2 一维变带宽存储 由于总刚矩阵的每行带宽均不一样，每一行内的元素个数也并不相同。因此 按照等带宽的存储方法进行存储的话，所需要的存储空间还是比较大的，特别是 总纲矩阵中带宽变化较大时，这种现象更加的明显。而一维变带宽存储就是将变 化的带宽内的元素按一定的顺序存储在一维数组中。由于其没有按照最大带宽存 储，因此比二维等带宽存储节省更多的内存。按照解法可分为按行一维变带宽压 缩存储和按列一维变带宽压缩存储。 按列一维变带宽压缩存储是按列依次存储元素，每