论文相关方法-基于Kriging模型的桁架结构可靠性分析-CNKI知网查重网站

论文相关方法-基于Kriging模型的桁架结构可靠性分析

2021-05-26 17:11:29
作者:杭州千明

  

  在结构可靠性分析方面,模型法被广泛使用,因为它可以大大减少技术分析中的负担。模型法通常的原理是采用构建如之相近的模型来掌握变量和结果二者间的相互关系,该方法可以代替可靠性研究中的极限函数。本文是基于Kriging和可靠性计算两者方法的整合,设计出了以Kriging方法为基础的MATLAB程序,采用DACE工具箱构建Kriging模型和预测函数,通过预测函数预测各点的响应值,计算结构的可靠度,从而进行结构的可靠性分析。桁架作为一种广泛应用的结构,结构简单高效,优点十分明显。首先,本文简要分析编写3个桁架结构的节点位移的有限元程序,得到了计算功能函数的程序,用最普遍的Monte Carlo模拟求解了其失效概率,随后介绍了结构可靠度的基本概念,介绍了一种半参数化的插值技术——Kriging方法,利用MATLAB的DACE工具箱建立模型并预测功能函数值,结合一种主动学习的方法(AK-MCS),更加高效地求解结构的失效概率。

  进入二十一世纪后,我国经济正在逐步向好发展,综合国力也在日益增强,然而国内虽向好发展,但从国际的角度来看,依旧存在较大的压力,经济在发展的同时存在很大的不确定性。现如今,消费者在购买和使用电子产品时,一方面要追求良好的性能和外观设计,另一方面对产品的可靠性也是日益加剧,所以,中国倘若要在国际上崭露头角,这就要求在相关产品的开发和设计中,严格对产品的可靠性研究进行更为系统的挖掘,只有这样方能使相关产品在可靠性上有一定的成效,同时在设计和开发时严控成本。早在1930年左右,可靠性相关的理论研究就已经产生,它的发展也是较快,此种方法产生的背景是,那个时期没有办法运用理论方法去说明某些产品失效。在20世纪40年代左右,《适航性统计学注释》[1]一书在航空行业深受欢迎,此书的编写单位是英国航空委员会,这本书第一次提出了飞机的故障率应该低于次/h的观点,这个标准后来认定为最早的飞机可靠度衡量指标;在二战期间,美国军队的各种电子设备大规模的发生失效事故,此时美国军队遭到重创,鉴于此,美军便逐对相关设备在可靠性方面有了一定的关注,并开始研究,由此组建了“电子产品可靠性联系小组”,这个组织的常规任务通常是对电子商品及军队器件采取完整的可靠性研究和分析。大约在1960年,可靠性研究相关的期刊《电子设备可靠性报告》已经出现,时间的推移,科技的力量正在壮大,这也将可靠性理论研究带入各行各业,也就是由最开始的电子产品、航空方面等转向机械及土木等行业,而且他们的发展也是相当快速。以此为节点,可靠性研究方法发生了较大的变化,已整体在社会生产中广泛运用,同时给社会生产生活带了极大地便利。另一方面,可靠性研究在结构上和思路方面的发展亦是如此。在上世纪五十年代左右,Freudenthal A.M.就发表了《结构的安全》[2]的文章,首次将结构安全度的思想进行了发表,这为可靠性研究在结构上的研究奠定了基石;70年代时,Cornell等人也是提出了可靠度指标的概念,他们以此作为结构可靠性研究的指标,这种方式为可靠性在结构上的分析有利,以此可以构建二阶矩模型;80年代左右,Rackwitz等人在考查变量数据不是正态分布时,得出可以运用转换的思路来去处理可靠性研究在结构上的分析。正是由于这些奠基人提出的结构可靠性方法,这促使了这种方法在每一个行业的应用越来越广泛。

  实际工程结构具有复杂性,为了保证在各种不确定因素影响下,如结构受载、材料特性不同以及实际安装存在误差等,能够保证结构的安全使用,一般需要专业的可靠性分析,对其进行定量的分析。加之工程结构的大多数功能具有高度的非线性和低失效概率,需要对有限元进行复杂的分析。更广泛的MCS方法需要大量的调用运行,大量的计算,低效率,并且在实际应用方面有局限性。

  1.2国内外研究现状

  因此,目前最有效的可靠性分析,主要有反应面(rsm)方法、矢量支持机、克里金模型和神经网络。与其他方法相比,克里金模型在建立代理模型方面具有更高的精确度,较为繁琐的结构的分析效率也较高。现如今推出了较多的新方法,其目标通常是实现获取某些少许样本点,同时构建一种较为准确的模型,该模型可以增强准确性,在下一次迭代中增加样本点,并逐步提高代理模型的准确性。即根据学习功能选择最佳样本点。相关研究[3]指出一种EGO(Efficient Global Optimization)方法,并用该法处理较为繁琐函数的完善方面,得出的结果是,可运用EI(Efficient Improvement)函数去获取Kriging的最优样本点。根据EGO法,研究[4]推出了EGRA(Efficient Global Reliability Analysis)的方法,并将其应用于求解线性关系不好的隐式状态方程的可靠性研究中,同理得出,若运用EFF函数,其获取的样本点会增加。文献[5]结合Kriging模型和学习函数U,提出了采用AK-MCS方法及AK-MCS方法求解结构的可靠性。Yan Gang Zhao和Tetsuro Ono给出了影响一次二阶矩、二次二阶距方法精度的三个要素的范围,和这两种方法的使用场合和具体使用的步骤[6];A1TIleDDer Kiureghian等人[7]在文中推出结构可靠性研究方法,采用的是抛物面函数处理极限状态函数的二阶近似;李炜和康海贵亦是推出了一种结构可靠性分析的方法,该法也称二次二阶矩法,这种方法可以有效处理曲率计算时存在的困难,另一方面,运用这种思路,还可以构建计算曲率的公式,从而是整个过程变得简单[8];赵维涛等人也是在研究中,得出如果删除泰勒级数展开式中含有的变量的耦合项,把结构功能函数转变为上述矩阵方法,将会使最终的结果准确性得以显著提高[9];20世纪50年代后期,又出现了一种新的方法,即蒙特卡洛数值模拟法,在科技力量的不断促进下,实践工程方面的要求的设计对象十分的繁琐,原始的研究方法不能实现,然而应用上述方法,则可以发挥它的优势,也正是这样,这种方法也越来越受到人们的推崇。上述方法采用的原理和随机抽样原理一致,运用此法进行实验,可以获取对象函数的统计分布参数值,这对很多领域问题的解决都有帮助,此外,这个过程对设计变量几乎无任何其他的要求。然而,该法的原理因为比较简单,在计算机相关的指令下即可达到编程目的,因而在结合数学理论后形成了一系列的方法。如果将该法应用于结构可靠性研究中,也是非常恰当的,我们将其称之为蒙特卡洛可靠性分析方法,这种方法得到的结果准确度高,而且非常简便,因为备受人们的欢迎。在蒙特卡洛方法基础上发展起来的方法有很多,例如重要抽样法[10]、子集模拟法[11]等。S.K.Au等人推行了一种自适应可靠性计算方法,它的原理是重要抽样[12],然而,该法也与马尔科夫模拟算法相关,获得样本点遵循马尔科夫链规则;万越等人在上述方法的基础上将其灵敏度进行优化,得到一种新的方法[13],它的主要意思是采取标准正态分布中无效部分分布在一个球上,该球的球心是坐标原点,如果超出了球的半径,那么只需要分析超出球外的重要抽样样本点相关的数据,以实现灵敏度加强的一种评估方式;袁修开等人也是模拟了重要抽样,设计出比较新颖的可提高灵敏度的方法[14],此法显示,如果在分析准确度一致时,该法会比蒙特卡洛模拟法的分析效率更高;宋述芳等人研究了一种基于线抽样的可靠性灵敏度分析方法[15],这种方法拥有精度高、收敛快且适用于高维及多模式情况等优点;PradlwarterH.J.等人研究了线抽样模拟方法在可靠性问题中的应用情况[16],通过一些具有代表性的实际问题验证了这种方法的高效性和准确性;Song S.F.等人在研究中,得到了一种可优化线抽样的可靠性研究的方法[17],同时还特别分析了此种方法的使用领域,该法的原理是以马尔科夫链在无效域内可以抽样,同时抽样中最关键的方向均由相关样本点控制;S.K.Au等人以高维状态为背景,采取子集模拟法去分析较低失效机率的疑问,得出的结论是该法这种方法的处理效率非常高[18]。最初的基于替代模型的可靠性分析方法主要是多项式响应面法,以及一些它的变换形式,多项式响应面法按照某种指定的实验设计方法进行抽样[19],结果是运用迭代获得一准确度较好的多项式回归模型,可以取代原来的隐式极限状态函数,同时按照此多项式通过响应面的方式进行可靠性分析。在这以后,人们意识到这种方法在实际使用中会受到多方面因素的干扰,造成的结果便是不能很好地满足实践需求,为此,科研工作者便开始其他的可靠性方法,期望能有所改善,在这过程中,人们也发现了一些类似的模型,于是便将其进行充分的整合,并投入到可靠性分析的方法之中,这其中最为典型几种方法除了响应面之外,还有神经网络、支持向量机和Kriging等几种模型。神经网络法最开始是由Jan Deng等人推出的[20],解决了一次二阶矩、二次二阶距方法等方法无法解决的工程问题,也就是当功能函数为隐式的情况下的应用。程进等人提出了一种改进的响应面法[21],该方法很好地融合了其他可靠性方法的优势,对计算分析的效率和准确度有了较大的提升。孟广伟等人推出了一种亏优化的神经网络法[22],该法能有效攻克功能函数非线性程度高同时为隐式的可靠性问题;谭晓慧等人推出了一种可优化的响应面法,同时考察了该方法在可靠性分析中的运用状况[23],这种方法在所有求解中仅仅只要拟合一个响应面,此种方法在分析的准确性和效率上比一般的方法都要好;史妍妍等人也是探索得出一种根据响应面法的隐式极限状态函数可靠性灵敏度分析方法[24],同时根据实际应用证明了此法的普适性;陈学前等人在对响应面法在结构参数灵敏度及可靠性分析中的实际应用研究[25],整个过程运用的是Box.Behnken矩阵设计方法,按照分析出的设计点的响应,采用最小二乘法构建响应面函数;阎宏生等人对根据神经网络响应面的结构可靠性分析方法进行了进行了一定探索,针对二次多项式响应面法暴露出的不足,构建了神经网络响应面,从而推出一种根据神经网络响应面的结构可靠性分析方法[26];吕震宙等人在研究中,发现了一种基于神经网络的可靠性分析的新方法[27],同现有方法进行比较,所提方法选择近似极限状态方程而不是近似极限状态函数的策略,并通过筛选训练样本实现这一策略,从而提高可靠性分析的精度;B.Echard等人提出了一种将重要抽样与Kriging模型相结合的可靠性分析方法[28],他们是把Kriging模型和主动学习函数进行了充分结合,在这个过程中采用重要抽样法来降低计算量;Tong Cao等人便是据此,研究得出一种新的混合型可靠性分析方法[29];黄晓旭等人亦是如此,将两种关联的方法结合[31],同时进行综合比较,分析并证实了上述方法存在的优势和弊端。

  随着时间的推移,上述各类方法也是在不断的发展,现如今,基于已确定性优化设计的方法,已推出了新的方法,即结构可靠性优化设计方法。这种方法是把可靠性分析方法与数学上的优化方法进行了有机的结合,其目标函数或者是约束函数的考查主要是通过以可靠度指标来实现,整个过程中,并运用最佳的求解方法去求解,最后的结果是结构在符合可靠性要求下的最优化的设计变量值的数值方法。在1960-1980间,一些比较古老的结构确定性优化设计方法的发展十分迅速,但由于传统方法没有考虑结构设计中不确定性因素的影响,因此,寻求更加合理的结构可靠性优化设计方法成为目前学者们积极研究的重点和热点之一。在结构设计的过程中,会存在很多的不确定性因素,例如制造过程中的不确定性、村料属性的不确定性、操作环境的不确定性等,因为在整个设计中将与之对应的不确定性因素进行了思考,所以使得该方法能够较好地反映机械成品的费用和可靠性两者之间是否均衡的关系。不难得出,在本次提出的设计方法中,有着2个循环流程,它们分别是外层和内层循环,前者主要把确定性优化循环这一操作进行优化以便得到最优设计变量,然而后者则是真正的可靠性分析循环。本文提供的方法可以很好地解决不确定性的不足,因其有着非常高的效率以及精确度,所以在实际应用中得到了很好的应用。另一方面,为了保证当前拟需要的变量能达到可靠度指标提出的标准,通常情况下,能够使用的方法有:性能测度法(PMA)、可靠度指标法(RIA)、改进均值法(AMV)和混合均值法(HMV)四种等,这四种方法利弊不一,但是都是非常传统的计算方法。

  1.3本文研究的主要内容和方法

  针对实际工程中的隐式功能函数的小失效概率问题,本文基于提出的一种主动学习的Kriging代理模型与Monte Carlo模拟方法相结合的可靠度分析方法—AK-MCS。针对工程实际中大量存在的小失效概率问题,AK-MCS方法有求解失效概率和主动学习的Kriging模型代替真实功能函数的优势。该方法首先釆用Kriging模型代替真实功能函数,通过主动学习方法逐步扩充实验设计点,逐步改善Kriging模型的精度;然后利用Monte Carlo模拟方法的基本思路,通过引入合理的随机变量计算小失效概率。结果表明,AK-MCS方法在保证结果精度的同时减少了功能函数的评估次数,对于工程实际中具有隐式功能函数的小失效概率计算问题具有较强的应用前景。

  1.4本章小结

  第1章主要介绍了该课题的研究背景,国内外的发展及现状,以及第2章主要介绍桁架结构的有限元原理和MATLAB软件的基本操作和桁架的程序编写,通过三个实际的结构对桁架进行分析,求解位移,同时用ANSYS软件对程序的结果进行验证。第3章主要介绍可靠性的基本理论和几种可靠性的分析方法,针对第2章的三个实例进行蒙特卡洛法和AK-MCS法的可靠性分析,并比较两者的优劣。

  第2章桁架结构的有限元程序

  2.1桁架及有限元法简介

  对于实际的复杂工程结构,有限元分析是近代发展起来求解复杂结构位移的有效方法。其基本思想是将连续的求解区域离散为一组有限个单元、且按一定方式互相排列制约的单元组合体,由于单元能按照不同的方式组合在一起,单元本身也有不同的形状,可以把几何形状复杂的求解域模型化,在每一个单元内假设近似函数来分片求解。随着单元数的增加,也就是单元尺寸的缩小,解的近似程度将会不断收敛于精确解,直至满足精度要求。有限元方程以节点位移作为未知量,求解方程组从而得到总体连续的节点位移,最终求出单元应力等。以便继续进行后续的结构分析。

  桁架(truss)是一种由杆件彼此在两端用铰链连接而成,可以绕连接点转动的结构。它简单高效,应用广泛。桁架由是直杆组成的一般具有三角形单元的平面或者空间结构,桁架杆件主要承受的是轴向拉力或压力,从而能充分利用材料的强度,在跨度较大时比实腹梁节省材料,减轻自重同时增大刚度。各杆的受力为拉或压,因此通过合理的布置,可以适应结构的内部剪力分布和弯矩分布,整个结构对支座不产生水平推力。而且其受拉压截面两端,内里臂增大了所以在相同的材料用量下,桁架有更加优良的抗弯强度。通过合理的布置,剪力可以逐步的被传递给支座。综合看来,这样的桁架结构有良好的抗弯和抗剪能力。更重要的是,相对于梁而言,它把横弯条件下的实腹梁内部复杂的应力状态经过转化,成为杆件内部简单的拉压应力状态,更为直观。我国桁架结构发展迅速,在许多个方面都有应用,如大型建筑、机械、船舶工程等,在建筑领域得也到了许多青睐。结构力学中有节点法和截面法,在有限元问题里,把它按理想的铰接模型计算,桁架杆件内的位移分布函数可以得到精确的结果,所以不必再细分单元即可,直接以每一杆件为一单元。

  本文采用MATLAB作为编程平台,对平面桁架和空间桁架进行静力学分析。首先对整体结构进行分析,确定结构的输入参数,如横截面积、弹性模量、载荷等,然后对节点和单元进行编号,根据有限元理论方法写出单元刚度矩阵,依据每根杆单元的角度写出坐标转换矩阵,然后可得整体坐标系下的刚度矩阵。再进行刚度矩阵组装,得到的整体刚度矩阵可写出整体坐标方程,根据实际条件确定边界条件,代入方程中求解得出未知点的节点位移。

  具体的程序如后文的例子所示。

  2.2 MATLAB基本语言

  在我们国家,现如今应用MATLAB的一般人群通常是学生和一些科研院所,而商业化应用的并不多见。为何MATLAB深受人们的喜爱,经过分析,有以下点原因:

  1)具有非常高效的数值计算性能。单从这点就足矣,现如今其他编程语言和相似的数学平台都是不能与之匹敌的;

  2)有着非常齐全的计算结果和编程可视化功能。同理,以此也是很有竞争力的,难以被其他软件超越;:直接调用MATLAB的Figure潜入到WinForm中去,这和MATLAB的结果优秀的可视化功能密不可分;

  3)有着极其齐全的应用工具箱及Help系统,现如今它的工具箱数量已非常之多,可以说涵盖了包括数学、统计等所有领域。如果要说MATLAB的功能强,那其根本的原因还是其工具箱数量多,这其中的任何一个工具箱都有与之不同领域内比较常见的某些算法和应对方法。这可以使得时间得到节省,正因如此,对于科研工作者来讲,运用它可以很迅速的将自己的想法进行检验,它对实现算法和进行测试来说是十分关键的。

  4)它具有非常友好优化的编程设计条件,同时含有比较相近的数学计算公式的自然化m口令。学习时极易接受;MATLAB就是它的编程设计条件之一,自带的m口令极易掌握,特别是一些有编程设计的人员,更容易学会;

  工作界面介绍

  Current Folder:指的是当前存储的工作文件路径。同时也是当前MATLAB工作文件夹的路径,该路径通常情况下,一旦启动更改完成后,就不能轻易更改,但在特别的条件下,还是能够用在文件件位置的查找上;此时,便会加载用户一般使用的MATLAB目录,同时能够下拉,迅速的进行切换不同的MATLAB工作路径,实现对工作的便利;

  Current Folder:同上,它会将整个资源都在其中显示,如果用户要对其进行打开,用户只需要对相应的m文件进行双击操作;

  WorkSpace,表示的是工作变量空间,通常它所显示得是当前MATLAB中存在的变量的值,涵盖有变量的名称、值等内容。若这些内容是数组,那么它会以最大和最小值的方式进行显示,它的作用类似于调试,相当于VS中加断点后的局部变量的值,但这里显示更加直观,非常重要;

  Command Window主窗口,也是输入命令的地方,MATLAB中极为重要的一部分,做一些简单的测试,学习命令的时候都可以在这里面进行,也可以输入函数所需变量,但他的输入具有限制,修改变量等十分不方便,主要还是简单的操作,更复杂的仍然是在m文件中完成;

  Command History,历史命令窗口。这里可以选择要关闭,但在出错时会很方便,即使用者在命令的窗口将自己的命令进行输入,那么均会得到显示,此时可把键盘上的左右及上下键进行结合使用,上下键表示可以转至前面的命令中,这个时候如果进行测试,那么时间将会缩短,也会更易操作;

  单个m文件或者函数,编写脚本的位置,方便程序的引用、保存及调用等;

  常见的命令

  1.clear:清除内存变量和函数,包括矩阵等,也就是把工作区间WorkSpace的变量清空;

  2.clc:清楚当前MATLAB命令窗口的内容;重新开始,使界面更简洁一点,注意clc是不清除变量的。一般在编写m文件的时候,不是函数的话,前面一般都要加上clear;clc;目的就是在m文件运行的时候,把内存和屏幕都清空,以免同名称的变量影响以及使得屏幕容易观察;

  3.rand:随机数生成器,可以直接使用生成任意维度的矩阵,例如rand(2,3)产生一个0~1的2×3矩阵,其他的随机数生成还有normrnd,uniform等等不同的随机数分布规律,可以自行考究;

  4.zeros:创建一个都为0的矩阵;ones:创建1个都为1的矩阵。参数可以是多维的,例如zeros(2,3)创建一个2×3的0矩阵;

  5.size:可以计算矩阵的大小,同样相关的还有length,size可以计算不同的维上的大小,例如size(A,1)计算矩阵A的行数,size(A,2)计算矩阵的列数;

  6.help:在想要知道某个函数的相关说明的时候,可以使用help函数名来获取,当然也可以打开帮助文档。不过这种方法是比较快的。

  6.plot:绘制图像的操作,具体格式可以参考MATLAB文档帮助。plot(x,y)绘制以x为横坐标,y为纵坐标的二维图,另外还有三维图绘制,使用它的操作非常频繁,做科研的时候,经常要看趋势,绘图在MATLAB是非常常见的一件事情,subplot函数是将一个figure分割为多个块来操作。

  7.figure,hold on:在使用plot绘图的时候,基本的默认都是在figure上面,figure可以新建一个空白图像,如有需要,还可以新建多个空白图像,同时还可以使用hold on命令在同一个figure上面绘制多条曲线。

  基本操作

  MATLAB中变量名区分大小写,命名规则没有权威的规定,基本都是小写开头,以罗马数字和英文字母组成,注意,如果一行语句结束,后面加分号则只显示在工作区域不加分号,就会直接在Command Window显示变量的值,在命令行窗口也是类似的道理,MATLAB中,注释的符号是%,即百分号,但注意该软件中的所以语句都要用英文输入法完成,否则会变红无法被识别。变量无需定义,你可以给它一个空值,在使用的时候再赋其他值。它的长度不是固定的,你可以按你需要和需求进行增加或删减。矩阵操作很灵活简便,大部分的时候都要考虑不使用循环而是矩阵能否完成,因为在该软件下,一旦循环的层数较多,耗费的时间将会很长,因此不建议这样做。矩阵操作一般都可以直接一句话完成,然而免不了还是要用循环的。MATLAB也有for,while等循环语法。每一层的for或者while都要对应end才行,和其他的编程软件一样。编写函数使用的关键字是function,以end结尾,输入和输出参数都可以是多个,例如function[r1,r2,r3]=testfun(p1,p2,p3),但改函数在命令窗口是不能直接运行的,需写一个m文件调用所写的函数才能运行。

  2.3 ANSYS软件验证结果

  ANSYS是当今国际范围内发展速度最快的计算机辅助工程(CAE)软件,是美国ANSYS公司研发的大型通用有限元分析(FEA)软件,能与多数计算机辅助设计(CAD,computer Aided design)软件相互接口,实现数据的共享和交换,方便各种学术交流和验证求解等等。如Creo,、AutoCAD NASTRAN、Algor和I-DEAS等。是集合了融声场、结构、磁场、流体、电场、分析为一体的大型通用有限元分析软件。在铁道、能源、石油化工、航空航天、机械制造、汽车交通、土木工程、核工业、电子、造船、生物医学、地矿、水利、日用家电、轻工以及新兴的领域等领域有着广泛的应用。ANSYS的作用非常强大,操作也是非常的简单,目前是世界范围内比较受欢迎的一种有限元分析软件,在这几年的FEA环比中,都是佼佼者。现如今,国内100多所高等学府都运用了这种软件,正是如此,这种软件已经成为学术研究中必不可少的一种软件。该款软件含有一个功能较多的有限元法计算机设计程序,这种程序在处理求解结构、流体等相关领域的问题时,可以很好的求解出可靠的结果。在处理静力分析时,通常可以求解外载荷引起的位移、应力和力。这种分析方法在处理求解阻尼对结构的影响并不显著和惯性的难题时行之有效。该程序下的静力分析无论是对线性分析还是非线性分析都是能够实现的,如塑性、膨胀、蠕变、大变形、大应变及接触分析。

  该款软件系统通常含有3个部分,他们分别是前处理模块,分析计算模块和后处理模块。第一个模块可为使用者提供一个的实体建模及网格分割的工具,使用者可以很简单的使用它进行建模;第二个模块涵盖了结构分析、流体动力学分析、电磁场分析、压电分析、声场分析以及多物理场的耦合分析,能够对多种物理介质发生的作用进行模拟,它的灵敏度和优化性都是较好的;后处理模块可将计算结果以彩色等值线显示、立体切片显示、矢量显示、梯度显示、粒子流迹显示、透明及半透明显示(可看到结构内部)等多种图形方式显示清楚,甚至也可将计算结果以图表、曲线形式显示以及输出,以便继续进行更好的分析。同时,该软件也能提供多种单元类型和不同的属性定义,以便做到满足模拟材料和各种工程上结构,和实际结合很优良。因此它可应用于以下广泛的工业领域:汽车工业、电子产品、生物医学、桥梁、建筑、航空航天、重型机械、微机电系统等。如今,越来越多的的领域都可以用它来研究,呈现着很好的势头。

  利用ANSYS软件再次对MATLAB求解程序的正确性进行验证。

  2.4三个桁架结构的有限元程序

  2.4.1平面4杆桁架

  如图2-1,各杆的弹性模量和横截面积相同:均为E=29.5×N/,A=100,求解该结构的节点位移。

  图2-1 4杆桁架结构

  节点及单元编号按图2-1上所示,编写MATLAB有限元程序:

  function disp=FEM4bar(x)

  %x=[A E F1 F2]=[100 295000 20000-25000]

  A=x(:,1);Q=[x(3)0 x(4)];%x为随机变量

  A_t=[A A A A];%截面面积

  E_t=x(1,2)*ones(1,4);%弹性模量

  row_t=[1 2;2 3;1 3;4 3];%单元的节点号

  Node=[0 0;400 0;400 300;0 300];%节点的坐标

  for j=1:size(row_t,1)

  Node1=Node(row_t(j,1),:);

  Node2=Node(row_t(j,2),:);

  plot([Node1(1)Node2(1)],[Node1(2)Node2(2)])

  hold on

  end

  l_t=zeros(1,size(row_t,1));%单元长度

  sita_t=[0 90 45 0];%单元的角度--求转换矩阵

  for j=1:size(row_t,1)

  Node1=Node(row_t(j,1),:);

  Node2=Node(row_t(j,2),:);

  l_t(j)=norm(Node2-Node1);%单元长度

  end

  k_t=zeros(8,8);

  for ne=1:4

  A=A_t(ne);

  sita=sita_t(ne);

  l=l_t(ne);

  E=E_t(ne);

  row=row_t(ne,:);

  line=row;

  tt=[cosd(sita)^2 sind(sita)*cosd(sita);sind(sita)*cosd(sita)sind(sita)^2];%转换矩阵

  ke=[tt-tt;-tt tt]*A*E/l;%单元刚度矩阵

  k2=zeros(8,8);

  for i=1:2%刚度矩阵组装

  for j=1:2

  k2(row(i)*2-1:row(i)*2,line(j)*2-1:line(j)*2)=ke(i*2-1:i*2,j*2-1:j*2);

  end

  end

  k_t=k_t+k2;

  end

  K_t=k_t;%备份刚度矩阵求支反力

  boundary=[1 2 4 7 8];%边界条件

  k_t(boundary,:)=[];%划掉行

  k_t(:,boundary)=[];%划掉列

  K=k_t;

  invK=Keye(3);

  X=invK*Q’;;%求逆

  disp=X;%节点位移

  end

  平面四杆桁架的MATLAB程序运行结果如图2-2和表2-1所示:

  图2-2 MATLAB程序绘制结构

  表2-1平面4杆的节点位移

  节点x轴(mm)y轴(mm)

  1 0 0

  2 0.2712 0

  3 0.0598-0.2094

  4 0 0

  利用ANSYS对结构再次进行计算,以验证该程序正确与否。

  ANSYS软件内大致步骤:

  设置一个新的静力学分析,选择单元为LINK180。

  设置单元属性,包括横截面积。密度和泊松比可暂不输入,此处对机构本身的重力不作考虑。

  建立有限元模型。选取节点,按顺序输入节点的坐标,再根据节点连接单元完成有限元模型。

  划分网格。由于桁架结构不需要将每个杆单元另外划分成有限个单元,否则无法求解,因此将单元划分设置成1。

  添加约束并加载。对所有节点z向位移进行约束,限制节点1和4的x、y向位移,以及2节点的y向位移。

  求解并查看求解结果。选择绘制变形图,结构的受力情况和变形情况十分清楚。

  图2-3 4杆桁架ANSYS求解结果

  由图2-3可见,该结构的最大位移出现在2节点处,为0.271186mm,MATLAB结果0.2712mm与ANSYS几乎一致。最大位移仅相差0.000014mm。说明MATLAB程序求解的该4杆桁架结构是准确的。

  2.4.2平面61杆桁架

  如图,平面61杆桁架其弹性模量、均为E=2.1×Pa,所有单元横截面积均为1×,在10、12、14、16、18五个节点分别施加一个竖直向上1000N的载荷,求解节点位移。

  图2-4平面61杆桁架结构

  在MATLAB软件编写有限元程序:

  function disp=FEM61bar(x)%x=[A E l F]=[1e-4 2.1e11 1 1000]

  l=x(:,3);P=x(:,4);%x为随机变量

  A_t=x(1,1).*ones(1,61);%截面面积m^2

  E_t=ones(1,61)*x(:,2);%弹性模量Pa

  row_t=[1 2;1 4;2 3;1 3;2 4;3 4;3 6;4 5;3 5;4 6;5 6;5 8;6 7;5 7;6 8;7 8;7 10;8 9;7 9;8 10;9 10;9 12;10 11;9 11;10 12;11 12;11 14;12 13;11 13;12 14;13 14;13 16;14 15;13 15;14 16;15 16;15 18;16 17;15 17;16 18;17 18;17 20;18 19;17 19;18 20;19 20;19 22;20 21;19 21;20 22;21 22;21 24;22 23;21 23;22 24;23 24;23 26;24 25;23 25;24 26;25 26];%单元的节点号

  Node=[0 0;0 1;1 0;1 1;2 0;2 1;3 0;3 1;4 0;4 1;5 0;5 1;6 0;6 1;7 0;7 1;8 0;8 1;9 0;9 1;10 0;10 1;11 0;11 1;12 0;12 1]*l;%节点的坐标

  for j=1:size(row_t,1)

  Node1=Node(row_t(j,1),:);

  Node2=Node(row_t(j,2),:);

  plot([Node1(1)Node2(1)],[Node1(2)Node2(2)])

  hold on

  end

  l_t=zeros(1,size(row_t,1));%单元长度

  s=[90 45-45 0 0];

  sita_t=[s s s s s s s s s s s s 90];%单元的角度--求转换矩阵

  for j=1:size(row_t,1)

  Node1=Node(row_t(j,1),:);

  Node2=Node(row_t(j,2),:);

  l_t(j)=norm(Node2-Node1);%单元长度

  end

  k_t=zeros(52,52);

  for ne=1:61

  A=A_t(ne);

  sita=sita_t(ne);

  l=l_t(ne);

  E=E_t(ne);

  row=row_t(ne,:);

  line=row;

  tt=[cosd(sita)^2 sind(sita)*cosd(sita);sind(sita)*cosd(sita)sind(sita)^2];%转换矩阵

  ke=[tt-tt;-tt tt]*A*E/l;%单元刚度矩阵

  k2=zeros(52,52);

  for i=1:2%刚度矩阵组装

  for j=1:2

  k2(row(i)*2-1:row(i)*2,line(j)*2-1:line(j)*2)=ke(i*2-1:i*2,j*2-1:j*2);

  end

  end

  k_t=k_t+k2;

  end

  boundary=[1 2 49 50];%边界条件

  k_t(boundary,:)=[];%划掉行

  k_t(:,boundary)=[];%划掉列

  K=k_t;

  invK=Keye(48);

  Q=zeros(1,48);

  force=[18 22 26 30 34];

  Q(:,force)=P;%载荷向量

  X=invK*Q';%求逆

  disp=(reshape(X,2,[]))';%节点位移

  [disp,~]=max(disp);

  disp=max(disp);%输出最大节点位移

  end

  MATLAB软件中上述程序运行结果如图2-5和表2-2所示:

  图2-5 MATLAB程序结构绘制

  表2-2 MATLAB程序求解的61杆桁架节点位移

  节点编号横坐标(m)纵坐标(m)

  1 0 0

  2-0.00181325 1.80E-05

  3 0.000271351 0.002153528

  4-0.001795246 0.002081041

  5 0.000434208 0.004391268

  6-0.001647641 0.004328465

  7 0.000477146 0.006531561

  8-0.00138186 0.006465745

  9 0.000398896 0.008318756

  10-0.000999171 0.008272341

  11 0.00022314 0.009493148

  12-0.000523512 0.0094447

  13 3.45E-17 0.009900029

  14-4.80E-18 0.009851815

  15-0.00022314 0.009493148

  16 0.000523512 0.0094447

  17-0.000398896 0.008318756

  18 0.000999171 0.008272341

  19-0.000477146 0.006531561

  20 0.00138186 0.006465745

  21-0.000434208 0.004391268

  22 0.001647641 0.004328465

  23-0.000271351 0.002153528

  24 0.001795246 0.002081041

  25 0.00181325 1.80E-05

  26 0 0

  从表2-2可以看出,该桁架的变形呈对称状态,最大节点位移为0.0099m,出现在节点13和节点14处。

  利用ANSYS软件再次对结构进行求解分析,以验证所MATLAB软件求得的结果。

  ANSYS软件内的大致步骤:

  设置一个新的静力学分析,选择单元为LINK180。

  设置单元属性,包括横截面积、弹性模量。密度和泊松比可暂不输入,此处对机构本身的重力不作考虑。

  建立三维模型。选取节点,按顺序输入节点的坐标,再根据节点连接单元完成有限元模型。这种方法简洁直观,但该结构是61杆26节点桁架结构,操作起来较为耗时,也容易出错。在这里,可以考虑使用命令流建立模型,编辑好各节点及单元的程序,输入到命令流窗口即可。

  划分网格,建立有限元模型。由于桁架结构不需要将每个杆单元另外划分成有限个单元,否则无法求解,因此将单元划分设置成1。

  添加约束并加载。对所有节点z向位移进行约束,限制节点1和节点25的x、y向位移。

  求解并查看求解结果。选择绘制变形图,结构的受力情况和变形情况一目了然。

  图2-6 61杆平面桁架AYSYS结果

  从图2-6可见,最大位移是0.0099m,出现在中心的13和14节点处,且结构变形情况呈对称分布。ANSYS软件的结果与MATLAB程序结果一致,说明MATLAB求解程序是正确的。

  2.4.3空间132杆桁架

  图2-7 132杆空间桁架

  如图2-7,空间132杆桁架结构弹性模量为E=205.8MPa,横截面积为A=100,其六分之一的节点编号以及坐标图2-8和表2-3所示,在竖直位移大于零的37个节点处分别施加16kN的力,约束竖直位移为零的节点。求解节点位移。

  图2-8 132杆桁架1/6结构节点编号图

  节点编号x(mm)y(mm)z(mm)

  1 0 0 4618.8

  2 2390.87 0 4304.04

  3 1195.435 2070.554 4304.04

  4 4618.8 0 3381.2

  5 4000 2309.4 3381.2

  6 2309.4 3999.998 3381.2

  7 6531.97 0 1913.17

  8 6138.05 2234.07 1913.17

  9 5003.78 4198.67 1913.17

  10 3265.99 5656.85 1913.17

  11 8000 0 0

  12 7727.41 2070.55 0

  13 6928.2 4000 0

  14 5656.85 5656.85 0

  15 4000 6928.2 0

  表2-3 132杆桁架1/6结构节点坐标

  根据已知节点坐标求出所有的节点坐标,进行有限元程序编写。

  首先据已知节点坐标求出所有的节点坐标,进行MATLAB有限元程序编写。

  function[disp]=FEM132bar(x)

  m=132;%单元数

  n=61;%节点数

  X=zeros(1,112);

  for j=1:size(x,2)

  X(:,(3*(j-1)+1))=x(:,j);%扩展x、y方向的力

  end

  A_t=ones(1,m)*100;%横截面积100mm^2

  E_t=ones(1,m)*x(1);%弹性模量均值205.8Mpa

  row_t=[1,2;1,12;1,22;1,32;1,42;1,52;2,12;12,22;22,32;32,42;42,52;2,52;2,3;2,4;4,12;12,13;12,14;14,22;22,23;22,24;24,32;32,33;32,34;34,42;42,43;42,44;44,52;52,53;52,54;2,54;3,4;4,13;13,14;14,23;23,24;24,33;33,34;34,43;43,44;44,53;53,54;3,54;3,5;3,6;4,6;4,7;7,13;13,15;13,16;14,16;14,17;17,23;23,25;23,26;24,26;24,27;27,33;33,35;33,36;34,36;34,37;37,43;43,45;43,46;44,46;44,47;47,53;53,55;53,56;54,56;54,57;3,57;5,6;6,7;7,15;15,16;16,17;17,25;25,26;26,27;27,35;35,36;36,37;37,45;45,46;46,47;47,55;55,56;56,57;5,57;5,8;5,9;6,9;6,10;7,10;7,11;11,15;15,18;15,19;16,19;16,20;17,20;17,21;21,25;25,28;25,29;26,29;26,30;27,30;27,31;31,35;35,38;35,39;36,39;36,40;37,40;37,41;41,45;45,48;45,49;46,49;46,50;47,50;47,51;51,55;55,58;55,59;56,59;56,60;57,60;57,61;5,61];

  Node=[0,0,4618.80000000000;2390.87000000000,0,4304.04000000000;4618.80000000000,0,3381.20000000000;4000,2309.40000000000,3381.20000000000;6531.97000000000,0,1913.17000000000;6138.05000000000,2234.07000000000,1913.17000000000;5003.78000000000,4198.67000000000,1913.17000000000;8000,0,0;7727.41000000000,2070.55000000000,0;6928.20000000000,4000,0;5656.85000000000,5656.85000000000,0;1195.43500000000,2070.55415714610,4304.04000000000;2309.40000000000,3999.99813499957,3381.20000000000;0.000932500217913912,4618.80161513775,3381.20000000000;3265.98500000000,5656.85195675784,1913.17000000000;1134.26362616730,6432.74222969907,1913.17000000000;-1134.26488210761,6432.73559494850,1913.17000000000;4000.00000000000,6928.20323027551,0;2070.55610019413,7727.40836545791,0;-0.00161513775356070,7999.99720249935,0;-2070.55080539800,7727.40080539800,0;-1195.43500000000,2070.55415714610,4304.04000000000;-2309.40000000000,3999.99813499957,3381.20000000000;-3999.99906749978,2309.40161513776,3381.20000000000;-3265.98500000000,5656.85195675784,1913.17000000000;-5003.78637383270,4198.67222969908,1913.17000000000;-6138.04488210761,2234.06559494850,1913.17000000000;-4000.00000000000,6928.20323027551,0;-5656.85389980587,5656.85836545791,0;-6928.20161513775,3999.99720249935,0;-7727.40080539800,2070.55080539800,0;-2390.87000000000,2.92917068638737e-13,4304.04000000000;-4618.80000000000,5.65871568353193e-13,3381.20000000000;-4000.00000000000,-2309.40000000000,3381.20000000000;-6531.97000000000,8.00263295300946e-13,1913.17000000000;-6138.05000000000,-2234.07000000000,1913.17000000000;-5003.78000000000,-4198.67000000000,1913.17000000000;-8000,9.80118763926896e-13,0;-7727.41000000000,-2070.55000000000,0;-6928.20000000000,-4000.00000000000,0;-5656.85000000000,-5656.85000000000,0;-1195.43500000000,-2070.55415714610,4304.04000000000;-2309.40000000000,-3999.99813499956,3381.20000000000;-0.000932500219732901,-4618.80161513775,3381.20000000000;-3265.98500000000,-5656.85195675784,1913.17000000000;-1134.26362616730,-6432.74222969907,1913.17000000000;1134.26488210761,-6432.73559494850,1913.17000000000;-4000.00000000000,-6928.20323027551,0;-2070.55610019413,-7727.40836545791,0;0.00161513775037747,-7999.99720249935,0;2070.55080539800,-7727.40080539800,0;1195.43500000000,-2070.55415714610,4304.04000000000;2309.40000000000,-3999.99813499957,3381.20000000000;3999.99906749978,-2309.40161513775,3381.20000000000;3265.98500000000,-5656.85195675784,1913.17000000000;5003.78637383270,-4198.67222969907,1913.17000000000;6138.04488210761,-2234.06559494850,1913.17000000000;4000.00000000000,-6928.20323027551,0;5656.85389980587,-5656.85836545791,0;6928.20161513776,-3999.99720249935,0;7727.40080539800,-2070.55080539800,0];

  %%绘图

  %plot

  %for ne=1:m

  %row=row_t(ne,:);

  %Node1=Node(row(1),:);

  %Node2=Node(row(2),:);

  %Node_corrd=[Node1;Node2];

  %plot3(Node_corrd(:,1),Node_corrd(:,2),Node_corrd(:,3))

  %%view(2);%2D

  %hold on

  %mid_Node=(Node1+Node2)/2;

  %c=num2str(ne);

  %c=['(',c,')'];

  %text(mid_Node(1),mid_Node(2),mid_Node(3),c)

  %end

  %for n_node=1:n

  %c=num2str(n_node);

  %c=['',c];

  %text(Node(n_node,1),Node(n_node,2),Node(n_node,3),c)

  %end

  k_t=zeros(3*n,3*n);

  for ne=1:m

  row=row_t(ne,:);

  Node1=Node(row(1),:);

  Node2=Node(row(2),:);

  A=A_t(ne);

  l=norm(Node2-Node1);

  E=E_t(ne);

  line=row;

  cos_l=(Node1(1)-Node2(1))/norm(Node2-Node1);

  cos_m=(Node1(2)-Node2(2))/norm(Node2-Node1);

  cos_n=(Node1(3)-Node2(3))/norm(Node2-Node1);

  %tt=[cos_l*cos_l cos_l*cos_m cos_l*cos_n-cos_l*cos_l-cos_l*cos_m-cos_l*cos_n;

  %cos_l*cos_m cos_m*cos_m cos_m*cos_n-cos_l*cos_m-cos_m*cos_m-cos_m*cos_n;

  %cos_l*cos_n cos_m*cos_n cos_n*cos_n-cos_l*cos_n-cos_m*cos_n-cos_n*cos_n;

  %-cos_l*cos_l-cos_l*cos_m-cos_l*cos_n cos_l*cos_l cos_l*cos_m cos_l*cos_n;

  %-cos_l*cos_m-cos_m*cos_m-cos_m*cos_n cos_l*cos_m cos_m*cos_m cos_m*cos_n;

  %-cos_l*cos_n-cos_m*cos_n-cos_n*cos_n cos_l*cos_n cos_m*cos_n cos_n*cos_n

  %];

  %ke=A*E/l*tt;

  %%

  cos2=(cos_l^2+cos_m^2).^0.5+eps;

  tt=[cos_l cos_m cos_n;

  -cos_m/cos2 cos_l/cos2 0;

  -cos_n*cos_l/cos2-cos_n*cos_m/cos2 cos2];

  T=[tt zeros(3,3);zeros(3,3)tt];

  ke=zeros(6);

  ke(1,1)=1;ke(1,4)=-1;ke(4,1)=-1;ke(4,4)=1;

  ke=A*E/l*T'*ke*T;

  k2=zeros(3*n,3*n);

  for i=1:2

  for j=1:2

  k2(3*(row(i)-1)+1:row(i)*3,3*(line(j)-1)+1:line(j)*3)=ke(3*(i-1)+1:i*3,3*(j-1)+1:j*3);%刚度矩阵组装

  end

  end

  k_t=k_t+k2;

  end

  %%加载求位移

  Load=zeros(n,3);

  Load([1:7,12:17,22:27,32:37,42:47,52:57],1:3)=(reshape(X(2:112),[3,37]))';%分布荷载

  Re=ones(n,3);

  Re([1:7,12:17,22:27,32:37,42:47,52:57],1:3)=zeros(37,3);

  Re=reshape(Re',[1,3*n]);

  boundary=find(Re==1);

  Q=reshape(Load',[1,3*n]);%返回1*(3*n)矩阵

  Q(boundary)=[];

  K=k_t;

  K(boundary,:)=[];%划掉行

  K(:,boundary)=[];%划掉列

  X=KQ';

  glob_disp=zeros(3*n,1);

  glob_disp(boundary)=0;

  glob_disp(find(Re~=1))=X;%返回Re不等于1,得到所有节点位移

  Node_disp=reshape(glob_disp,3,n);

  disp=max(abs(Node_disp(3,:)));%返回节点最大位移

  %%绘制变形图

  %Node_def=Node+50*Node_disp';

  %for ne=1:m

  %row=row_t(ne,:);

  %Node1=Node(row(1),:);

  %Node2=Node(row(2),:);

  %Node_corrd=[Node1;Node2];

  %plot3(Node_corrd(:,1),Node_corrd(:,2),Node_corrd(:,3))

  %%view(2);%2D

  %hold on

  %end

  %for ne=1:m

  %row=row_t(ne,:);

  %Node1=Node_def(row(1),:);

  %Node2=Node_def(row(2),:);

  %Node_corrd=[Node1;Node2];

  %plot3(Node_corrd(:,1),Node_corrd(:,2),Node_corrd(:,3),':')

  %end

  end

  在MATLAB软件中运行上述程序,得到的结果如图2-9和图2-10,每个节点的位移见附录1。

  图2-9 132杆桁架节点及单元编号

  图2-10位移放大50倍结构变形图

  根据MATLAB的运行结果,其节点位移见附录1,可知,最大位移是中心处的1号节点,z向位移是17.32245mm。

  利用ANSYS软件再次分析结构,验证MATLAB程序的结果。

  ANSYS软件内大致步骤:

  设置一个新的静力学分析,选择单元为LINK180。

  设置单元属性,包括横截面积、弹性模量。密度和泊松比可暂不输入,系统自动设置为零,此处对机构本身的重力不作考虑。

  建立三维模型。选取节点,按顺序输入节点的坐标,再根据节点连接单元完成有限元模型。这种方法简洁直观,但该结构是132杆61节点空间桁架结构,操作起来较为耗时,也容易出错。在这里,可以考虑使用命令流建立模型,编辑好各节点及单元的程序,输入到命令流窗口即可。注意检查结构是否正确,然后进行下一步。

  划分网格,建立有限元模型。由于桁架结构不需要将每个杆单元另外划分成有限个单元,否则无法求解,因此将单元划分设置成1。

  添加约束并加载。对所有在平面的节点约束x、y、z向的位移。

  求解并查添加其余37个节点的z向分布载荷。

  查看求解结果。选择绘制变形图,结构的受力情况和变形情况一目了然。

  由图2-11可见,ANSYS软件的结果与MATLAB结果几乎一致,最大节点位移出现在最中心处的顶点,与MATLAB结果相差仅0.0482mm。

  图2-11 ANSYS求解变形图

  2.4本章小结

  本章介绍了相关有限元的理论知识和MATLAB基础的操作知识,通过MATLAB软件,对桁架结构进行了程序编写,绘制结构图形,建立模型并按照有限元的过程求解位移,同时利用ANSYS软件建模,对所编写的三个例子进行结果验证,最终两种方法的结果几乎一致。说明MATLAB程序是正确的。

  第3章桁架结构的可靠性分析

  3.1可靠性基本理论

  工程结构必须满足安全性、适用性和耐久性,安全性是指在正常施工的情况下,结构仍然可以保持所需的整体稳定性,在承受可能出现的各种外界作用下,如狂风,暴雨,轻微地震等,在预计到的这些事件发生后,机构应该维持正常工作状态。结构的可靠性是指在规定的时间和条件的情况下,结构能够完成特定的功能和指标的能力,适用性是指结构在正常使用时,应该具有优良的性能,使结构的变形、振动或裂缝性能等都不超过规定的限度,可以做到应对这些预计情况的发生。耐久性是结构在正常的使用和维护管理下,具有足够的长时间使用的性能。结构的可靠度是在特定的时间和空间,和一定条件作用下,如外力,自身的材料变化,环境影响等,结构可以完成相应功能的概率。换句话说,可靠性分析就是系统的响应量可以满足特定的概率条件,目前最普遍的方法是在可靠度即失效概率方面上进行分析,因而,本文也主要着重从这个方向研究桁架结构的可靠性。

  几个名词:

  (1)随机变量(basic variables)是在可靠性分析中,对系统的响应量有影响的不确定因素,包括材料的性能、载荷的不确定、以及几何结构的不确定性等。对结构进行可靠性设计,则应充分考虑与之有关的基本随机变量参数。在具体应用中,在进行设计前,所涉及到的每一个参数的特定数值是需要调研或者模拟得到的,否则无法进行下一步的设计和讨论。外界荷载变化,温度变化等等未知的条件,没有办法操作“待建”结构的材料强度为一特定值,所以使得每一个随机参数都是随机变量,开发者唯一可以操作便是运用外界的的数据内容来统计出实际的参数,经验足够,将它们总结成随机设计参数的相对标准的统计规则,这样也为后续的研究提供便利。以上所述规则,组成了结构可靠性分析和设计的基本结构和信息,那么该种分析和设计便应运而生。这些随机变量服从一定规律的分布,例如正态分布(高斯分布)、均匀分布等等。遗憾的是,现阶段,研究中是不好获得随机变量比较准确度的统计数据(如平均值、标准差等),因此在想要分析出精确的分布概率也是很难的,另一方面,某些随机变量的相关的参数,因为没有一些主要数据,那么就需要专家们靠经验去猜测,在这些方面仍然是需要不断地研究的。针对小概率事件的结构失效概率的困难,在处理统计相关参数和概型时,也表现出有差异性的灵敏度,所以运用数值积分法去分析计算结构的失效概率,得出的结果也不能符合研究的要求,特别是精确度方面。所需的分布条件,也可以根据实际生产中的经验取得,否则无法进行下一步的失效概率分析。

  (2)响应量(respones variables)r是在结构对于施加某些规定的条件的结构行为特性,是随机变量的函数,响应量可以是应力、位移、寿命等等,变量与响应量之间是遵循自然规律的极为复杂的函数,或者是没有固定关系式的函数,即隐式函数。可靠性分析的最终结果是得出响应量的统计分布规律,从而得到失效概率。

  (3)极限状态函数(limit state function),也称功能函数(performance function),我国的《工程结构可靠度设计统一标准》(GB 50153—92)对极限状态的定义为:某个结构或结构的一部分超过某一特定状态就不能满足设计规定的某一功能要求,此特定状态称为该功能的极限状态,并且把结构的极限状态分为两种,为承载能力极限状态和正常使用极限状态。它的概念是所有结构或结构中的部分高于某种条件不能实现设计需要的功能要求,那么我们把这种条件或是状态称之为该功能的极限状态。达到这种状,也就意味着结构接近于最大承载能力,或者无法再应用在继续承载的变形状态。当结构发生下列情况之一或其中几种同时存在时,即可以认为大于承载能力极限状态:结构或结构构件将会不稳定,无法保持平衡;所有的结构或者一些结构作为剐体也会出现无法平衡的情况;大于结构构件或其连接由于材质强度问题出现破坏,另外,由于经常性的塑性形变最终导致无法再继续承受载荷,此时,结构将会转变为机动体系。正常使用下极限状态对应于结构或结构构件达到正常使用或者耐久性能的某项要求的限值。当结构或构件发生以下其中一种情况时,就可以判定为超过了常规运用的极限状态:将会干扰其常规的使用或外观的变形;此外,也会干扰常规使用或耐受性的部分破损;将会严重干扰正常运转下的振动,也会干扰其他状态的下常规使用。通常条件下,结构设计可以把承载能力和常规使用的极限状态都要进行一定思量。用于描述结构系统是否处于失效的状态,一般来说是阀值与响应值的差,即,对于没有解析表达式的隐式极限状态函数可靠性分析较为困难,需要通过建模和编程建立相关的求解过程,工作量也较大。功能函数具有的特点:为多个随机变量组成的非线性函数;变量并不都服从正态分布或对数分布;分析结构可靠度变化时需要近似简化,即近似概率法。极限状态函数等于零时称为极限状态方程,是判断失效与否的界线。

  (4)失效概率(failure probability),结构系统的失效的可能性称为失效概率:

  (3-1)

  系统的安全概率,即可靠度:

  (3-2)

  那么,失效概率与可靠度之和等于1:

  (3-3)

  失效概率越小,结构的可靠度越高;反之失效概率越小,结构的可靠度越低。习惯上以失效概率来度量结构的可靠度。

  桁架结构本身具有复杂性,结构的材料特性差异、承载以及约束条件等不确定性影响结构的可靠性,对于这些相关的不确定因素,在利用模型分析计算时,设置其为输入的随机变量。在实际工程中,结构的安全程度的指标称为可靠度,而可靠性指标主要有失效概率、安全度等。

  可靠度分析方法:

  中心点法,也称均值一次二阶距法,它的基本思路是利用随机变量的均值(一阶原点距)和标准差(二阶中心距)模型,分析结构的可靠度,并把极限状态功能函数在均值(即中心点处)作泰勒级数展开,使之线性化,然后求解可靠度指标。中心点法计算简便,概念明确,但存在:基本变量的概率分布不是正态分布或者对数正态分布时,结构可靠度的计算结构会与实际情况比较,差别非常大,通常情况下不适用,然而对非线性函数而言,运用平均值处泰勒级数的方式进行展开也不理想,此外进行展开只有线性项,其他项缺失,因此带来的误差非常大;如果同样的方法云运用有差异的功能函数(不同极限状态方程),可靠指标计算值将明显变化,在一定程度上差别还很大。均值一次二阶矩方法因为它操作时非常简单,计算可靠度指标时只需要计算设计变量的一阶矩和二阶矩,所以在实际工程当中,仍然得到了广泛的应用,它作为一种近似方法,这种方法的使用范围是十分狭窄的,当功能函数的非线性程度很明显时,其计算精度以及结果准确性就不一定可信,于是在使用这种方法时,研究者第一步要确保功能函数的非线性程度是非常低的,这样才可以更好地进行后续的计算和分析。根据根本理论分析,优化一次二阶矩法与均值一次二阶矩法的原理是一致的,Hasofer-Lind率先在其科研成果发布,与均值一次二阶矩法类似,研究人员可以将功能函数进行线性化的拓展,才能进行后续的积分和分析,最后得到的结果基本上可以作为原来功能函数的可靠度结果。那么与均值一次二阶矩法不一样的点就是线性化的点不一样,改进一次二阶矩法是在设计点处将功能函数进行线性化,所以在线性化之前,我们需要求解设计点,设计点的求解一般来说是一个最优化求解的过程。优化一次二阶矩方法同样存在不足的地方,通常有以下几个方面:(1)当使用者处理问题对应的功能函数呈现强烈的非线性时,这时的功能函数本身可能对应有多个设计点,此时研究者在运用最优化方法去进行求解设计点时,肯定会产生问题;(2)当不同非线性程度的功能函数拥有相同的设计点时,此时的可靠度分析结果基本是一致的,不难得出,它是错误的;(3)在实际应用中,人们通常获取不到显式的功能函数计算式,此时运用上述方法将毫无意义。

  验算点法,即JC法,与中心点法不同的是,线性化近似选在失效边界上,而不是选在中心点(均值)处,即以极限状态方程上的一点的切平面作为线性近似,提高可靠指标的计算精度。与中心点法相比,验算点法能够考虑非正态分布的随机变量,在计算量增加不太多的条件下,可对可靠指标进行精度较高的近似计算,求出满足极限状态方程的验算点的响应值。当随机变量是正态分布时,且功能函数是线性方程时,验算法和中心点法的计算结果相同。

  响应面法,因为现阶段机械结构的设计是大家熟知的隐式状态,十分繁琐,它所对相应的功能函数通常也是不易获得,此外,通常会包含有限元分析的一些步骤,响应面的基本思路就是先通过一定形式的试验设计点,并且选取多项式响应面的形式,通过迭代更新策略使多项式响应面方法计算出的失效概率能够在概率上收敛于真实功能函数的失效概率。正是如此,为了使设计的效率得以提升,期望采用某种替代模型的方法来计算失效概率,即采用某种形式的显式替代模型来替换隐式的功能函数,此种方法即为现阶段通常使用的响应面法产生的原因。响应面法最开始是在可靠性分析之中出现,那时主要还是以线性多项式和二次多项式形式居多,对于它的可靠度的分析和计算,通常的方式是先采用一定迭代方式,采用逐步的迭代收敛规则,让多项式响应面对原始功能函数的近似程度符合结束条件。此外,响应面法也能够用在解决功能函数为隐式状态的问题,此时需要把它与有限元分析软件进行结合,同时应用于实际,结果是非常受欢迎的,这正是发挥了它的优势。响应面法因为比较简便,同时它的推导流程也是非常方便简单,因此在可靠性研究中应用的较好,此外,因为该方法不因其他功能函数为隐式的状态而变化,它的应用更是得到了普及,与其他方法相比,有着计算准确度大的优势,结合应用后,更是发展为一种比较系统的方法。响应面法中存在的缺点是,不同形式的响应面方法计算得到的可靠度,对结果的影响十分大。对响应面法来说基本形式是多项式,最重要的是一次以及二次多项式。针对已经出现的问题,可根据它的特点进行有选择性的选取响应面近似函数,然而在处理新生问题的求解时,将失去意义。另一方面,如果结构的极限状态方程非线性的情况非常高时,那么二次多项式方式的响应面法则基本难以实现实时显示出极限状态曲面的非线性的状况,由此会带来分析结果存在误差的问题。因为待确定的参数较多,那么要求越来越多的采样点,虽然使用了高阶次阶的多项式模型,运算起来也难合人心意,非常地不便,难以被广泛的运用。通常采用结构分析和可靠度分析迭代的方法,分析使用响应面法从而计算结构可靠度,这样的计算过程就对工程人员地要求比较高,需要自己另行结构有限元分析,然而根据工程工作者自身有限元分析程序时,无论是在分析精度上,还是它的工作效率和运行范畴方面,与已商业化的分析平台相比,是得不到较好的应用的。

  结构体系可靠度

  实际构件存在许多截面,结构由许多的构建组成,属于结构体系,,若需要其中若干个构件均失效时,结构体系才会失效,该结构为超静定结构;当其中任意杆失效时,结构体系也随之失效,该结构为静定结构。根据结构杆件失效与体系失效之间的关系,可以把实际结构体系分为三种类型:串联体系、并联体系和串并联体系。串联体系是任意构件失效,则结构的体系失效,要求所有的构件都不失效才能安全。并联体系是当中有一个或一个以上的构件失效,剩余构件或失效构件仍然可以保持整体机构的功能。串并联机构实际上是超静定结构,最终的失效形式不限于唯一一种,每一种的失效模式都可以用一个并联体系来模拟,并将其又组成串联体系,组成串并联体系。

  3.2可靠性的Monte Carlo法(MCS)

  Monte Carlo法又称为随机抽样法,统计试验法或者概率模拟法,蒙特卡洛方法的优点是,它不需要考虑结构极限状态曲面的复杂性,没有结构可靠度分析中不便,将所求解的问题同一定的概率模型相联系,从而得到需要的结构的响应,目前,计算机十分发达,借助它进行分析可以很有效的计算出失效概率,实际上,这种方法也是在计算机的发展上而发展起来的。蒙特卡洛模型它是以概率和统计理论方法为基础的一种计算方法,是一种随机模拟办法。为纪念这一方法的概率统计具有的特点,故借用赌城蒙特卡洛命名。1940年左右,S.M.乌拉姆和J.冯·诺伊曼两位科学家在研究核武器时,率先推出该概念。虽然1940年以前,这种方法就开始出现,但是没能得到很好发展。上述方法是基于1777年,法国Buffon推出的可以用投针实验的方式进行计算圆周率的思想,它以概率论和数理统计为基本,根据随机选择变量,然后计算得到待求的响应值,统计所有失效点与总样本点的比值得到,因此它需要试验的样本必须保证足够大,才能满足一定的精度,从而具有一定的可信度。从理论上来说,蒙特卡洛方法需要大量的实验。时间到了20世纪初期,通过千百次的实验,运用蒙特卡罗方法可以计算出圆周率值,但是仍然无法与公元5世纪祖冲之得到结果相媲美。所以该方法一致没有得到较好的应用。随着计算机技术的发展,蒙特卡罗方法在最近10年迅猛发展,计算机在现代生活里已经十分普遍了。现阶段的蒙特卡罗方法,根本不需要人工实验,充分发挥电脑的分析能力,就可以达到效果,使得该方法因比较笨拙的实验过程,不再是其绊脚石。这种方法开始逐步应用于解决各种比较繁琐的科学研究,同样在管理人员中也是深受欢迎。这种方法存在2个优势:第一个优势是它非常的简单,无需一些繁琐的推算过程,通常情况下,一般的人就可以掌握;第二个优势便是非常的快速,正是由于这两大优势使得该法在现代项目管理中得到了广泛的使用。

  除此之外,上述方法还有着一项特别的优势,那就是它的适应性非常强,需要研究的问题约束少。这种方法的收敛性是指在概率的条件下的收敛,如果维数在增加,它的收敛速度也不会改变,而且存储也不占空间。因而在实际大的工程难题里,蒙特卡洛具有很强的优点。

  概率论中,来自同一族体,并且具有相同分布规律的n个互相独立的样本、、…,有相同的均值和方差,则:

  (3-4)

  设随机事件A发生的概率是P(A),n次独立事件里频数是m则随机事件发生频率:

  (3-5)

  且

  (3-6)

  Monte Carlo法的理论依据即上述两条定律,样本的均值对于族体是收敛的、事件的频率也是收敛于事件的期望,在构造的概率模型中通过抽样,进行模拟,只要抽取的数据整体足够大,两条定律成立精确的条件是样本点是无穷多个,显然,我们是做不到的,只有尽可能多的抽取样本点,而这又造成了计算量和时间的耗费,该方法仍然有需要改进的地方。

  Monte Carlo法的基本思路是根据联合密度函数由随机变量产生N个随机样本,将其代入功能函数,得到响应值与阀值的差,统计小于零的即失效的数据个数,计算出失效概率,注意,得到的失效概率只是估计值,并不完全准确。具体步骤如下:

  (1)确定随机变量的类型,它的分布形式以及相关参数,产生的N组随机向量。

  (2)把样本代入功能函数,求出功能函数的值。

  (3)求出失效概率估计值。

  (4)估计失效概率的变异系数cov。

  该方法应用十分广泛,没有对功能函数有特殊的要求,因此限制也较少,在不同方法求失效概率时一般把它作为比照的依据,验证其他的求解是否正确。从Monte Carlo法的基本思路也可以看出,编程能够较为轻松的实现,代码通常比较简单。只要保证抽样次数达到数万甚至几百万,它的精度在实际生产生活中是可以认可的。Monte Carlo法是基于概率论中的大数定理,但由于计算量大,所以在理论研究中应用更多,可以解决功能函数是非线性的隐式函数。

  下面用MCS对本文中的3个例子进行可靠性分析。

  3.2.1平面4杆桁架

  选取随机变量:杆的横截面积、弹性模量、承载

  失效极限位移:0.32mm

  随机变量分布:正态分布,标准差为均值的5%

  将其按照上述条件和蒙特卡洛法的步骤编写MATLAB程序:

  function[uu,disp,g,pf]=FP4bar(x)

  %x=[A E F1 F2]=[100 295000 20000-25000]随机变量

  N=1e4;%随机取样

  u=zeros(N,4);uu=u;%初始化

  for j=1:size(x,2)

  u(:,j)=normrnd(0,0.05,N,1);

  uu(:,j)=u(:,j).*x(j)+ones(N,1).*x(j);%返回原空间

  end

  disp=zeros(1,N);

  for j=1:N

  disp(j)=FEM4bar(uu(j,:));%调用节点位移

  end

  g=0.32-disp;%功能函数

  n=sum(g<0);%统计失效的次数

  pf=n/N;%失效概率

  end

  在MATLAB软件中运行得到失效概率:pf=3.04%

  3.2.2平面61杆桁架

  选取随机变量:杆的横截面积、弹性模量、杆长、承载

  失效极限位移:0.012m

  随机变量分布:正态分布,标准差为均值的5%

  将其按照上述条件和蒙特卡洛法的步骤编写MATLAB程序:

  function[uu,disp,g,pf]=FP61bar(x)

  %x=[A E l F]=[1e-4 2.1e11 1 1000]随机变量

  N=1e4;%随机取样10000个

  u=zeros(N,4);uu=u;%初始化

  for j=1:size(x,2)

  u(:,j)=normrnd(0,0.05,N,1);

  uu(:,j)=u(:,j).*x(j)+ones(N,1).*x(j);%返回原空间

  end

  disp=zeros(1,N);

  for j=1:N

  disp(j)=FEM61bar(uu(j,:));%调用节点位移

  end

  g=0.012-disp;%功能函数

  n=sum(g<0);%统计失效的次数

  pf=n/N;%失效概率

  end

  在MATLAB软件中运行得到:失效概率pf=2.732%

  3.2.3空间132杆桁架

  选取随机变量:杆的横截面积、弹性模量、杆长、承载

  失效极限位移:18.8mm

  随机变量分布:正态分布,标准差为均值的2%

  将其按照上述条件和蒙特卡洛法的步骤编写MATLAB程序:

  function[uu,disp,g,pf]=FP132bar(x)

  %x=[205.8-16.*ones(1,37)]

  N=1e4;%随机取样

  u=zeros(N,38);uu=u;%初始化

  for j=1:size(x,2)

  u(:,j)=normrnd(0,0.02,N,1);

  uu(:,j)=u(:,j).*x(j)+ones(N,1).*x(j);%返回原空间

  end

  disp=zeros(1,N);

  for j=1:N

  disp(j)=FEM132bar_pf(uu(j,:));%调用节点位移

  end

  g=18.8-disp;%功能函数

  n=sum(g<0);%统计失效的次数

  pf=n/N;%失效概率

  end

  在MATLAB软件中运行得到:失效概率pf=1.953%

  3.3.AK-MCS法

  3.3.1克里金理论

  克里金方法是一个先进的统计学过程,从一组z值的分散点中产生面积估计数。与内插工具集中的其他内插方法不同,克里金工具的有效使用意味着研究z值现象的空间行为之间的相互作用,然后选择最佳估计方法来产生输出面积。Idw插值工具和spline插值工具被称为确定性插值方法,因为这些方法直接基于周围的测量或一个特定的数学公式来确定产生的表面的柔度。第二种插值方法是地理统计学方法(如克里金方法),它基于包括自相关性(即测量点之间的统计关系)的统计学模型。因此,地理统计方法不仅可以产生预测的表面,而且还可以提供预测的确定性或准确性的某种衡量标准。

  Kriging作为一种半参数化的插值模型,该模型可对某一研究点相关的信息实现该点的响应等的模拟,这种模拟方式适宜加权的线性组合来完成的。而该模拟方式是采用最小化模拟值的误差的均方差进行确认的。该Kriging模型在线性无偏估计中,一直认为是最好的。同响应面法等模型方相比较,它无需构建特殊的数学模型,使得该模型的应用极为的方便。克里金方法假定取样点之间的距离或方向可以反映空间相关性。在地质领域,这种方法的应用是相当宽的。

  现阶段,通常比较典型的基于Kriging模型的可靠性分析方法,其原理基本是运用了Kriging模型预测的随机这一特点,由此发展成基于Kriging模型和学习函数的可靠性分析方法。现如今,越来越多的科学工作都对其进行了研究,也得到了不同类型的学习函数的概念,这些当中,比较多的是学习函数有UlSS]和EFF[s61函数。该种模型方法有着非常高的效率,在相关领域的应用得到肯定。Kriging模型是一种精确地插值模型,多应用于上世纪九十年代计算机领域。建立的Kriging模型可以预测出均值和方差等。它假设通过随机过程实现功能函数,Kriging模型包括一个回归模型和一个随机误差:

  (3-7)

  其中,是均值等于零的平稳高斯过程。由协方差函数确定:

  (3-8)

  式中l是相关系数的参数向量。非均匀的平方指数模型:

  (3-9)

  正态随机变量分布:

  (3-10)

  其中,

  (3-11)

  是根据得到的相关矩阵,试验点和预测点之间相关系数的向量,,是回归矩阵,。广义最小二乘法得到:

  (3-12)

  克里金模型的构建可分为两个阶段,分别是初始模型的构建和后续模型的更新,第一个阶段可以通过抽样如拉丁抽样、随机抽样、正交试验等获得初始的样本点,构建初始的克里金模型,接下来,则可根据原始的样本构建原始Kriging模型。第二阶段可根据不同加点原则向原始样本中逐步添加有非常关键作用的样本,并对建成Kriging模型进行实时更新,直到符合相关的收敛要求。据此在逐步添加样本的不同意义,此时可把加点规则分为两类。一通常是着重说明代理模型的精度,同时这个原则表明只要准确度好,由此得到的模型分析的结果失效概率即可符合相关要求。在第一类加点原则中,常用的学习函数有EGO(Efficient Global Optimization)、EFF(Expected Feasibility Function)、U和H等。EGO函数被用来描述后续增加样本对Kriging模型改进精度的期望。EFF函数通常用来评判添加样本处的真实功能函数值在某一范围内趋于零的趋势。U函数的功能是判断最有可能被错误判定的样本点。H函数的功能则是判断样本预测值的存在的不确定性。二是直接注重失效概率精度。然而两类加点原则相比,后者学习的相关的函数相对较少,已提出的学习函数有LIF(Least Improvement Function)和Zhu等推导的学习函数。LIF可以用来衡量增加样本对失效概率精度的改善程度。由Zhu等提出的学习函数可用来识别对失效概率误差影响最大的样本。

  3.3.2 DACE工具箱

  现阶段,若要构建Kriging模型,可以直接运用MATLAB工具箱下的DACE,本论文得到程序就是按照这种思想去编写的。这里简明而要的介绍一下运用该工具箱建立Kriging模型的一些重要流程:首先根据试验设计方法(常用的为拉丁超立方抽样)产生出设计变量的样本集S及其响应值L设定参数的初始值(同时给定一个上下限),然后设定出回归多项式函数的形式和相关函数的形式。这个工具箱的典型的用途是,根据计算机经验建立一个基于克里金理论的近似模型。近似模型是计算机模型通常用到的另一个别称。而计算机经验则是通过计算机模型的两个输入及响应系列的整合获取的,该经验通常情况下是高维数的,而该模型则是确定的,因此模型的答案会产生相同的错误,也就是说,模型的答案是相同的。选择模型往往是解决设计问题必不可少的,因此需要更好地选择计算机模型,例如,计算机建立的模型是否适合于实际数据,计算机模型中能更好的预测问题。

  模型函数

  给定一组设计点及其响应值,试验点满足相应的正态分布。

  (3-13)

  式中,为矩阵X第j列给出的向量,和分别是均值和方差。

  我们采用模型y来表示n维输入的确定性响应y(x),作为回归模型F和随机函数(随机过程)的实现,

  (3-14)

  用一个p个函数组成的线性结合的回归模型:

  (3-15)

  式中,是回归参数,假设随机过程z的均值为零,协方差为零。

  (3-16)

  在z(w)和z(x)之中,是过程方差,组件的响应和相关模型与参数是。真实值可以表示为:

  (3-17)

  式中是近似误差。

  调用格式:[dmodel,perf]=dacefit(S,Y,regr,corr,theta0)[dmodel,perf]=...

  dacefit(S,Y,regr,corr,theta0,lob,upb)

  其中,输入参数:

  S,试验点m×n的矩阵。

  Y,试验点S的响应值或响应矩阵m×q。

  theta0,如果lob和upb并不存在,那么theta0应持有相关函数参数。否则theta0应持有初始猜测。

  lob,upb可选填的。如果存在,那么他们应该长度相同的向量分别theta0和应持有低和的上界。

  如果他们不存在,则由theta0中的值。

  预测函数

  对于设计点的集合S,我们有扩展的设计矩阵,

  (3-18)

  预测函数predictor可以预测单个输出或者多个输出,可以一次性的预测若干组的输入变量,可以节省大量的时间。

  调用函数的格式:

  y=predictor(x,dmodel)

  [y,or]=predictor(x,dmodel)

  [y,dy,mse]=predictor(x,dmodel)

  [y,dy,mse,dmse]=predictor(x,dmodel)

  如果m=1和n>1,那么行向量和列向量都可以作为输入x。否则x必须是一个m×n的数组,其中的位置按行存储。

  其中,y是预测的响应值,mse是预测值的方差。

  回归模型调用格式:

  f=regpoly0(S)

  [f,df]=regpoly0(S)

  得到一阶多项式的值。

  f=regpoly1(S)

  [f,df]=regpoly1(S)

  得到二阶多项式的值。

  f=regpoly2(S)

  [f,df]=regpoly2(S)

  得到三阶多项式的值。

  相关模型

  工具箱提供了七个实现模型的函数。

  调用格式:r=correxp(theta,d)

  [r,dr]=correxp(theta,d)

  theta,相关函数中的参数。允许使用标量值。这对应于一个各向同性模型:所有等于。否则,元素的个数必须等于d中的维数n。

  3.3.3学习函数

  Kriging模型输入量得到对应预测响应值,以及预测响应值的方差,这个方差可以用来表征预测响应值的局部不确定性,一般方差值越大,说明预测响应值越不确定。假若人们将某个拥有高方差的预测响应值对应的样本点添加到试验设计样本集中,这种模型的拟合精度将会得到很好的改变,上述过程换句话说就是主动学习过程。现如今,此类根据主动学习函数的加点准则已经得到了普遍的接受,也有了一定的发展,而且它通常的功能即确定需要添加的最佳样本点,接着将其加至试验设计中,对Kriging模型进行实时更新,经过循环,一直到符合相关精条件的Kriging模型。当经过试验得到该模型后,便开始对它的精度进行分析验证,如果原始样本点较少,将会使得该模型的精度偏低,不能满足设计要求,为了解决这中问题,常规的方法是在试验设计时就生成尽量多的初始样本点,这样在一次试验设计之后得到的Kriging模型的模拟精度基本就能满足要求,很明显,由于所取初始样本点过多,这种方法会是一种效率很低的方法。在这之后,相关研究者通过对序列迭代法进行探索,一种新的方法便推出,新方法可以用较少的样本点就能够去模拟需要符合精度要求的Kriging模型,通过一定的验证,该方法行之有效,此外该方法对原始样本点的限制条件并不高,所以,该方法的应用也是越来越普及。这种序列迭代方法的基本思想是首先根据试验设计产生出很少的初始样本点,并且计算出其对应的响应值,然后构建出Kriging模型,之后对模型进检验,假若不满足迭代终止要求,则根据加点准则,产生出一个添加样本点集,将这个添加样本点集添加到原始试验设计之中,再次构建Kriging模型,如次迭代更新下去,直至满足终止要求。不难发现,此种序列迭代方式对该模型的更新,是根据序列添加样本点来实现的,依据加点规则,只需要增添对该模型拟合精度起关键意义的样本点,一些不关键的样品本点可以直接忽略,因此它的效率会更显著。基于性能函数表示法在蒙特卡罗仿真中很重要的主张,可以提出基于不同于EFF概念的学习函数。事实上,在AK-MCS中,符号的精度只需要在蒙特卡罗族的点之间,根据随机变量的分布,极限状态可以在低概率密度的配置中近似。应在实验设计中增加高风险点,并对性能函数进行评估。事实上,这些点的不确定性会导致预测值从正变为负(或从负变为正)。然后失效概率发生变化。潜在的“危险”点可能有三个特征:接近极限的状态,具有重要的不确定性(高克里金方差),或者两者兼而有之。为了识别它们,我们提出了学习功能U(x):

  (3-19)

  表示预测极限状态和估计极限状态之间的距离,用Kriging标准差表示。这是一个可靠性指标,它考虑了如果使用相同的符号就会发生的符号错误的风险。在可靠性分析中,蒙特卡洛样本点分为两类功能函数为负数,即结构失效和功能函数是正数,即结构处于安全区域。学习函数的值越小说明该处点的不确定性很大,即很有可能不安全,那么分辨极限函数的正负号的预测很有可能使可靠性分析结果有偏差;学习函数的值越大,说明这个指数U(x)可能与Cox和John提出的优化的下置信限函数(LCB)有关。事实上,LCB的定义是:

  (3-20)

  首先,选项b等于2或2.5。这个选择取决于对局部改进或更完整的改进的期望。然后计算并最小化LCB。由于这是一个可靠性问题,而不是优化问题,所以执行是相反的。LCB不只是作为预期事件的符号计算,因此定义为LCB=0。然后计算整个总体(第5阶段)的U(x)=b(第5阶段)。对于xS,学习标准(步骤5)是。

  阶段6需要一个停止条件。对于U,定义为xS的最小值。作为指标的可靠性,U=2对应符号中相应的概率。停止条件。

  Cox和John提出了优化算法的值为2。此外,如下面的验证示例所示,这种停止条件保证了高度的准确性。

  可以看出,学习函数U对接近预测极限状态的点的权重大于克里金方差较大的点。所以,此种学习函数与EFF函数相比,在近似估计故障概率方面是相对较慢的,然而当收敛到终止的情况下,则会相对较快。学习函数U的挑选可参照自然规律,非常简单的理解U值与克里金方差和预测直接相关。除了以不同的方式进行优化外,Juang等人还提出了一种近似的概念。

  3.3.4 AK-MCS法

  AK-MCS法是一种结合了蒙特卡洛和克里金模型的主动学习方法,它能有效减少调用函数的次数,减少计算时间。Kriging模型同相关模型进行比较时,可以发现Kriging模型在“统计性”方面更加突出。另一方面,这种模型在有效性方面,可完全不考虑随机误差,换句话说,现有信息中有没有噪声信息对该模型的有效性程度,是不会产生干扰的。然而在可靠性分析中,一些相一致的结构数据肯定会实现对应一致的结构响应,在处理这种情况时,该模型的精度则会更有优势。此外,Kriging模型若是半参数化的,它根本不用再次建立新的特殊数学模型,综合对比,该模型无论是在灵活性还是简便性方面,都有着独特的优势。Kriging是由一个参数模型和一个非参数随机过程联合构成的。它比单个的参数化模型更具有灵活性,同时又克服了非参数化模型处理高维数据存在的局限性,比单个的参数化模型具有更强的预测能力。与通常的插值技术相比较,Kriging模型有很多优势。首先,它是基于已有信息的动态构造,也就是仅运用了估计点附近的相关信息,而非用整个的信息对不知的信息进行模拟。其次,它兼有局部和全局两种统计优势,可以使得该模型能够计算分析现有信息的方向和情况。该模型是线性回归分析中一种优化模型,涵盖了线性回归和非参数两个部分,然而这里面非参数的部分可以以随机分布的方式来实现。受限制较少。

  利用DACE工具箱建立Kriging模型并编写MATLAB程序。它没有传统方法中可靠度计算和结构分析交替进行的过程。是一个先结构分析后可靠度计算的过程,可以最大程度的利用现有的大型工程分析软件,更加高效的计算相应的难题。

  AK-MCS的具体步骤是:

  (1)在设计空间中生成蒙特卡洛随机样本点作为候选点。将它名为S,由个点组成。在这一阶段,没有对任何性能函数进行评估。如果主动学习需要,他们表示待评估的候选点。除非达到第9阶段,否则在AK-MCS的整个学习过程中,这一总体始终保持不变。

  (2)初步实验设计(DoE)的定义。为了执行克里金法,需要进行实验设计。在此,它候选点S中的个点的随机选择组成。这些点在性能函数上进行评估,并用作Kriging模型的实验的初始设计。根据经验,十几个点即可。实验的初始设计最好定义得小一些,并且逐步地只添加最能改进元模型的点。随着问题规模的增大,实验的初始设计要大得多。

  (3)根据实验设计计算克里金模型。此阶段由工具箱DACE执行,因为该工具已在多个Kriging参考文章中使用。选择相关模型为高斯模型,并将回归模型取常数(普通克里金模型)。

  (4)通过克里金法进行预测并估计失效概率。首先,式(10)中的克里金预测值(i=1,……,)由DACE得到。然后,根据这些预测的符号来估计失败的概率。通过负或零克里金预测的S中点数与式(1)中S中点数总数的比值得到:

  (3-21)

  (5)确定S中下一个最佳点以评估性能函数,这一阶段需要实现我们称为学习函数的东西。该函数针对S的每个点进行计算。由于它取决于Kriging预测和方差,因此G的Kriging方差(i=1,……,)首先需使用式(11)计算。在本文中,我们使用学习函数。我们实现了一个命名为U的新函数。然后,然后根据学习函数值确定下一个评价性能函数的最佳点。识别的准则称为学习准则,它取决于学习函数的选择。

  (6)学习停止条件。一旦用学习准则确定了S中要评估的下一个最佳点,就将其对应的学习函数的值与停止的值进行比较。这称为停止条件。同样,这取决于讨论的学习函数的选择。

  (7)以最佳点更新先前的实验设计。如果不满足阶段6的停止条件,则继续学习,并根据真实性能函数评估S中的最佳点。之后,将其添加到实验设计中:。然后,该方法返回到第3阶段,以使用由N个点组成的更新实验设计来计算新的Kriging模型。

  (8)计算失效概率的变异系数。如果满足阶段6的停止条件,则学习将停止,并且在性能函数的点符号上,元模型被称为足够准确。下一阶段包括检查蒙特卡洛族S是否足够大,以便在失效概率的Kriging估计上有较低的变异系数(阶段4)。在这里,允许低于5%的变化系数。它的计算方法与公式(2)相同。但关于克里金结果:

  (3-22)

  (9)更新蒙特卡洛随机样本点。如果发现变异系数太高,则S会随着来自另一个蒙特卡洛族的新点而增加(类似于阶段1)。AK-MCS返回到阶段4,以预测新的总体,并且主动学习方法继续进行,直到再次满足停止条件为止。重要的是要注意,不会丢失有关先前性能函数评估的任何信息。

  (10)AK-MCS的终止。如果的变化系数足够低,则AK-MCS停止,并且将失效概率的最后估计视为该方法的结果。

  图3.1 AK-MCS流程图

  3.3.5平面4杆桁架结构实例

  该结构的随机变量仍然是截面积、所受载荷和弹性模量,与蒙特卡洛方法中一样,并且其分布规律仍是正态分布。

  将4杆桁架的待测点从蒙特卡洛抽样的到的点作为输入,并选择12个点作为试验点,求出功能函数值,添加到DACE工具箱中,开始建立Kriging模型,进而对预测点进行预测,从得到的预测的点中筛选出最佳的下一个点并把该点逐渐加入试验点重新建立模型直至学习函数U的值大于或等于2,计算出失效概率。

  编写MATLAB程序如下:

  function[pf]=AK_MCS_4bar

  clear;clc

  %%

  load 4bar_pf

  s=uu;%待测点

  load 4bar_lhs%12个试验点

  varNum=4;

  Y=0.32-y;

  S=X;

  %%建立模型并预测

  for i=1:1000

  theta=10.*ones(1,varNum);lob=(1e-10).*ones(1,varNum);upb=5.*ones(1,varNum);

  [dmodel,perf]=...

  dacefit(S,Y, regpoly0, corrgauss,theta,lob,upb)

  [YX,MSE]=predictor(s,dmodel);

  N=size(s,1);

  pf=sum(YX<0)/N;

  %学习函数

  MSE=MSE.^0.5;

  U=abs(YX)./MSE;

  [m,p]=min(U);

  best_x=s(p,:);

  if m>=2%终止条件

  break

  end

  S=[S;best_x];

  Y=[Y;0.32-FEM4bar(best_x)];

  cov(i)=((1-pf)/(pf*N))^0.5;

  end

  end

  AK-MCS方法与蒙特卡洛方法的对比如表3-1所示:

  表3-1 AK-MCS方法与蒙特卡洛方法对比

  Monte Carlo AK-MCS

  失效概率3.039%3.15%

  调用次数100000 18

  变异系数1.786%1.753%

  图3-2 AK-MCS方法调用的功能函数

  从表3-1可知,AK-MCS方法仅调用函数18次即得到了失效概率,而蒙特卡洛方法则调用了十万次函数,且AK-MCS方法的失效概率精确度不亚于蒙特卡洛方法得到的失效概率,AK-MCS方法更加高效快捷。图3-2是AK-MCS方法调用的功能函数图,其中,前12个点是初始设计点,之后的几个点均在响应值为0附近,AK-MCS方法的学习过程正是在失效临界处模拟的,具体的变量与响应值见附录2。

  3.3.6平面61杆桁架结构实例

  与上一节类似,在编写程序如下:

  function[pf]=AK_MCS_61bar

  clear;clc

  %%u待测点

  varNum=4;

  load X_61bar%试验点

  load uu_61bar%待测点

  load y_61bar

  N=10000;

  Y=0.012-y;S=X;

  %%建立模型并预测

  for i=1:1000

  theta=10.*ones(1,varNum);lob=(1e-10).*ones(1,varNum);upb=5.*ones(1,varNum);

  [dmodel,perf]=...

  dacefit(S,Y, regpoly0, corrgauss,theta,lob,upb)

  [YX,MSE]=predictor(uu,dmodel);

  pf=sum(YX<0)/N;

  %学习函数

  MSE=MSE.^0.5;

  U=abs(YX)./MSE;

  [m,p]=min(U);

  best_x=uu(p,:);

  if m>=2%终止条件

  break

  end

  S=[S;best_x];

  Y=[Y;0.012-FEM61bar(best_x)];

  cov(i)=((1-pf)/(pf*N))^0.5;

  end

  end

  AK-MCS方法与蒙特卡洛方法的对比如表3-2所示:

  表3-2 AK-MCS方法与蒙特卡洛方法对比

  Monte Carlo AK-MCS

  失效概率2.732%2.730%

  调用次数100000 33

  变异系数1.887%1.89%

  图3-3调用的功能函数值

  从表3-2可知,AK-MCS方法仅调用函数33次即得到了失效概率,而蒙特卡洛方法则调用了十万次函数,且AK-MCS方法的失效概率精确度不亚于蒙特卡洛方法得到的失效概率,AK-MCS方法更加高效快捷。图3-3是AK-MCS方法调用的功能函数图,其中,前12个点是初始设计点,之后添加的训练点均在响应值为0附近,AK-MCS方法的学习过程正是在失效临界处模拟的,具体的变量与响应值见附录3。

  3.3.7空间132杆桁架结构实例

  基本的思路过程仍与前两个实例相同,程序如下:

  function[pf]=AK_MCS

  clear;clc

  %%

  load uu_mcs

  s=uu;%待测点

  load u_lhs%试验点

  load disp_lhs

  varNum=38;

  Y=18.8-disp;

  S=X;

  %%建立模型并预测

  for i=1:1000

  theta=10.*ones(1,varNum);lob=(1e-10).*ones(1,varNum);upb=5.*ones(1,varNum);

  [dmodel,perf]=...

  dacefit(S,Y, regpoly0, corrgauss,theta,lob,upb)

  [YX,MSE]=predictor(s,dmodel);

  N=size(s,1);

  pf=sum(YX<0)/N;

  %学习函数

  MSE=MSE.^0.5;

  U=abs(YX)./MSE;

  [m,p]=min(U);

  best_x=s(p,:);

  if m>=2%终止条件

  break

  end

  S=[S;best_x];

  Y=[Y;18.8-FEM132bar_pf(best_x)];

  cov(i)=((1-pf)/(pf*N))^0.5;

  end

  AK-MCS方法与蒙特卡洛方法的对比如表3-3所示:

  表3-3 AK-MCS方法与蒙特卡洛方法对比

  Monte Carlo AK-MCS

  失效概率1.953%1.90%

  调用次数100000 112

  变异系数2.241%2.272%

  图3-4调用的功能函数

  从表3-3可知,AK-MCS方法仅调用函数112次即得到了失效概率,而蒙特卡洛方法则调用了十万次函数,且AK-MCS方法的失效概率精确度不亚于蒙特卡洛方法得到的失效概率,AK-MCS方法更加高效快捷。图3-4是AK-MCS方法调用的功能函数图,其中,前12个点是初始设计点,之后添加的训练点均在响应值为0附近,AK-MCS方法的学习过程正是在失效临界处模拟的。

  3.4本章小结

  本章介绍了可靠性的基本理论和一些求解可靠度的方法并比较了优缺点,详细介绍了蒙特卡洛法和基于克里金模型的主动学习方法(AK-MCS),针对蒙特卡洛法和AK-MCS方法求解了前文的三个桁架结构,二者的结果精度均较好,蒙特卡洛法的精度和调用次数有关,在相近的精度下,AK-MCS的的对有限元求解模型的调用次数低得多。

  第4章总结与展望

  本论文利用MATLAB软件编写了桁架结构的可靠性分析,根据一种结合了蒙特卡洛和克里金模型的主动学习方法,分析了3个桁架结构的可靠性。结果表明,该方法在可靠性的精确度上是可以被信任的,其稳定性也较好。AK-MCS避免了广泛的蒙特卡洛带来的计算耗时,工作量大的问题。它选取了数十个点调用隐式功能函数,建立初始模型并预测总体中的其他点的均值和响应值,通过学习函数U筛选出最佳的点加入Kriging模型的构造中,不断如此往复,直至达到停止条件。其结果的准确度仍然可观,与蒙特卡洛的结果仅相差很小。虽然本文中的例子不能完全概括所有的桁架结构,但可以看出对不具有解析式的非线性隐式函数的桁架,这种方法有很好的实用性,减少了实际中不必要的耗费资源的情况。

  本文在进行论证时,在构建Kriging模型时运用到了MATLAB中提供的DACE工具箱,因为它受到一定的约束,所以只对表明的方法可以解决现阶段可靠性分析方法受限于功能函数非线性程度以及功能函数为隐式的问题进行了论证,未思考现阶段可靠性分析方法的存在的其他问题,也就是对设计变量的维度的干扰,因此可作为后续研究的方向。此外,AK-MCS方法的收敛条件较为严苛,在实际运用中受限制较多,目前已有更为先进的方法分析可靠性,对收敛条件的要求也更低。综上所述,期望在后续的研究中,对其存在的不足,通过一定的方式进优化同时完善设计出一种可靠性分析及可靠性优化方法,此外,可以通过工程实践,采用实验的思路对所涉及的理论方法进行实践论证。

  致谢

  本论文是在导师杨xx老师的悉心指导下完成的。导师具有深厚的专业知识技能,在我遇到困难时给予耐心的指导,很快的回复答疑的问题,在疫情期间虽然相隔甚远,仍让我深深感受到了老师的敬业和无私。最终,我能够顺利完成毕业论文,在此表示向老师我最诚挚的谢意。通过本次论文的完成,我体会到了周围人对我的无私帮助,在平时生活里也许感受不深,在此,向你们致以我深深的感激之情。

  参考文献

  [1]吕震宙,宋述芳,李洪双,袁修开.结构机构可靠性及可靠性灵敏度分[M].北京科学出版社,2009.

  [2]B.Echard,N.Gayton,M.Lemaire.AK-MCS:An active learning reliability method combining Kriging and Monte Carlo Simulation[J].Structural Safety,2011,33(2).

  [3]A.M.Freudenthal.The safety of structures[J].Transitition ASCE,1947.

  [4]Comell C.A.A probability-based structural code[J].Journal of ACI,1969,66(2):15-25.

  [5]Rackwitz R,Fiessler B.Structural reliability under combined random load sequences[J].Computers and Structures,1978,9(5):489-494.

  [6]万越,吕震宙,袁修开.改进的重要抽样可靠性灵敏度估计及其方差分[J].机械设计,2008,25(1 2):69-78.

  [7]H.J.Pradlwarter,G.I.Schueller.Application of line sampling simulation method to reliability benchmark problems[J].Structural Safety,2007(29):208-221.

  [8]程进.基于响应面法的几何非线性结构概率响应分析[J].同济大学学报(自然科学版),2006,34(9):1147-1151.

  [9]谭晓慧,王建国,刘新荣.改进的响应面法及其在可靠度分析中的应用[J].岩石力学与工程学报,2005.

  [10]吕震宙,杨子政.基于神经网络的可靠性分析新方法[J].机械强度,2006,28(5):699-702.

  [11]Echard B,Gayton N,Lemaire M.A combined importance sampling and Kriging reliability method for small failure probabilities with time-demanding numerical models[J].Reliability Engineering and System Safety,2013,111:232-240.

  [12]Makoto OHSAKI.Structural Optimization for Specified Nonlinear Buckling Load Factor Japan[J].Indust.Appl.Math.,19(2002),163-179.

  [13]Makoto Ohsaki,Kiyohiro Ikeda.Stability and Optimization of Structures Generalized Sensitivity Analysis[M].February 2007,143-151.

  [14]徐宇豪.Kriging模型在可靠性及可靠性优化设计中的应用研究[D].东北大学,2017.

  [15]黄晓旭,陈建桥.基于主动学习Kriging模型的可靠性分析[J].固体力学学报,2016,37(02):172-180.

  [16]孙志礼,李瑞,闫玉涛,王健.一种用于结构可靠性分析的Kriging学习函数[J].哈尔滨工业大学学报,2017,49(07):146-151.

  [17]唐恺.基于Kriging方法的结构可靠性优化设计[D].沈阳航空航天大学,2017.

  [18]刘阔,李晓雷,王健.一种基于Kriging模型的机械结构可靠性分析方法[J].东北大学学报(自然科学版),2017,38(07):1002-1006.

  [19]米彩盈.铁道机车车辆结构强度[M].西南交通大学出版社,2007:7-26.

  [20]徐斌,,高跃飞等.MATLAB有限元结构动力学分析与工程应用[M].清华大学出版

  社,2009:69-76.

  [21]徐斌,高跃飞,余龙.MATLAB有限元结构动力学分析与工程应用[M].北京清华大学学研大厦:清华大学,0912.69-76.

  [22]Byeng D.Youn,Kyung K.Choi.A new response surface methodology for reliability-based design optimization[J].Computers and Structures,2004,82:24 1.256.

  [23]J.Tu,K.K.Choi.A new study on Reliability—Based Design Optimization[J].Journal of Mechanical Design,1999.

  [24]Liang Zhao,K.K.Choi,and Ikjin Lee.Response Surface Method Using Sequential Sampling for Reliability-Based Design Optimization[J].ASME,2009.

  [25]阎宏生,胡云昌,牛勇.基于神经网络响应面的结构可靠性分析方法研究[J].海洋工程,2002,20(2):1-6.

  [26]Fan Li,Teresa Wu,Adedeji Badiru.A single-loop deterministic method for reliability-based design optimization[J].Engineering Optimization,2012,45(4):435-458.

  [27]Chen Zhenzhong.An adaptive decoupling approach for reliability-based design optimization[J].Computers and Structures,2013,117:58-66.

  [28]Barton J.Bichon,Sankaran Mahadevan,Michael S.Eldred.Reliability-Based Design Optimization using Efficient Global Reliability Analysis[J].American Institute of Aeronautics and Astronautics,2009.

  [29]Harish Agarwal,John Renaud.Reliability based design optimization using response surfaces in application to multidisciplinary systems[J].Engineering Optimization,2004,36(3):291-311.

  [30]史妍妍,孙志礼,闫明.基于响应面法的隐式极限状态函数可靠性灵敏度分析方法[J].机械科学与技术,2009,28(5):648-651.

  [31]程进.基于响应面法的几何非线性结构概率响应分析[J].同济大学学报(自然科学版),2006,34(9):1147-1151.

  [32]张崎,李兴斯.结构可靠性分析的模拟重要抽样方法[J].工程力学,2007,24(1):24-37.

  [33]Li G Cheng GD.Optimal decision for the target value of performance based structural svstem reliability,Structural Multidiseiplinary Optimization,2001.22:261-267

  [34]Costa JP,Pronzato L and Thierry E.A comparison between kriging and radial basis function networks for nonlinear prediction.Darts Proc.NSIP'99,Antalya,June 1999.

  [35]Lee Tw'Kwak BM.A reliability based optimal design using advanced first order second moment method.Mechanics Structural Machines,1987,15(4):523-542.

  [36]Kirjner-Neto C,Polak E,Kiureghian AD.An outer approximations approach to reliability based optimal design of structures.Journal of Optimization Theory and Applications,1998,98(1):1-16.

  [37]李兴斯,钱令希.基于概率极限状态的结构优化设计,计算结构力学及其应用,1996,13(4):379-384.

  [38]李为吉,基于可靠性的结构优化设计方法,西北工业大学学报,1989,7(2):147-154.

  [39]Sakata S,Ashida F'Zako M.Structural optimization using kriging approximation.Computer Methods inApplied Mechanics and Engineering,2003,192:923-939.

  [40]Bennea RM,Structural analysis methods for system reliability,Stmctural Saf色咄1990,7:109-114.

  [41]孙海虹.结构可靠性分析改进的重要抽样法.武汉交通科技大学学报,1994,Vol 18,No.3:24l-246.

  [42]Baker E Probability estimation and information principles,Structural Safety,1990,9:97-116.

  [43]thames A,Cooper WW.Chance constrained programming.Management Science,1959,6:73-79.

  [44]Kiuraghian AD,Polak E.Reliability Based Optimal Design:A Decoupled Approach,Reliability and Optimization of Structural Systems,A.S.Nowak ed.Book Crafiers,Chelsea.Mich.197-205.

  [45]吕震宙,冯蕴雯.结构可靠性问题研究的若干进展,力学进展,2000,30(1):21.28

  [46]许林,基于可靠度的结构优化研究,博士学位论文,大连理工大学,2004

  [47]李兴斯,钱令希.基于概率极限状态的结构优化设计,计算结构力学及其应用,1996,13(4):379-384.

  [48]Kiuraghian AD,Polak E.Reliability Based Optimal Design:A Decoupled Approach,Reliability and Optimization of Structural Systems,A.S.Nowak ed.Book Crafiers,Chelsea.Mich.197-205.

  [49]Moses F'Kinser DE.Analysis of structural reliability.Journal of Structural Division,1 967,93(5):147-164.

  [50]徐军,张利民,郑颖人.基于数值模拟和BP网络的可靠度计算方法,岩石力学与J一程学报,2003,22(3):395.

  [51]杨林德,徐超.Monte Carlo模拟法与基坑变形的可靠度分析,岩土力学,1999,20(1):15-18.

  附录

  附表1 132杆空间桁架的MATLAB有限元程序节点位移

  节点编号x(mm)y(mm)z(mm)

  1-4.86E-15 1.68E-15-17.32245

  2-2.22243-5.88E-07-15.86256

  3-3.772694 3.51E-06-11.44666

  4-1.874815-1.082427-11.68494

  5-2.69572 1.19E-05-4.913355

  6 0.5028519 0.1872396-3.468763

  7 0.4135441 0.3418438-3.468808

  8 0 0 0

  9 0 0 0

  10 0 0 0

  11 0 0 0

  12-1.111215-1.924681-15.86256

  13-1.88635-3.267247-11.44666

  14 1.85E-06-2.164851-11.68494

  15-1.34787-2.334556-4.913355

  16 0.0892717 0.5291023-3.468763

  17-0.089273 0.5290616-3.468808

  18 0 0 0

  19 0 0 0

  20 0 0 0

  21 0 0 0

  22 1.1112156-1.924681-15.86256

  23 1.886344-3.267251-11.44666

  24 1.8748171-1.082424-11.68494

  25 1.3478497-2.334568-4.913355

  26-0.41358 0.3418627-3.468763

  27-0.502817 0.1872178-3.468808

  28 0 0 0

  29 0 0 0

  30 0 0 0

  31 0 0 0

  32 2.2224302 5.88E-07-15.86256

  33 3.772694-3.51E-06-11.44666

  34 1.8748152 1.0824272-11.68494

  35 2.69572004986617-1.18784495782069e-05-4.91335544781928

  36-0.502851896745587-0.187239556186615-3.46876302946966

  37-0.413544116542003-0.341843817163346-3.46880791224364

  38 0 0 0

  39 0 0 0

  40 0 0 0

  41 0 0 0

  42 1.11121458036445 1.92468128810960-15.8625598358957

  43 1.88635004942334 3.26724709972182-11.4466574629593

  44-1.85102103884907e-06 2.16485122409777-11.6849423301073

  45 1.34787031197218 2.33455610545037-4.91335544781928

  46-0.0892717361218448-0.529102295016124-3.46876302946972

  47 0.0892733715191155-0.529061619092640-3.46880791224363

  48 0 0 0

  49 0 0 0

  50 0 0 0

  51 0 0 0

  52-1.11121559950924 1.92468069970611-15.8625598358958

  53-1.88634396408833 3.26725061309138-11.4466574629590

  54-1.87481708099302 1.08242400901763-11.6849423301073

  55-1.34784973789402 2.33456798390000-4.91335544781932

  56 0.413580160623662-0.341862738829489-3.46876302946975

  57 0.502817488061144-0.187217801929291-3.46880791224361

  58 0 0 0

  59 0 0 0

  60 0 0 0

  61 0 0 0

  附表2 4杆桁架调用的设计点

  试验点编号试验点功能函数值

  1 113.5507 313461.8 19172.62-27141.9 0.1045399

  2 102.4358 332968.4 21786.43-25282.5 0.0644998

  3 103.9908 327305.3 17697.6-28124.5 0.1120177

  4 93.11279 292262.2 22639.56-28457.4-0.012771

  5 92.42501 277325.8 21197.06-23080.8-0.010793

  6 98.81337 257511.6 19753.94-25753.7 0.0094716

  7 96.17996 307864.6 20561.89-21802.8 0.0422342

  8 110.3609 317171.6 17242.24-22001 0.1229647

  9 86.12373 280549 20490.89-24958-0.019226

  10 88.88934 270185.9 22435.77-26382.9-0.05367

  11 106.5916 258514.3 18505.87-23476.7 0.0513653

  12 109.1203 299222.3 18281.7-23944.2 0.0960365

  13 97.13405 277596 21706.31-22470.4-0.002004

  14 100.5596 254151.4 20449.97-23566.2-6.38E-05

  15 94.76197 267409.3 20265.73-26569.2 0.0001019

  16 90.45244 272807.9 19746.27-24800-8.68E-05

  17 96.63871 260446.5 20133.42-23556.6 3.10E-05

  18 86.38171 299422.1 20635.49-23559.8 0.000869

  附表3 61杆桁架调用的设计点

  试验点编号试验点功能函数值

  1 0.000107 2.04E+11 1.061385 1010.227 0.00181028

  2 0.000103 2.06E+11 0.983297 1052.345 0.0018763

  3 9.85E-05 2.09E+11 1.015563 1020.39 0.00153853

  4 0.00011 2.06E+11 1.046837 1101.372 0.00141075

  5 0.000104 2.01E+11 0.980435 1003.629 0.00224311

  6 9.78E-05 2.17E+11 1.17299 1042.502-7.94E-06

  7 9.86E-05 2.08E+11 0.981419 1049.702 0.00157375

  8 0.000104 2.1E+11 1.027917 901.7384 0.00320491

  9 9.89E-05 2.27E+11 0.991739 932.1579 0.00345446

  10 0.000104 2.13E+11 1.05892 1003.057 0.00203881

  11 0.000102 2.07E+11 1.022306 1113.653 0.00076919

  12 9.56E-05 1.95E+11 1.058238 1013.256 9.64E-06

  13 9.29E-05 2.03E+11 0.979207 1118.208-8.93E-05

  14 9.77E-05 1.75E+11 1.011371 967.7657 7.71E-05

  15 9.23E-05 1.95E+11 1.026752 1008.488 9.23E-06

  16 9.42E-05 1.87E+11 0.979755 1032.225 9.34E-05

  17 9.56E-05 2.03E+11 1.108939 1009.916 1.02E-05

  18 9.77E-05 1.93E+11 0.982119 1106.91-2.56E-06

  19 9.36E-05 2.08E+11 1.080059 1038.49-4.52E-06

  20 0.000104 1.95E+11 1.066189 1092.828-2.30E-05

  21 8.23E-05 1.99E+11 0.932491 1013.458 1.19E-05

  22 9.99E-05 1.92E+11 1.073997 1031.553-1.56E-06

  23 9.07E-05 2.05E+11 1.09415 981.2604-3.79E-06

  24 8.83E-05 1.97E+11 1.031835 971.4163 2.07E-06

  25 9.77E-05 1.88E+11 1.027742 1033.927-1.90E-07

  26 9.21E-05 1.97E+11 1.071481 976.3448-1.70E-08

  27 9.38E-05 1.92E+11 0.971696 1067.306-3.95E-07

  28 9.63E-05 2.08E+11 1.038676 1113.234-2.28E-06

  29 0.000106 1.8E+11 1.049526 1049.043 3.03E-06

  30 0.000104 2.17E+11 1.154667 1135.017-4.92E-06

  31 0.000101 1.66E+11 0.930796 1039.462-5.87E-07

  32 9.55E-05 2.02E+11 1.047389 1061.714 1.73E-06

  33 9.93E-05 2.1E+11 1.02616 1172.924-4.81E-07