matlab 模特卡罗计算节点可靠度

clear; clc mu=[29.1;33.33;565;200000;100];sigma=[3.783;4.0;39.55;4000;30];%均值标准差向量 sLn=sqrt(log(1+(sigma(1:4)./mu(1:4)).^2));%对数分布的正态化还原 mLn=log(mu(1:4))-sLn.^2/2; aEv=sqrt(6)*sigma(5)/pi; uEv=-psi(1)*aEv-mu(5); level1=0;level2=0; level3=0;level4=0; for i=1:5 CUnconfFc=-lognrnd(mLn(1),sLn(1)); BUnconfFc=CUnconfFc; PCUnconfFc=-lognrnd(mLn(2),sLn(2)); PBUnconfFc=PCUnconfFc; CFy=lognrnd(mLn(3),sLn(3)); BFy=CFy; CEs=lognrnd(mLn(4),sLn(4)); BEs=CEs; D=-evrnd(uEv,aEv); fid=fopen('csj.tcl','wt'); fprintf(fid,'%s\n','model BasicBuilder -ndm 2 -ndf 3'); fprintf(fid,'%s','set fName "file');fprintf(fid,'%d',i);fprintf(fid,'%s\n','";'); fprintf(fid,'%s\n','source ProcMKPC2.tcl'); fprintf(fid,'%s\n','puts "a"'); fprintf(fid,'%s\n','source ProcRC2.tcl'); fprintf(fid,'%s\n','puts "a"'); fprintf(fid,'%s\n','set Strfactor 145; set Lenfactor [expr 1/25.4];'); fprintf(fid,'%s\n','set colY 600; set colZ 600;'); fprintf(fid,'%s\n','set bmY 950; set bmZ 600;'); fprintf(fid,'%s\n','set colCov 35; set bmCov1 30; set bmCov2 30; set bmCov $bmCov1;'); fprintf(fid,'%s\n','set JointWidth [expr $colY]; set JointHeight [expr $bmY]; set JointDepth $colZ ;'); fprintf(fid,'%s\n','set BeamLength1 500; set BeamLength2 1850; set BeamLength3 700; set BeamLength4 800; set ColumnLength1 300; set ColumnLength2 700; set ColumnLength3 600;'); fprintf(fid,'%s\n','set JointVolume [expr $JointWidth*$JointHeight*$JointDepth];'); fprintf(fid,'%s','set CUnconfFc ');fprintf(fid,'%d\n',CUnconfFc); fprintf(fid,'%s\n','set CUnconfEc -0.002;'); fprintf(fid,'%s\n','set CTSmispace 50; set CTSspace 200; set CTSlength 2120; set CTSFy 568.7; set CTSarea 50.24;'); fprintf(fid,'%s','set CFy ');fprintf(fid,'%d\n',CFy); % fprintf(fid,'%s','set CEs ');fprintf(fid,'%d\n',CEs); fprintf(fid,'%s\n','set CsHratio 0.004216; set CAs1 615.4; set CAs2 490.6'); fprintf(fid,'%s\n','procMKPC $CUnconfFc $CUnconfEc $colY $colZ $colCov $CTSmispace $CTSlength $CTSFy $CTSarea $Strfactor $Lenfactor'); fprintf(fid,'%s\n','set CUnconfFcu [lindex $concreteProp 2]; set CUnconfEcu [lindex $concreteProp 3];'); fprintf(fid,'%s\n','set CConfmiFc [lindex $concreteProp 4]; set CConfmiEc [lindex $concreteProp 5];'); fprintf(fid,'%s\n','set CConfmiFcu [lindex $concreteProp 6]; set CConfmiEcu [lindex $concreteProp 7];'); fprintf(fid,'%s\n','procMKPC $CUnconfFc $CUnconfEc $colY $colZ $colCov $CTSspace $CTSlength $CTSFy $CTSarea $Strfactor $Lenfactor'); fprintf(fid,'%s\n','set CUnconfFcu [lindex $concreteProp 2]; set CUnconfEcu [lindex $concreteProp 3];'); fprintf(fid,'%s\n','set CConfFc [lindex $concreteProp 4]; set CConfEc [lindex $concreteProp 5];'); fprintf(fid,'%s\n','set CConfFcu [lindex $concreteProp 6]; set CConfEcu [lindex $concreteProp 7];'); fprintf(fid,'%s','set BUnconfFc ');fprintf(fid,'%d\n',BUnconfFc); fprintf(fid,'%s\n','set BUnconfEc -0.002;'); fprintf(fid,'%s\n','set BTSmispace 50; set BTSspace 100; set BTSlength 2860; set BTSFy 568.7; set BTSarea 50.24;'); fprintf(fid,'%s','set BFy ');fprintf(fid,'%d\n',BFy); % fprintf(fid,'%s','set BEs ');fprintf(fid,'%d\n',BEs); fprintf(fid,'%s\n','set BAs1 490.6; set BAs2 380; set BAs3 201; set BsHratio 0.002322;'); fprintf(fid,'%s\n','procMKPC $BUnconfFc $BUnconfEc $bmY $bmZ $bmCov $BTSmispace $BTSlength $BTSFy $BTSarea $Strfactor $Lenfactor'); fprintf(fid,'%s\n','set BUnconfFcu [lindex $concreteProp 2]; set BUnconfEcu [lindex $concreteProp 3];'); fprintf(fid,'%s\n','set BConfmiFc [lindex $concreteProp 4]; set BConfmiEc [lindex $concreteProp 5];'); fprintf(fid,'%s\n','set BConfmiFcu [lindex $concreteProp 6]; set BConfmiEcu [lindex $concreteProp 7];'); fprintf(fid,'%s\n','procMKPC $BUnconfFc $BUnconfEc $bmY $bmZ $bmCov $BTSspace $BTSlength $BTSFy $BTSarea $Strfactor $Lenfactor'); fprintf(fid,'%s\n','set BUnconfFcu [lindex $concreteProp 2]; set BUnconfEcu [lindex $concreteProp 3];'); fprintf(fid,'%s\n','set BConfFc [lindex $concreteProp 4]; set BConfEc [lindex $concreteProp 5];'); fprintf(fid,'%s\n','set BConfFcu [lindex $concreteProp 6]; set BConfEcu [lindex $concreteProp 7];'); fprintf(fid,'%s\n','set bs_fc [expr -$BUnconfFc]; set bs_fs $BFy; set bs_es $BEs; set bs_fsu 711.7; set bs_dbar 25; set bs_esh [expr $BsHratio*$BEs];'); fprintf(fid,'%s\n','set bs_wid $colZ; set bs_dep $bmY;'); fprintf(fid,'%s\n','set bsT_nbars 4; set bsB_nbars 6;'); fprintf(fid,'%s\n','set bs_ljoint $colY;'); fprintf(fid,'%s\n','set cs_fc [expr -$CUnconfFc]; set cs_fs $CFy; set cs_es $CEs; set cs_fsu 711.7; set cs_dbar 28; set cs_esh [expr $CsHratio*$CEs];'); fprintf(fid,'%s\n','set cs_wid $colZ; set cs_dep $colY;'); fprintf(fid,'%s\n','set cs_nbars 5;'); fprintf(fid,'%s\n','set cs_ljoint $bmY;'); fprintf(fid,'%s\n','node 1 0.0 0.0'); fprintf(fid,'%s\n','node 2 0.0 $ColumnLength1'); fprintf(fid,'%s\n','node 3 0.0 [expr $ColumnLength1+$ColumnLength2]'); fprintf(fid,'%s\n','node 4 0.0 [expr $ColumnLength1+$ColumnLength2+$ColumnLength3]'); fprintf(fid,'%s\n','node 5 [expr -$BeamLength1-$BeamLength2-$BeamLength3-$BeamLength4-$JointWidth/2] [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight/2]'); fprintf(fid,'%s\n','node 6 [expr -$BeamLength2-$BeamLength3-$BeamLength4-$JointWidth/2] [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight/2]'); fprintf(fid,'%s\n','node 7 [expr -$BeamLength3-$BeamLength4-$JointWidth/2] [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight/2]'); fprintf(fid,'%s\n','node 8 [expr -$BeamLength4-$JointWidth/2] [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight/2]'); fprintf(fid,'%s\n','node 9 [expr -$JointWidth/2] [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight/2]'); fprintf(fid,'%s\n','node 10 [expr $JointWidth/2] [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight/2]'); fprintf(fid,'%s\n','node 11 [expr $JointWidth/2+$BeamLength4] [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight/2]'); fprintf(fid,'%s\n','node 12 [expr $JointWidth/2+$BeamLength4+$BeamLength3] [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight/2]'); fprintf(fid,'%s\n','node 13 [expr $JointWidth/2+$BeamLength4+$BeamLength3+$BeamLength2] [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight/2]'); fprintf(fid,'%s\n','node 14 [expr $JointWidth/2+$BeamLength4+$BeamLength3+$BeamLength2+$BeamLength1] [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight/2]'); fprintf(fid,'%s\n','node 15 0.0 [expr $ColumnLength1+$ColumnLength2+$ColumnLength3+$JointHeight]'); fprintf(fid,'%s\n','node 16 0.0 [expr $ColumnLength1+$ColumnLength2+2*$ColumnLength3+$JointHeight]'); fprintf(fid,'%s\n','node 17 0.0 [expr $ColumnLength1+2*$ColumnLength2+2*$ColumnLength3+$JointHeight]'); fprintf(fid,'%s\n','node 18 0.0 [expr 2*$ColumnLength1+2*$ColumnLength2+2*$ColumnLength3+$JointHeight]'); fprintf(fid,'%s\n','uniaxialMaterial Elastic 1 10000000000.0'); fprintf(fid,'%s\n','uniaxialMaterial Concrete01 10 $BUnconfFc $BUnconfEc $BUnconfFcu $BUnconfEcu'); fprintf(fid,'%s\n','uniaxialMaterial Concrete01 20 $BConfFc $BConfEc $BConfFcu $BConfEcu'); fprintf(fid,'%s\n','uniaxialMaterial Concrete01 30 $BConfmiFc $BConfmiEc $BConfmiFcu $BConfmiEcu'); fprintf(fid,'%s\n','uniaxialMaterial Steel02 40 $BFy $BEs $BsHratio 18.5 0.925 0.15 0.0 0.4 0.0 0.5'); fprintf(fid,'%s\n','uniaxialMaterial Concrete01 50 $CUnconfFc $CUnconfEc $CUnconfFcu $CUnconfEcu'); fprintf(fid,'%s\n','uniaxialMaterial Concrete01 60 $CConfFc $CConfEc $CConfFcu $CConfEcu'); fprintf(fid,'%s\n','uniaxialMaterial Concrete01 70 $CConfmiFc $CConfmiEc $CConfmiFcu $CConfmiEcu'); fprintf(fid,'%s\n','uniaxialMaterial Steel02 80 $CFy $CEs $CsHratio 18.5 0.925 0.15 0.0 0.4 0.0 0.5'); fprintf(fid,'%s\n','set z [expr $colZ/2.0]; set y [expr $colY/2.0];'); fprintf(fid,'%s\n','section Fiber 1 {'); fprintf(fid,'%s\n','patch rect 70 8 1 [expr $colCov-$y] [expr $colCov-$z] [expr $y-$colCov] [expr $z-$colCov]'); fprintf(fid,'%s\n','patch rect 50 2 1 [expr -$y] [expr $colCov-