I would like to make surface plot of my probability of failure curve?

6 views (last 30 days)
Below is the main file and my Function Code an the bottom, i wanna plot ptobability of failure with changing the muv (velocity) and mub (pier width) in the same and plot in surface plot. I attached the function code (scourdepth1) and main code file (POFPOF)
clear all;
mue=0.82; %mean of unbiased factor of study case
sige=0.189; %SD of unbiased factor of study case
muk1=0.9; %lower limit of K1 of study case
sigk1=1.1; %upper limit K1 of study case
muk2=1; %lower limit of K of study case
sigk2=1; %upper limit of K2 of study case
muk3=1.1; %lower limit of K3 of study case
sigk3=1.3; %upper limit of K3 of study case
mud50=0.01; %Mean of D50
sigd50=0.0195; %SD of D50
muv=2.48; %mean value of stream velocity
sigv=0.28; %SD of stream velocity COV=0.28 due to Time Variant
muy=4.042; %mean value of stream depth of study case
sigy=0.4042; %SD of stream depth of study case COV=0.1 due to Time Variant
%mub=1.067; %mean of pier width of study case
sigb=0.0534; %SD of pier width of study case
i=1;
for mub=0.5:0.1:5
[Ds] = scourdepth1(muv,sigv,muy,sigy,mub,sigb,muk1,sigk1,muk2,sigk2,muk3,sigk3,mue,sige,mud50,sigd50);
count=nnz((3-Ds(:))<0);
total=size(Ds,2);
b1(i)=mub;
p(i)=count/total;
i=i+1;
end
plot(b1,p,'red')
title('Probability of Failure Vs. Flow Velocity');
xlabel('Flow Velocity (m/sec)');
ylabel('Probability of Failure');
grid on
function [Ds] = scourdepth1(muv,sigv,muy,sigy,mub,sigb,muk1,sigk1,muk2,sigk2,muk3,sigk3,mue,sige,mud50,sigd50)
for i=1:10000
e(i)=normrnd(mue,sige);
k1(i)=unifrnd(muk1,sigk1); %correction factor for pier shape
k2(i)=unifrnd(muk2,sigk2); %correction factor of angle of attack
k3(i)=unifrnd(muk3,sigk3); %correction factor of bed configuration
D50(i)=lognrnd(mud50,sigd50); %Median Grain Size
v(i)=normrnd(muv,sigv); %flow velocity
y(i)=lognrnd(muy,sigy); %flow stream depth
b(i)=lognrnd(mub,sigb); %effective pier width
k4(i)=0.35*((b(i)/D50(i))^0.19); %Correction Factor for soil properties grain size factor
f(i)=v(i)./(((9.81.*y(i))^0.5)); %Froude number
Ds(i)=e(i).*2.*y(i).*k1(i).*k2(i).*k3(i).*k4(i).*((b(i)./y(i))^0.65).*(f(i)^0.43); %maximum Predicted scour depth
end
end

Accepted Answer

Torsten
Torsten on 9 Oct 2022
clear all;
mue=0.82; %mean of unbiased factor of study case
sige=0.189; %SD of unbiased factor of study case
muk1=0.9; %lower limit of K1 of study case
sigk1=1.1; %upper limit K1 of study case
muk2=1; %lower limit of K of study case
sigk2=1; %upper limit of K2 of study case
muk3=1.1; %lower limit of K3 of study case
sigk3=1.3; %upper limit of K3 of study case
mud50=0.01; %Mean of D50
sigd50=0.0195; %SD of D50
%muv=2.48; %mean value of stream velocity
sigv=0.28; %SD of stream velocity COV=0.28 due to Time Variant
muy=4.042; %mean value of stream depth of study case
sigy=0.4042; %SD of stream depth of study case COV=0.1 due to Time Variant
%mub=1.067; %mean of pier width of study case
sigb=0.0534; %SD of pier width of study case
j=1;
for muv = 0.5:0.5:5
i=1;
for mub=0.5:0.5:5
[Ds] = scourdepth1(muv,sigv,muy,sigy,mub,sigb,muk1,sigk1,muk2,sigk2,muk3,sigk3,mue,sige,mud50,sigd50);
count=nnz((3-Ds(:))<0);
total=size(Ds,2);
b1(i,j)=mub;
b2(i,j)=muv;
p(i,j)=count/total;
i=i+1;
end
j=j+1;
end
surf(b1,b2,p)
function [Ds] = scourdepth1(muv,sigv,muy,sigy,mub,sigb,muk1,sigk1,muk2,sigk2,muk3,sigk3,mue,sige,mud50,sigd50)
for i=1:10000
e(i)=normrnd(mue,sige);
k1(i)=unifrnd(muk1,sigk1); %correction factor for pier shape
k2(i)=unifrnd(muk2,sigk2); %correction factor of angle of attack
k3(i)=unifrnd(muk3,sigk3); %correction factor of bed configuration
D50(i)=lognrnd(mud50,sigd50); %Median Grain Size
v(i)=normrnd(muv,sigv); %flow velocity
y(i)=lognrnd(muy,sigy); %flow stream depth
b(i)=lognrnd(mub,sigb); %effective pier width
k4(i)=0.35*((b(i)/D50(i))^0.19); %Correction Factor for soil properties grain size factor
f(i)=v(i)./(((9.81.*y(i))^0.5)); %Froude number
Ds(i)=e(i).*2.*y(i).*k1(i).*k2(i).*k3(i).*k4(i).*((b(i)./y(i))^0.65).*(f(i)^0.43); %maximum Predicted scour depth
end
end

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!