工程结构可靠度计算的Matlab实现

　 计算机应用 　

工程结构可靠度计算的Matlab实现

章慧健

(西南交通大学土木工程学院,四川成都610031)

【摘　要】　从基本原理和相关算例分析两方面,阐述利用功能极其强大的科学和工程计算数学软件系统Matlab的优化工具箱和随机数发生器实现工程结构可靠度的多种计算方法。　　【关键词】　结构可靠度;　优化工具箱;　蒙特卡罗法;　随机变量;　Matlab软件　　【中图分类号】　TP319　　　　　　　　　　　　【文献标识码】　B

工程结构可靠度计算除了如何确定极限状态功能函数

和各变量的概型分布以外,很大程度上是一种数学计算(本文不讨论如何确定功能函数和变量的概型分布,只对已知功能函数如何计算可靠度指标β的情况作一论述)。Matlab是一种功能极其强大的科学和工程计算数学软件系统,大量数学、统计、科学和工程所需的函数、C程语言相比,Matlab运算功能强等特点,可以大大提高编程效率。

2　基于Matlab优化工具箱的最优化方法

211　,1,,n是n,,,见下例),由这些随机变量表示的结构极限状:Z=g(X1,X2,…,Xn)=0,根据可靠度指标β的几何意义(在标准正态坐标系中,原点到极限状态曲面的最短距离)可知:

β=

1　　　可靠度计算常用方法有解析法和模拟法。

解析法中普遍采用的是一次二阶矩法,包括中心点法、验算点法(JC法)、映射变换法、实用分析法等。用Matlab实现这些方法的编程计算是相当方便快捷的,因为它有很多现成可用的计算程序,如求反函数、概率分布函数,概率密度函数,函数求导等等,而且采用矩阵运算可以避免循环结构。与同样方法的C语言实现,Matlab的实现语句要少的多。笔者从可靠度指标β的几何意义出发,建立求解可靠指标的优化模型,利用Matlab优化工具箱计算可靠度指标和验算点。

模拟法这里主要指蒙特卡罗法。蒙特卡罗法在目前可靠度计算中,被认为是一种相对精确法,它不受极限状态方程非线性的限制。对那些大型复杂结构,相应的极限状态功能函数往往是非线性的,如果采用解析法中求导的方法计算可靠度,当考虑因素较多时,极限状态功能函数将变得难以处理。蒙特卡罗法回避了可靠度分析中的数学困难,不需要考虑极限状态曲面的复杂性。此外,若基本变量相关,可利用条件概率密度,把多维问题化为一维问题来解决。因此从理论上来说,该方法的应用几乎没有什么限制。但在实际问题中,变量连续型分布是很复杂的。有的只能给分布函数的解析表达式,但给不出其反函数的解析表达式,如著名的β分布;有的则连分布函数的解析表达式都给不出。所以通常情况下对连续型分布采用直接抽样是有一定困难的。笔者利用Matlab的强大数值计算功能,实现了在Matlab中采用蒙特卡罗直接抽样计算结构可靠度,较好地解决了上述问题。将Matlab用于蒙特卡罗法的一个显著优点是它拥有功能强大的随机数发生器指令。

[1]

∑[(X

3

i

σXi]-μXi)/

2

开始时验算点未知,把β看成极限状态曲面上点P(X1,

X2,…,Xn)的函数,通过优化求解,找到β最小值,即可得可

333

靠度指标β和验算点P3(X1,X2,…,Xn)。求解可靠度指标

可以归结为以下约束优化模型:

β=min　

2

∑[(X

i=1

n

3

i

σXi]-μXi)/

2

333

S1t1　Z=g(X1,X2,…,Xn)=0

如果极限状态方程中某个变量(Xj)可用其他变量表示,即:

Xj=h(X1,X2,…,Xj-1,Xj+1,…,Xn),则上述模型可转化

为无约束优化模型:

minβ=

3

3

3

2

i=1,i≠j

∑[(X

n

3

i

33

σXi]2+{[h(X1-μ,X2,…,Xi)/

σXj}2Xj-1,Xj+1,…,Xn)-μXj]/

212　算　例

例:已知非线形极限状态方程Z=g(f,r,H)=567fr-2

015H=0,f服从正态分布,μ服从正f=016,σf=010786;r态分布,μr=2118,σr=010654;H服从对数正态分布,μH=

3218,σH=01984。f,r,H相互独立,求可靠度指标β及验算

[2]点(f3,r3,H3)。

μ解:先将H当量正态化:(变异系数为δH=σH/H=

0103)

[收稿日期]2006-05-25[作者简介]章慧健,硕士研究生。

154

　 计算

H服从对数正态分布,则h=lnH服从正态分布,且

机应用 　

μ′=H

μ+2

σH=3149,　′=2

ln(1+δH)=0103

他两个变量表示

CC=sqrt(((x(1)-016)/010786)^2+((x(2)-2118)/010654)^2+…((h-3149)/0103)^2);%β表达式

