（工程力学专业论文）基于响应面方法的二维连续体形状优化理论及程序开发.pdf

摘璺摘要本文将响应面方法引入结构优化设计，解决了二维连续体结构形状优化问题并进行了程序实现。利用响应面方法构造目标和约束函数的近似函数，将原问题中的隐函数显式化并构造二次规划近似模型。通过反复构造近似模型并求解序列二次规划问题，求解结构在位移和应力约束条件下的最优形状。利用响应面方法构造近似模型无需进行灵敏度分析，适于解决灵敏度难求的情况。响应面的构造不涉及结构分析具体过程，可以解决多种类型问题及其耦合问题。然而由于采用最小二乘原理，响应面不能达到精确最优点；此外其整体逼近效果只能利用统计指标模糊地指出，无法构造合理的运动极限，造成近似模型失真。针对以上问题，本文对响应面方法提出改进，构造通过精确点的响应面，结合适当的试验设计方法，在减少试验点的同时保留了约束函数的基本曲率信息，并利用该信息构造二次评价函数。通过线性响应面对评价函数的相对误差确定具有给定逼近精度的响应面范围，形成具有自适应能力的理性运动极限。文中算例对改进响应面和理性运动极限进行了分析，说明了这种优化策略能够明显提高优化质量，验证了它的有效性和稳定性。在程序设计方面，本文实现了基于视窗系统的二维连续体形状优化程序。该程序采用利用v i s u a lc 抖语言开发程序主体结构和图形界面，利用m a t l a bc c + + 数学函数库进行科学计算，利用s t r a n d 7 有限元分析系统进行结构分析。程序能够完成应力和位移约束下的形状优化问题，能够解决带有热载荷的多工况优化问题。关键词形状优化；响应面方法：运动极限：二次规划t h i sa r t i c l ea p p l i e dr e s p o n s es u r f a c em e t h o d o l o g y ( r s m ) t os t r u c t u r a lo p t i m i z a t i o nw i t hc o d i n gp r a c t i c ei ns h a p eo p t i m i z a t i o no ft w o d i m e n s i o n a lc o n t i n u u ms t r u c t u r e s b a s e do nr s m ，t h eo b j e c t i v ef u n c t i o na n dc o n s t r a i n tc o n d i t i o n sa r ea p p r o x i m a t e l ye x p l i c i t e dt oc o n s t r u c tq u a d r a t i cp r o g r a m m i n g ( q p )m o d e l s e q u e n t i a lq u a d r a t i cp r o g r a m m i n g ( s q p ) m o d e lw a sb u i l tb yr e p e a t e d l yc o n s t r u c t i n gr e s p o n s es u r f a c ea p p r o x i m a t i o na n ds o l v e dt of i n do p t i m a ls h a p eu n d e rd i s p l a c e m e n ta n ds t r e s sc o n s t r a i n t s t oc o n s t r u c ta p p r o x i m a t em o d e lw i t hr s m ，n os e n s i t i v i t ya n a l y s i sw a sr e q u i r e d ，t h u sa p p l i c a b l ef o rp r o b l e m sw i t hs e n s i t i v i t yd i f f i c u l t y r e s p o n s es u r f a c ec o n s t r u c t i o ni n v o t v e dn oi n f o r m a t i o ni n s i d es t r u c t u r a ia n a l y s i sp r o c e d u r e h e n c ea b l et ob ea p p l i e dt oav a r i e t yo fp r o b l e m sa n dt h e i rm i x t u r e h o w e v e r , r e s p o n s es u r f a c ea p p r o x i m a t i o nh a dd i f f i c u l t yt of i n dp r e c i s eo p t i m ad u et oi t sl e a s ts q u a r es o l u t i o n b e s i d e s ，i t sa p p r o a c he f f e c tc o u l do n l yb ee v a l u a t e db yo b s c u r es t a t i s t i c a le s t i m a t e s ，w h i c hm a d ei ti m p o s s i b l et oc o m p u t er e a s o n a b l em o v el i m i t sf o ra p p r o x i m a t em o d e la n dc o n s e q u e n t i a l l yl e dt ot h el o s eo fa c c u r a c y t oo v e r c o m et h e s ed i s a d v a n t a g e s ，a t li m p r o v e dr e s p o n s es u r f a c ew i t he x a c tv a l u eo nt h ed e s i g np o i n tw a sc o n s t r u c t e d b yu s i n gp r o p e rd e s i g no fe x p e r i m e n t s ，c o m p u t a t i o n a le f f o r t sw e r er e d u c e d ，w h i l eb a s i cc u r v a t u r ei n f o r m a t i o no ft h eo r i g i n a lc o n s 打m n t sw a sp r e s e r v e d a c c o r d i n gt ot h i si n f o r m a t i o n ，as e c o n d o r d e re v a l u a t i o nf u n c t i o nw a ss e tu pf o re a c hc o n s t r a i n t r a t i o n a lm o v el i m i t sf o re a c hd e s i g nv a r i a b l ew e r ec a l c u l a t e di nt e r m so ft h er e l a t i v ee r r o r sb e t w e e nc o n s t r a i n tr e s p o n s es u r f a c e sa n de v a l u a t i o nf u n c t i o n s e x a m p l e sw e r eg i v e nt oa n a l y s i st h ei m p r o v e dr e s p o n s es u r f a c ea n dr a t i o n a lm o v el i m i t s e v i d e n c ew a sf o u n dt h a tt h en e wm e t h o dw a se f f e c t i v e ，s t a b l ea n dc o n t r i b u t i v et ot h eq u a l i t yo f o p t i m i z a t i o np r o c e d u r e i na d d i t i o nt ot h et h e o r e t i c a lw o r k ，t h i sa r t i c l ea l s oa c h i e v e ds o f t w a r ei m p l e m e n t a t i o nf o rt w o d i m e n s i o n a lc o n t i n u u ms t r u c t u r es h a p eo p t i m i z a t i o n t h ep r o g r a ma d o p t e dv i s u a lc + + f o rw i n d o w sa p p l i c a t i o nf r a m e w o r k m a t l a bc c + +m a t h e m a t i c a ll i b r a r yf o rs c i e n t i f i cc o m p u t i n g ，a n ds t r a n d 7f i n i t ee l e m e n ta n a l y s i ss y s t e mf o rs t r u c t u r a la n a l y s i s t h ep r o g r a mw a sp r o v e dc o m p e t e n tf o rs h a p eo p t i m i z a t i o nu n d e rd i s p l a c e m e n ta n ds t r e s sc o n s t r a i n t s ，a n da p p l i c a b l et ot h e r m a ll o a dp r o b l e m sw i t hm u l t i p l e1 0 a d c a s e s k e y w o r d ss h a p eo p t i m i z a t i o n ；r e s p o n s es u r f a c em e t h o d o l o g y ；m o v el i m i t s ；q u a d r a t i cp r o g r a m m i n g独创性声明本人声明所呈交的论文是我个人在导师指导下进行的研究工作及取得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得北京工业大学或其它教育机构的学位或证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示了谢意。关于论文使用授权的说明兰! 堕：弓本人完全了解北京工业大学有关保留、使用学位论文的规定，即：学校有权保留送交论文的复印件，允许论文被查阅和借阅；学校可以公布论文的全部或部分内容，可以采用影印、缩印或其他复制手段保存论文。( 保密的论文在解密后应遵守此规定)签名：蝉导师签名徘旦麴i皇芝正：弓第l 章绪论第1 章绪论1 1 引言随着科技的进步，工程结构日益复杂。原始的结构设计方法已经不能满足工业发展的需要。在人们对结构安全性和经济性要求越来越苛刻的时候，优化设计成为神重要的手段，帮助人们实现灵活、高效、准确的设计过程。结构优化的发展经历了从不建立优化模型到建立优化模型的过程。优化模型的建立使得高效的数学规划方法得以应用。响应面方法作为近年来兴起的一种建模方法，受到普遍的关注。简便和灵活的特点使其能够被广泛地应用于传统方法难以进入的领域。形状优化是结构优化的主要研究领域之一，具有重要的实际意义。本文尝试利用响应面方法解决形状优化问题中出现的一些难点，对响应面方法进行了改造，同时结合试验设计方法和运动极限的改进探讨了它在形状优化过程中的应用策略。1 2 结构优化概述1 2 1 结构优化的意义结构设计作为一门综合决策的科学l t l 可以根据已有的条件最佳地确定结构的参数以满足一定的工程要求。满足同样要求的设计不是唯一的，各种设计之间有优劣之分。最优设计方案可以尽可能地减少制造成本和提高安全性能。因此寻找最优的设计方案成为很重要的过程。与依赖于工程人员的头脑和经验的传统设计模式不同，结构优化设计将结构分析和改进过程统一在完整的数学模型下，通过数学规划方法寻找满足约束条件的最优方案。数学模型和优化理论的引入使得结构设计的效率极大提高，并且保证了结构的安全性和经济性。随着计算机技术的快速发展，结构优化已经可以应用于大规模复杂结构中，并在航空航天、交通、建筑和机械领域发挥着十分重要的作用。1 2 2 结构优化的要素和研究领域结构优化问题具有目标、约束条件和设计变量三个基本要素。目标是优化设计所希望达到的目的，例如重量最大、刚度最大、造价最低等。一个优化问题可以只有一个目标，也可以有多个目标。多目标优化问题在一定条件下可以转化为单目标优化问题。在追求目标最大或最小的同对，多数工程问题要求结构必须满足一定的约束条件，如应力约束条件( 结构不能发生破坏) ，位移约束条件( 限制结构的变形量) ，频率约束条件( 避开频率禁区) 等等。设计变量是用以控制设计方案变化的参数可以是结构尺寸，如梁的截面积、板壳的厚度等；也可以是形状参数如圆弧的曲率、边界点的坐标等：还可以是结构的拓扑形式。在传统的结构设计中，没有明确划分设计变量的类型。设计的改进往往同时涉及尺寸、形状和拓扑形式的变化。然而在结构优化中，这三类问题有着截然不同的处理方法，须分别加以讨论。因此在结构优化中存在三个主要的研究领域，即尺寸优化、形状优化和拓扑优化。1 2 3 结构优化的基本过程无论是何种类型的结构优化问题，都遵循着一些基本的方法。优化的过程首先从结构模型的建立开始，常用的方法是利用c a d 技术建立数值模型，并对其划分网格建立有限元模型或边界元模型。之后在数值模型上定义设计变量并给出约束条件，建立优化模型。在优化模型中，约束函数往往是没计变量的隐函数，难以直接求解，因此需要对原始的优化模型进行简化。通过将隐函数近似显式化，建立近似优化模型，然后利用线性规划或其它基于导数的数学规划方法进行求解，得到近似最优解。由于近似模型的逼近一般不是全局有效的，因此这个解不是真实的最优解。这时需要反复构造近似模型并求解，逐步逼近真实最优解。另外一种方法不建立近似优化模型，采用随机试探的方法在原始模型上直接寻优。虽然这种方法可用于寻找全局最优解，但是它迭代次数过多，结构分析计算量过大，在实际应用中受到很大限制。1 3 形状优化概述1 3 1 形状优化的两种类型形状优化是结构优化的重要组成部分，在工业生产中有大量的需求。形状优化以表达结构外形的参数为设计变量。根据所处理的结构类型，可以分为骨架类结构形状优化和连续体结构形状优化。骨架类结构形状优化通过改变桁架的节点位置和框架的几何尺寸到达优化目的，多以节点位置或杆长为设计变量。这类问题的形状表达方法简单，设计变量定义明确，数量少，结构分析时计算量小，是较好解决的一类问题。在高压输电塔架和桁架建筑结构的设计中经常能够见到这类问题的应用。连续体结构形状优化通过改变结构的边界形状寻找最优结构，多以边界曲线参数、几何形状参数或边界节点坐标作为设计变量。连续体形状优化问题较骨架类形状优化问题困难。首先，形状的表达十分复杂，准确性和设计变量的个数常常是一对矛盾。其次，形状改变后需要相应地调整有限元或边界元网格。特别是有限元网格，不但要调整边界上的网格，还要调整结构内部的网格。网格发生变化后可能导致有限元单元发生畸形，导致分析结果精度严萤下降。最后，形状优化的设计变量形式多样，设计变量和目标函数以及约束条件之间的函数关系比较复杂，难以如尺寸优化那样根据结构特性采用倒变量或其它映射方法进行简化，属于高度非线性问题。由于这些问题的存在，连续体结构的形状优化一直是较难解决的问题，也是结构优化研究的热门问题。1 3 2 连续体结构形状优化的历史连续体结构形状优化问题可以追溯到2 0 世纪中叶。1 9 7 3 年z i e n k i e w i c z 和c a m p b e l l 2 】研究了水坝的优化问题。他们以有限元节点坐标作为设计变量，使用序列线性规划得到水坝的最优形状。此后大量的研究相继展开，并且在形状表述、网格变动、有限元分析、灵敏度分析和数学寻优方法等几个领域取得了成果。在形状表述方面，1 9 7 5 年v a n d e r p l a a t s 3 使用节点坐标作为设计变量，使用近似技术建立设计变量与结构外形之间的关系，从而减小了结构分析的计算量。但是计算量过大是以节点坐标作为设计变量方法的固有缺陷，因此人们开始寻求其它办法。1 9 7 6 年k r i s t e n s e n 和m a d s e n 4 l 利用正交函数的线性组合描述边界形状。1 9 8 4 年p r a s a d 和e m e r s o n 5 l 也做了类似的工作。多项式的采用减小了问题规模，但同时也引入了形状约束，使结构只能在固定模式下改变形状。另外，第1 带绪论使用高阶多项式容易出现边界的震荡现象。进入8 0 年代，l u c h il 6j 干hw e e k - 1 等人使用三次样条曲线代替多项式。b r a i b a n t 和f l e u r y 8 9 l 则使用b 样条和b e z i e r曲线描述边界形状。样条曲线具有照好的光滑性，而且使用b 样条曲线还能方便的得到灵敏度的解析表达式。1 9 8 2 年i r i l a m 【lo 】提出了设计单元的概念。设计单元是结构的若干子区域。它的形状由其边缘上的关键点控制，而关键点的位置由设计变量决定。设计单元内包含若干有限元单元和节点，通过形函数映射关系来改变有限元网格形状。变形后网格只是形状上发生变化，而拓扑结构没有发生变化。这使得总刚度矩阵不变，有利于减少结构重分析的计算量。在网格变动方面，b e n n e t t 和b o t k i n j ，s c h i e m e i e r 和s z a b o 【1 2 1 ，以及k i k u e h i 1 3 1 等人研究了自适应网格重划分方法。这些方法采用增加额外单元、提高原有单元阶数和改变节点位置来重新划分网格。其间许多衡量网格划分质量的方法被提出，其中最著名的是能量误差变化率法则，即计算应变能密度的期望值与计算值的差值( 即误差) ，并进行规格化处理，据此确定网格改进的区域并进行网格的自动改进。1 9 9 7 年王各、隋允康【1 4 j 提出了二级控制法，在第一级控制中，自然设计变量通过边界形状函数来决定关键点的坐标；在第二级中，关键点通过网格节点在设计变量中的参数坐标来决定网格节点的坐标。在有限元分析方面，n o o r 和l o w d e r i ”1 以及k i r s c h t ”】使用t a y l o r 级数展开式近似节点位移，以减少计算量。该方法虽然可以节省计算时间，但是精度较低。隋允康【1 7 】提出使用解析的方法，利用计算机符号计算技术得到结构的响应数学表达式，从而避免在迭代过程中使用昂贵的结构分析。y a n g 【18 】使用了类似方法对曲柄等小规模简单问题进行结构优化。该方法将结构重分析规模过大的难题转移到了形式计算问题上，但在实施中还有诸多困难，目前尚处于研究的初步阶段。灵敏度分析是结构优化中最重要的部分之一，这方面的研究十分活跃。k 五s t e n s e n 4 1 ，w a n 9 1 1 9 ，b r a i b a n t l 8 ，9 2 ，c h o il 2 1 , 2 2 1 ，l i e f o o b e ，s h y 3 和f l e u r y 2 3 l等人以及其他许多学者都做了深入研究。灵敏度的求解基本可以分为两种方法，其一是连续方程微分法b ，其二是离散系统微分法。连续方程微分法使用本质导数的概念，对连续方程微分，再对边界积分，从而获得灵敏度，最后才将结构离散为有限元模型。离散系统微分法则先将结构离散，然后以有限元方程为出发点求导数。以位移约束灵敏度分析为例，将有限元方程置u = p( 1 1 )式中置结构刚度矩阵；，节点位移向量；j p 体力载荷向量。左右两端对设计变量求导数( 1 2 )式中x ，设计变量第i 个分量。整理后得到位移对设计变量的导数孚：f 罢一篓u 1 ( 1 - 3 )眺l 眠缸。j离散系统微分法具有简单易行的优点，是目前比较流行的做法。竺挑=u堡钆到翌瓠得n北京 业人学r 学坝，【学位论文若以解析手段求解公式( 1 3 ) 中的芸和，称为灵敏度分机的解析法：若e x d x j以差分手段求解，则称为差分法。解析法通过数学推导得到目标和约束函数对设计变量灵敏度的准确表达式。通过解析法可以求得设计点上的准确导数值，但是在很多情况下，解析表达式的建立非常困难。因此考虑使用差分法获得近似导数。差分法可对任意函数求灵敏度，但是精度受差分步长的影响较大，并且会出现差分步长过小而导致不能求解的困难。另外差分法的计算量随着设计变量规模增长很快，难以应用于大规模问题。将以上两种方法结合起来，就形成了半解析法。半解析法可以有效地解决解析灵敏度难求的问题，但是它仍然存在差分法不精确和解析法通用性差的缺陷。在求解方法方面，1 9 6 0 年s c h m i t 【”j 将线性规划引入结构优化。1 9 7 9 年f l e u r y i ”】提出了凸规划的对偶方法，并指出优化准则法是数学规划的特例。隋允康【i l 从方法论的角度建立了关系映射反演理论，利用该理论推广了数学规划中的一些基本原理，并引入结构优化领域。他利用该理论对线性规划、二次规划、广义几何规划和对偶规划进行了改进，提出了更加有效的序列映射求解方法。在传统算法之外，数学界的最新优化理论研究成果也被引入结构优化中，如模糊优化理论、神经元理论和遗传算法等。在其它方面，b o t k i n 研究了多工况问题 27 】和子结构法在三维实体形状优化中的应用【2 8 】。k u h n 和t u c k e r 最早研究了多目标问题。张卫红 2 9 , 3 0 研究了运用加权法解决多目标优化问题，提出多目标优化的灵敏度分析技术。为处理大规模优化问题，a z a r r n 和l 1 改进了由s o b i e s k i 提出的二级优化方法，提出了更完善的灵敏度分析方法。1 3 _ 3 连续体结构形状优化的新进展近些年连续体结构形状优化出现了许多新进展。首先，形状优化的研究对象从弹性结构转向弹塑性结构、高塑性结构。涉及的问题也从简单的静力问题转向热力耦合、流体固体耦合、接触、疲劳等复杂问题。这些问题的研究仍然以灵敏度分析为重点。顾元宪i 2 】研究了热应力约束下的形状灵敏度分析。l u n d ”】研究了大位移下稳态流体固体接触问题。k j m ”】等人研究了高塑性结构接触问题的灵敏度分析。l i n t ”j 在其博士论文中详细讨论了非线性结构在疲劳载筒下的形状优化问题。其次，研究的手段也发生了变化，不再完全依赖有限元作为主要的计算方法，而是发展了基于边界元和无网格法的形状优化理论。边界元法只在结构边界上划分单元，并将域内积分转移到边界上，避免了优化过程中有限元网格发生畸形的问题，同时对简化形状描述和设计变量定义有所帮助。近几年快速多极算法【3 6 1 的引入消除了边界元法计算效率低的缺陷，使得这种方法可以被应用于大规模问题。l e e 和k w a k t ”j 给出了基于边界元法的二维热力问题的形状灵敏度分析。无网格法是近年兴起的结构分析方法。它采用节点周围范围可变的域代替有限元网格，同样解决了网格重划分的问题，非常适合于大变形结构分析问题和形状优化问题。b o t k i n l 3 8 和k i m 3 9 , 4 0 1 等人在此方面进行了研究。第l 章绪论1 4 本课题的目的和内容1 4 1 课题目的如前文所述形状优化需要建立近似模型。近似模型的建立往往通过结构灵敏度分析得到目标或约束函数的导数值，然后利用泰勒展式构造一阶或二阶近似函数。灵敏度分析是建立近似优化模型的关键步骤。解析法灵敏度分析可以得到当前设计点的精确导数，从而构造较高精度的近似函数，但是解析法需要对结构的状态方程或有限元方程进行求导，过程繁琐，计算量大，并且对于不同问题需要建立不同的灵敏度表达式，缺乏灵活性。对于弹塑性、热传导等经典问题，灵敏度解析表达式的建立比较容易，而对于摩擦、粘弹性、接触和多场耦合等复杂问题，解析灵敏度表达式的建立是十分困难甚至不可能的。在解析灵敏度不可求的情况下，可以采用差分方法求解设计点的近似导数值。然而差分法的精度受到差分步长影响很大，若要提高精度，需减小步长，但是步长过小将会导致截断误差起主导作用，同样会降低精度，甚至出现溢出错误。为了解决灵敏度分析的种种困难，r o u x 4 1 1 和s h y y l 42 j 等学者将响应面方法 4 3 】引入结构优化设计作为构造目标或约束近似函数的方法。利用响应面方法建立优化模型可以避开解析灵敏度分析的问题，同时也不存在差分步长困难。而且响应面方法可以通过同一组试验点集合得到多组响应值，特别适合于处理带有多个约束或多种约束的结构优化问题。然而响应面方法也存在诸如不能表达精确点、所需试验点过多和无法准确评价有效范围等缺陷。许多学者对响应面方法进行了改进，并探讨了它在结构优化中的应用i 4 4 。“j 。这些研究的焦点大多集中在响应面的构造和试验设计方法的改进上，然而在结构优化问题中，响应面构造、试验设计方法和响应面评价是密切相关的问题。它t 1 的改造和相互协调的问题，即优化策略问题，同样值得关注。本文希望通过研究寻找响应面方法在形状优化中的应用方法，并针对响应面方法的缺陷进行改进，同时研究运动极限在优化过程中发挥的作用，探索更好的运动极限计算方法，寻找合理、有效而稳定的优化策略。1 4 2 课题内容本课题涉及理论推导、程序设计和数值实验分析方面三个方面的工作。在理论推导方面，本文分析了原始响应面的缺陷，并针对这些缺陷提出了改进，推导了改进后的响应面公式。同时研究了中心复合、单纯形、拉丁超立方等试验设计方法，并予以结合和拓展。研究了响应面逼近的有效范围，分析了传统响应面评价方法的缺陷。在改进响应面的基础上提出合理的响应面评价方法，并推导了根据此方法计算运动极限的公式。在程序设计方面，本课题开发了二维平面结构优化程序p l a t e s h a p e r 。程序拥有完整的前处理功能，能够直接在c a d 模型上定义设计变量、划分设计单元并施加约束。程序界面设计符合规范，内容丰富，操作灵活。优化算法采用模块化设计，参数设置灵活，并提供多种算法，可供比较。程序具有优化结果评价功能，可以对优化结果、响应面逼近效果和运动极限做出评价。详细的优化结果数据以e x c e l 表格形式出现，方便查阅。程序还提供动画功能，可以显示优化时形状的变化过程。程序整合了s t r a n d 7 有限元分析系统、m a t l a b 函数库和e x c e l 表格，具有强大的计算、分析和数据处理功能。程序的丌发符合专北京r 业人学t 学f ! l 学位论文业丌发标准，采用面向对象技术，层次划分清晰，具有良好的可扩展性。在数值实验分析方面，本课题利用所开发的优化程序对应力、位移约束下的多种结构进行了分析，验证了改进响应面方法及相应运动极限计算方法的效果。对带有自重、热载荷等灵敏度分析困难结构进行了优化，并针对实际问题进行了多工况下的形状优化。1 5 开发平台简介二维连续体形状优化软件的开发采用v i s u a lc + + 语言，构造程序框架，采用m a t l a bc c + + 数学函数编写数值计算模块，采用s t r a n d 7 有限元分析系统作为结构分析引擎。v i s u a lc + + 是当前最流行的面向对象程序设计语言之一。它全面支持面向对象技术，使用灵活，功能强大，能够容易地实现程序间通讯，非常适于做软件的二次开发工作。此外它提供图形化工具和代码自动化工具，配合使用微软基础类库( m i c r o s o f tf o u n d a t i o nc l a s s e s ) 使开发者不必在程序界面上花费过多的精力。虽然c + + 语言功能强大，但它并非针对科学计算而设计，因此在数学计算方面有很多缺陷。对于响应面方法中所需的复杂矩阵运算和数学规划程序，使用c + + 语言编写费时费力，且难以保证计算质量。m a t l a b 是m a t h w o r k s 公司开发的功能强大、内容丰富的数学计算软件，具有丰富的数学函数和高质量的矩阵计算算法，可以作为计算引擎，弥补c + + 语言缺陷。有限元结构分析是形状优化过程中重要的组成部分。本课题采用澳大利亚g + dc o m p u t i n g 公司开发的通用有限元分析程序系统s t r a n d 7 作为有限元分析引擎。s t r a n d 7 可以广泛地应用于航空、土木、机械、造船、结构、岩土等领域。它能够求解线性和非线性静力和瞬态动力问题，能够进行线性曲屈分析，振动分析，简谐响应分析，谱响应分析，以及线性和非线性稳态和瞬态热传导分析。s t r a n d 7 有一套应用程序接口函数，允许开发者在其它应用程序中使用它所提供的功能。开发者可以在应用程序中提取s t r a n d 7 模型的几何数据、有限元数据和分析结果，同时也能直接改变模型中的参数。开发者还可以调用s t r a n d 7 求解器求解自己的有限元模型。该接口能够支持d e l p h i 、c c + + 、b a s i c 、p a s c a l和f o r t r a n 等多种语言，为二次开发提供了有力的支持。第2 节7 够状描述和州格控制第2 章形状描述和网格控制2 1 形状描述和设计变量定义2 1 1 形状描述为了优化结构形状，首先需给出对结构边界形状的描述。如第1 章所述，形状描述的方法多种多样，但是它们均和设计变量的定义有关。不合理的描述不但引入额外约束，还导致过大的设计变量数目。在建立近似优化模型时，计算量与设计变量个数息息相关，因此在选择形状描述方法时应当结合近似模型的建立方法，尽可能减少计算量。样条曲线和分段函数是常用的描述边界曲线的方法。样条曲线可以表达较复杂的曲线并且保证曲线的光滑性。但是样条曲线的形状往往通过曲线参数或位于曲线外的控制点控制，这使得设计变量的定义脱离直观几何意义。另外曲线参数或控制点坐标变量数目较多，无益于减少设计变量个数。分段函数的控制点少，且均位于曲线上，具有直观的几何意义。这一点可以为利用自然设计变量进一步减少设计变量数目提供条件。采用分段二次函数，尽管不能保证边界的光滑性，且具有一定的逼近误差，但能够满足一般的形状要求。由于本文主要在于讨论优化问题而非形状表达问题，因此拟结合网格控制使用分段二次函数表达结构形状。2 1 2 自然设计变量若将结构边界曲线分为”段二次函数，则控制点个数为2 n ( 曲线封闭) 或2 n + 1 ( 曲线不封闭) 。当结构形状较为复杂时，控制点可能会很多。若赢接选取控制点坐标作为设计变量，设计变量数目过多。在实际应用中，结构的设计往往需要满足多方面的准则，因此不能让结构形状随意变化，而应当遵循一定方式。结构的形状必须首先满足一定的几何功能，因此工程人员倾向于根据其几何特征改变形状。几何特征是用来描述几何结构的自然参数，比如长、宽、高、角度、直径、曲率半径等等。描述结构形状的最简便方法就是利用这些参数作为设计变量，称为自然设计变量( 为方便起见，自下节起简称设计变量) 。由于结构的边界形状由分段二次曲线表达，因此必须通过一些函数将自然设计变量与二次曲线的控制点联系起来，称这些函数为控制函数。由于自然设计变量个数一般较分段二次曲线的控制点少很多，因此自然设计变量和控制函数的引入起到了减少设计变量数目的作用。同时控制函数为用户指定的函数关系，其形式灵活，为添加设计变量之间的附加约束关系提供了方便。2 1 3 控制函数控制函数是指定分段二次曲线的控制点与设计变量之间关系的函数。在平面上，每个控制点具有两个坐标，每个坐标均需要一个控制函数。若定义设计变量为x ee ”，关键点为异( 口。，肛) “= 1 ，m ) ，则第i 个关键点的控制函数为卜。4 ( 。)协1 1i 声，= e ( x )、北京l 业人掌t 学坝i 学位论义式中a ，只关于工的任意函数，其形式由用户根据关键点和没汁变量的关系给出，具体操作详见2 3 4 节内容。2 _ 2 网格控制2 2 1 网格的变形由于采用有限元方法进行结构分析，当结构形状发生改变后，需改变有限元网格。有限元网格的改变主要依赖于网格映射方法和网格重划分技术。网格重划分在每次形状发生改变后根据新的几何体重新划分网格，出现与原有网格节点数量、单元数量甚至形式完全不同的新有限元模型。虽然重划分可以保证网格的质量，但是每次形状改变都要重新计算所有网格，计算量过大。另外，当网格的拓扑结构变化后结构的刚度矩阵亦发生变化，这给依靠刚度矩阵进行的敏度分析带来影响。因此最好采用映射方法调整网格形状，并对由于变形而产生畸形的网格予以调整以保证网格质量。2 2 2 二级控制的基本原理二级控制是一种通过映射方法调整有限元网格的技术。其基本思路为：首先，在结构上定义一些描述形状的关键点( 如二次曲线的控制点) ，然后利用关键点将结构划分成多个子块，称为设计单元。设计单元内部包含有限元节点。在初始形状上建立有限元节点与设计单元的映射关系。当设计变量发生变化时，通过控制函数改变关键点位置，设计单元的形状随之发生改变。之后，通过先前建立的映射关系修改设计单元内部的有限元节点坐标，从而达到调整有限元网格的目的。其数学原理可以描述如下。第一级控制：若设计变量为x e “，第i 个可动关键点坐标为口。和卢，可以表示为设计变量的函数i q = 4 ( x )【店= b j ( 石)第二级控制：若有限元节点的坐标是( j ，y j ) ，有限元节点坐标可以表示为，z ，= m 如) = ( f ，咖。1 7y ，= m 舻) = f ( f ，即) 层f = l、jj 7设计变量图2 - i 二级控制f i g u r e2 - 1t w ol e v e lc o n t r o l( 2 2 )并位于设计单元k 上，则式中i 设计单元k 中关键点的总数；m ( 手，7 7 ) 内部序号为i 的关键点在设计单元k 中的形函数孝，7 7 有限元节点在设计单元k 中的参数坐标。( 2 - 3 )2 2 3 关键点和设计单元的定义关键点是结构形状的基本控制点，往往定义在能够表达结构形状特征的位置上，如边界、角点、圆心等。当结构的某些部分不需要发生变化时，这些部位的关键点可以定义为固定关键点，而其它关键点定义为可动关键点。关键点除了包括结构边界上的二次曲线控制点外，还包括结构内部适当位置的一些点。这些点虽然不表达结构外形，但是参与构造设计单元。设计单元是具有规则形状的子结构。设计单元的构造采用与有限元单元类似的方法。利用关键点的位置通过形函数插值计算其上各处的变形量。为了便于实施，采用8 关键点( 1 = 8 ) 的四边形单元作为设计单元，正好符合使用分段二次曲线表达形状边界的要求。设计单元的形式如图2 2所示。实心圆点代表关键点，关键点下方的序号是该关键点在设计单元上的内部序号，括号内的坐标是该点在设计单元内的参数坐标。设计单元的形函数采用s e r e n d i p i t y 族，具有如下形式月，t ，1 )5 ( 口，1 )b ( 1 ，1l ( 1 ，口)：0 ，口)；( i ，厅1 ( 1 。1 )i ( o 1 )( 1 ，1图2 - 2 设计单元f i g u r e2 - 2d e s i g ne l e m e n tn i = ( 1 一掌) ( 1 一叩) ( 一f r - 1 ) 4n 2 = ( i f ! ) ( 1 一呷) 23 = ( 1 + f ) ( 1 7 7 ) ( f 一叩一1 ) 4a = ( 1 + f ) ( 1 一叩2 ) 2 ( 2 - 4 )5 = ( 1 + f ) ( 1 + 7 7 ) ( 乎+ 7 7 - 1 ) 4n 6 = ( 1 一孝二) ( 1 + 叩) 2n 7 = ( 1 一f ) ( 1 + 叩) ( 一掌+ 1 7 1 ) 4 i = ( 1 一掌) ( 1 一叩2 ) 22 2 4 初始网格映射的建立建立设计单元后，将有限元节点的初始参数坐标代入式( 2 3 ) 即可得到有限元节点的坐标。有限元节点的初始参数坐标来自于设计单元建立时有限元节点在设计单元上的参数坐标。因此在建立设计单元之后，需要求解各个有限元节点的初始参数坐标。建立初始网格、关键点和设计单元后，由式( 2 。3 ) 可以反解出有限元节点，在设计单元k 上的参数坐标。将式( 2 4 ) 所示形函数代入式( 2 3 ) ，得到l x = ( c 7 4 2 1 7 + c 6 4 r 2 + c s 孝2 + c 4 r 1 2 + c 3 f 叩+ c 2 f + c 1 ，7 + c o ) 4，1ni y = ( d 7 孝2 r l + d 6 善玎2 + d g 。+ d 4 r 1 2 + d 3 善7 7 + d 2 掌+ d l ，7 + d o ) 4式中c o i ，o ，d 0 ，一，d ，有限元节点所在设计单元上的关键点横、纵坐标的线性组合；若以( a 。屈) ( f = 1 ，8 ) 表示关键点坐标，有北京t 业大学f 掌蝤11 “学位论文= 一口i + 2 口! 一口3 + 2 0 ：4 一口5 + 2 a 6 一口7 + 2 口rq = 一2 a 2 十2 , z 6。2 = 2 a 4 2 三= 。0 口。- - + t 口2 3 ，+ 一0 5 - - 。+ t 2 口+ 2 a ；+ 口，一：口。( z 一6 )0 = 口i + 口3 4 + 口5 + 口7 2 口8、一叫c 5 = 口i 一2 口2 + 口3 + a 5 2 口6 + 口7c 6 = 一口l + 口3 2 口4 + 口5 一口7 + 2 c 7 = 一q + 2 口2 一口3 + o r 5 2 a 6 + 口7d o = 一乒l + 2 p 1 一盆l + 2 $ 一p s + 2 a p 1 + 2 _ b 。d l = 一2 岛+ 2 成d 2 = 2 p , 一2 p 8瓣：会：笏袅+ p 7 蝴，d 。= 届+ 屈一2 尻+ 店+一2 风d 5 = 移t 一2 p 2 + p 3 + ；b s 一2 p , + p 1d 6 = 一0 1 + 8 ，一2 z , + p s 一8 1 + 2 ；b d 1 = 一芦t + 2 p t p l + p s 一2 p b 十母1式( 2 5 ) 可以化为关于f 和叮的非线性方程组 簟叩) = 。j i ：? + c s q r f 。+ c s 善2 7 。t 2 。，掌r 1 + c z 掌+ c t 町+ c o - 4 x = o ( 2 8 )( 手，7 7 ) = d 7 孝2 r + d 6 亭r 2 + d 5 f 2 + d 4 r 2 + 以翻+ d 2 善+ d l r + d o 一4 y = 0 式( 2 8 ) 可以使用牛顿法求解。对于上述问题，牛顿法的迭代公式如下点+ i = 茧+ _ 1j 1 巩+ 万k 1l 工( 磊，仉) =戮( 薯，1 7 )a f瓠皤，仉)a fz ( 厶，仉)五蛾，仉)利用式( 2 - 8 ) 可以求得以上各式中的导数表达式重量! i j 堕：2 c ，孝刁+ 吒7 7 ：+ 2 q f + c 3 j 7 十c ja !”“”2至挈：c ，善z + 2 c 6 善刁+ 2 c , r + c ，+ c 】d 门( 2 9 )( 2 1 0 )( 2 - 1 】)( 2 - 1 2 )( 2 1 3 )鞯2 带啦：描述和嘲格控制掣：2 ( f 1 手_ + d 。叩2 + 2 d 5 善+ d ，q + d ：( 2 - t 4 )a 与掣掣：d ，善! + 2 d 。勃+ 2 d 。，7 + 蟛+ d 。( 2 - 1 5 )u 选定初始值( 岛，7 7 。) ，代入式( 2 - 9 ) 、( 2 - 1 0 ) i 塞代求解，若l 0 ，则可以得到方程( 2 8 ) 的近似解。若两次迭代的解满足收敛准则m a x ( 1 f 。+ 。一 ，1 7 。一1 7 1 ) r a i ns f g f ( x ) s 0j = 1 ，n 。( 3 - 1 )g ：( x ) 0k = l ，心兰：x i ii = 1 ，- 一，”式中s ( x ) 目标函数。取结构重量，对于密度均匀的材料在二维问题中简化为面积；g j ( x ) 应力约束函数；或( x ) 位移约束函数；。x 。，暑设计变量各个分量的下、上限。优化模型给出了形状优化问题的完整描述，但是由于s ( x ) 是设计变量的复杂函数，而g j ( 工) 和戤( x ) 往往是设计变量的隐函数，直接求解该模型十分困难。因此必须将以上函数转化为具有简单形式的近似显函数，将以上模型转化为近似模型，才能利用已有的数学规划方法进行求解。二次规划是目前最有效的规划方法之一，为了采用二次规划作为解决问题的手段，需要将优化模型化为目标函数为二次函数、约束函数为线性函数的二次规划形式球e 45 ( _ ) = x 7 h x 2 + ，7 x _ m i ns t a o x + b os 0( 3 - 2 )a 。x + b 。0 x x x式中a o x + b o 应力约束近似函数；一。x + 玩位移约束近似函数；小融吲驴妒 吲耻l n 荒，“乏j1 j第3 聋优化模型的建口3 1 2 近似化方法将优化模型转化为近似模型需要对隐函数进行显式化近似。形状优化中应力和位移往往都是设计变量的复杂隐函数，非线性程度较高，因此不能在很大的范围上进行近似，而只能在一个较小的局部区域对其进行近似。局部近似最常用的方法是泰勒展开。根据高维空间下的一阶和二阶泰勒展式，在x 印的局部邻域上，函数f ( x ) 可以近似地表示为厂( 工) f ( x o ) + v 1 f ( x o ) ( x x o )( 3 3 )f ( x ) f ( x o ) + v 7 f ( x o ) ( 上一x 哪) +