Formulating a Monte carlo simulation for corrosion initiation time

Hi, I am currently working on formulating a Monte Carlo simulation for predicting corrosion initiation time whereby the formula being applied is

T = X^2 / [4D(erfc^-1(Cth/Co)^2]

where X, D, Cth, Co are random variables with lognormal distributions for X,Co and Cth is a uniform distribution and D a distinct value.

Below is some of the code for the variables required:

X = normrnd(40,0.1,1,1000);
Co = lognrnd(3.5,0.5,1,1000);
D = 5*(10^-5);
Cth= unifrnd(0.6,1.2,1,1000);
E= erfc(Cth/Co);

Where the issue arises is when I try to implement the equation as the random variables are vectors so the power function becomes inapplicable as it requires scalar quantities. More information and assistance is kindly welcome.

Answer 1

[Image Analyst]

I don't know about the erfc -1 versus the erfc, but how about this:

X = normrnd(40, 0.1, 1, 1000);
Co = lognrnd(3.5, 0.5, 1, 1000);
D = 5*(10^-5);
Cth = unifrnd(0.6, 1.2, 1, 1000);
E = erfc(Cth ./ Co);

T = X .^ 2 ./ [4 * D .* E .^ 2] 
plot(T, 'b-');
grid on;

Simply use dot before any * / or ^ operations that use a vector.

[Barry Akenga]

Thank you for your response. To provide further details on what I want to achieve is noted on the plot on the image whereby the result is a probability density function that predicts corrosion initiation time. The error function erfc-1 is applied in the equation however as mentioned earlier the random variables are X,Co,Cr whereas Dc or D is an assumed constant.

[Barry Akenga]

The plot arrived at is a Lognormal distribution plot for the probability density function. Which is also the plot I want to arrive at for the pdf.

[Barry Akenga]

Also I think I have made a mistake whereby, the random variables generation does not have a sigma(std) but has a COV(coefficient of variation). How are you able to generate the relative random variables for the code on matlab?

Answer 2

[Pablo Benítez]

Dear Barry,

I hope you are well. Interestingly, I am working on a similar research to yours, however, the time of corrosion initiation that I am looking for is for carbonation. Have you achieved to get the PDF for the initiation time?

I will be grateful to be in contact with you,

Regards,

Pablo

[Barry Akenga]

Hi Pablo,

Apologies for the late response, I only managed to stumble across your query while perusing through outstanding mail on my inactive account. I did manage to conduct this research and obtain data through matlab. Is there anything in particular I may help you out on?

Regards, Barry

[harsh Brar]

hi barry,

Hope you are doing well.

Lately I have been working on predicting the service life of concrete on the basis of chloride induced corrosion using monte carlo. I have been trying my best to get results from the code but its not working. I was hoping if you could help me with your code. That will be of great help.

looking forward to hear from you soon.

thank you.

Answer 3

[Pablo Benítez]

Hi Barry,

I hope you are going well. Fortunately, I managed to solve the problem of obtaining the time of corrosion onset by carbonation. I have had some problems to perform the Monte Carlo simulation at the beginning, but the problem has been fixed. Anyway, thanks for your response.

Best regards,

Pablo

Answer 4

[Barry Akenga]

Hi Pablo,

Thanks for getting back to me on this, are you able to convey the code you used to obtain the results required for your Monte Carlo simulation. I did this during my final year for my thesis, however I went a step further and refined my data obtained from the Monte Carlo simulation to obtain variable specific data. Thanks in advance.

Regards, Barry.

Answer 5

[Pablo Benítez]

Hi Barry,

I will need to talk to the colleague who helped me with the programming. The code for MCS comprises a part of a larger code that was made for the maintenance planning of structures. Thus, the complete code is currently being reorganized to make a commercial package that will be used for decision-making. So, I do not know if there is any confidentiality to convey a part of the code. If there is no problem, I will send you.

Cheers,

Pablo

Answer 6

[Abdelrahman Abdallah]

hi barry,

I am working on the same equation and I just wanted to know

did you use the erfinv in matlab and why didnt you use lognrnd(log(3.5),0.5) for surface cholride