运行如下:

〉〉noncon

可靠度指标为bata=11964254

验算点为[01456167,21158971,331418877]

(1)采用约束优化法:

求解该问题可靠度指标可归结为如下约束优化数学模型:

23

σf]2+[(r3-μr)/σr]2+[(h3-minβ=[(f-μf)/

2

μ)/σH′′]H

s1t1　567fr-015HMatlab源程序如下:

3332

=567fr-015e

332h3

=0

3　在Matlab中实现蒙特卡罗法

311　基本原理

functionrelia　%定义主函数

x0=[016,2118,3149];%初始迭代点,这里取均值%调用优化工具箱求解

通常在用蒙特卡罗直接抽样法时,必须解X=F-1(r)以求得服从相应分布类型的随机变量X,再代入功能函数求

解。而Matlab(610版)提供了23种随机变量分布类型的随机数发生器(如正态分布、对数正态分布、泊松分布、威布尔分布等),可直接产生变量X,省去了可能会,极大地提高了效

2

率。r=norrnd(N(u,s)分布的m

);options=optimset(′LargeScale′,′off′

[x,fval]=fmincon(@obj,x0,[],[],[],[],[],[],@st,options);

fprintf(′可靠度指标为bata=%f\n′,fval);fprintf(′验算点为[%f,%f,%f]\n′,x(1),x(2),exp(x(3)));

functionCC=obj(x)　%CC=sqrt(((x(1)-1+(x)-2118)/0+…(x3)-3149)/0103)^

2);%即为β

function[c,ceq]=st(x)　%约束条件子函数c=[];

ceq=5673x(1)3x(2)-0153exp(x(3))^2;

r。目前版本尚无极值,故对极值型分布,则需先用均匀分

),产生一个m′r=rand(m′,n′行n′列的均布随机变量数组r,然后用式X=F-1(r)解出服从极值分布的随机变量数组。

设功能函数为g(X1,X2,…,Xn),Xi(i=1,2,…,n),服从各相应概率分布函数的随机变量。首先用相应随机数发生器指令产生m′×n′的随机变量数组Xi(k,l),(i=1,2,…,

)。然后,将各数组中各元素n;k=1,2,…,m′;l=1,2,…,n′

运行结果如下:

〉〉constr

可靠度指标为bata=11964254

验算点为[01456165,21158960,331418735](2)采用无约束优化法

求解该问题可靠度指标也可归纳为如下约束优化数学模型:

23

σf]2+[(r3-μr)/σr]2+[(h3-minβ=[(f-μf)/

2

μ)/σH′′]H

一一对应代入功能函数。最后得到功能函数结果数组,统计结果数组中小于等于0的元素个数j,可得失效概率:

Pf=

m′×n′

312　算　例

题目同212算例。

Matlab程序如下:

functionmengte(x)　%定义带参数主函数

n=x;　%设定抽样次数

f=normrnd(016,010786,1,n);%产生f变量的随机数r=normrnd(2118,010654,1,n);%产生r变量的随机数H=exp(normrnd(3149,0103,1,n);%产生H变量的随

33

h=ln(2×567×f×r)

2

3

将h3代入目标函数就是无约束优化问题了。

Matlab程序只需对约束优化程序稍作修改,约束条件子函数取消即可。

functionnoncon　%定义主函数

x0=[016,2118];%初始迭代点,这里取均值%调用优化工具箱求解

[x,fval]=fminsearch(@obj,x0);

h=0153log(11343x(1)3x(2));%计算验算点h值fprintf(′可靠度指标为bata=%f\n′,fval);fprintf(′验算点为[%f,%f,%f]\n′,x(1),x(2),exp(h));

functionCC=obj(x)　%目标函数子函数

h=0153log(11343x(1)3x(2));%随机变量h用其

机数

z=5673f13r-0153H1^2;%计算功能函数的值,注意

是矩阵计算

zz=find(z<0);%找出功能函数小于零的值,存放于zz中

k=length(zz);%求出功能函数小于零的个数

pf=k/n;%失效概率

fprintf(′可靠度指标bata=%f\n′,norminv(1-pf));

56分别进行104、10、10次抽样,计算结果如下:〉〉mengte(10000)

可靠度指标bata=11946456〉〉mengte(100000)

155

　 计算机应用 　

56

分别进行104、10、10次抽样,计算结果如下:〉〉mengte2(10000)

可靠度指标bata=11964694

验算点为[01455022,21160885,331391712]〉〉mengte2(100000)

可靠度指标bata=11956893

〉〉mengte2(1000000)

可靠度指标bata=11950253

常规的蒙特卡罗法只能计算出结构失效概率,并以此换算出可靠度指标,而无法得到验算点。而基于最优化原理的蒙特卡罗法,可同时得出结构失效概率、可靠度指标和验算点。改进程序如下:

functionmengte2(x)　%定义带参数主函数n=x;　%设定抽样次数f=normrnd(016,010786,1,n);%产生f变量的随机数r=normrnd(2118,010654,1,n);%产生r变量的随机数hh=11343(f13r);H=hh1^015;　%由极限状态方程计算变量H%随机变量f、r、H分别映射变换为y1、y2、y3fx1=normcdf(f,016,010786);y1=norminv(fx1,0,1);

fx2=normcdf(r,2118,010654);y2=norminv(fx2,0,1);

fx3=normcdf(log(H),3149,0103);y3=norminv(fx3,0,1);b1=sqrt(y11^2+y2^y3b=min(b1);c=find(b1==位置

fprintf(′可靠度指标bata=%f\n′,b);fprintf(′验算点为[%f,%f,%f]\n′,f(c),r(c),H(c));

(上接第153页)　强大的控制工具,如单元大小和形状的控

可靠度指标bata=11964308

验算点为[01455747,21159526,331407770]〉〉mengte(1000000)

可靠度指标bata=11964268

验算点为[01455950,21159056,331411588]

4　结束语

Matlab的强大功能为结构可靠度计算提供了便利,科研人员可迅速编出科学高效的计算程序,大大提高了效率。与其它编程语言相比,Matlab编程可以以比其它语言少得多的语句实现同样的功能。因此,深信语言在可靠度计算中的应用,。

参献

1].Matlab实

[].四川建筑科学研究,2004(2).

[2]　李清富.工程结构可靠性原理[M].郑州:黄河水利出版社,

1999.

[3]　冯晓波,杨桦.用MATLAB实现蒙特卡罗法计算结构可靠度

[J].中国农村水利水电,2002(8).

[4]　李志华.基于Matlab优化工具箱的工程结构可靠度计算[J].

四川建筑科学研究,2005(3).

制、网格的划分类型以及网格的清除和细化,还可对实体模

型图直接划分网格[8]。

FLAC-3D和ANSYS所采用的单元体形状大致相同。但其每一单元节点编制的规则和节点坐标,即单元数据,却有一定的差别。因此,根据这两种软件单元形状及其单元数据的关系编写FLAC-3D-ANSYS接口程序是建模方法的关键。FLAC-3D模型的自动生成主要由以下步骤组成:ANSYS模型的建立、ANSYS和FLAC-3D的数据转换以及

[8]

FLAC-3D调用生成的模型数据文件。

根据对FLAC-3D与ANSYS单元数据关系的分析,可利用计算机程序设计语言编写了FLAC-3D-ANSYS接口程序包。最后,该程序直接产生FLAC-3D计算所需的数据文件,通过FLAC-3D命令“call”调入由接口程序输出的数据文件,并加入边界条件、初始条件以及岩土体的力学参数,即可生成数值模型。

间、水平有限,一些相关的技术问题并没有详细阐明和解决。希望本文能够抛砖引玉,激发大家在这方面的研究热情。

参考文献

[1]　ItascaConsultingGroup,Inc.FLAC-3D(FastLagrangianAnalysis

ofContinualin3Dimensions)UserManuals,Version

2111Minneapolis,Minnesota,USA,200216.

[2]　刘波,韩彦辉.FLAC原理、实例与应用指南[M].北京:人民交

4　结束语

笔者在岩土工程研究工作中多次使用FLAC-3D软件,

在学习和使用的过程中也参考了大量的文献(以论文为主),鉴于国内介绍该程序的书籍较少,本文对上述纷繁芜杂的文献资料进行了总结与归纳,同时阐述了笔者使用FLAC-3D程序的一些感受,希望能为正在学习、使用FLAC-3D的工程技术人员、研究人员提供一些借鉴与参考。由于笔者时

通出版社,2005.

[3]　姜海燕,李拥军.两种边坡稳定性分析的方法优缺点的探讨

[J].中国科技信息,2006(2).[4]　邹栋,郑宏.快速拉格朗日法及其在边坡稳定性分析中的应用

[J].矿业研究与开发,2005,25(5).[5]　潘鹏飞,梁峥祥,李军才.边坡治理数值模拟的研究[J].矿业

工程,2005,3(4).[6]　王向东,文江泉.用FLAC-3D进行土质高边坡稳定性分析

[J].西华大学学报,2005,24(3).[7]　谢斌,钟敏,陈广平.基于AutoCAD和FLAC的边坡稳定性分析

[J].岩土工程技术,2005,19(2).[8]　廖秋林,曾钱帮,刘彤,等.基于ANSYS平台复杂地质体FLAC

-3D模型的自动生成[J].岩石力学与工程学报,2005,24

(6).

[9]　魏继红,吴继敏,孙少锐.FLAC-3D在边坡稳定性分析中的应

用[J].勘察科学技术,2005(2).[10]　丁秀美,黄润秋,刘光士.FLAC-3D前处理程序开发及其工

程应用[J].地质灾害与环境保护,2004,15(2).

156

本文来源：https://www.bwwdw.com/article/mwyi.html