-
Notifications
You must be signed in to change notification settings - Fork 3
/
MCC_main.m
140 lines (108 loc) · 3.11 KB
/
MCC_main.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
%% 5R5 EP2 Q4 and Q6
% Modified (Burland's) Cam Clay
% Strategy: model everything in effective stress
% To impose undrained shear, implose e2 = e3 = -0.5 e1, gammas = 0
close all;
clearvars;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input these values before running Initialise_MCC
% Material properties
kappa = 0.05;
lambda = 0.20;
nu = 0.20; % Only used in elastic
phics = 33.5 * pi() / 180; % Critical state friction angle
patm = 100;
vatm = 2.12;
% State variables: Initial conditions
vstart = vatm; % specific volume at start, dimensionless 2.12
pstart = patm;
pc = pstart; % kPa
s = [100 99.9 99.9 0 0 0]'; % kPa. Note dqds will blow up if s1==s2==s3 at start
% Targets
pcompact = 400;
pswell = 100;
% Iteration setting
destep = 1e-5; % Absolute strain, dimensionless. Strain step.
% End of things to set before running Initialise_MCC
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Initialise_MCC;
Update_State_MCC;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do whatever you want here with the soil
% The Save_State routine also creates the storage variables
% if they don't already exist
Impose_Strain_MCC;
Update_State_MCC;
%
% Isotropically compress
while (p < pcompact)
% Try to impose some strain increment
de = destep .* [1 1 1 0 0 0]';
Impose_Strain_MCC;
Update_State_MCC;
pt = p;
u = 0;
Save_State;
end
ti %Print ti
% Isotropically swell
while (p > pswell)
% Try to impose some strain increment
de = -destep .* [1 1 1 0 0 0]';
Impose_Strain_MCC;
Update_State_MCC;
pt = p;
u = 0;
Save_State;
end
tiStartShear = ti
%% Undrained shear to failure
for eq = 0:destep:0.2
de = destep .* [1 -0.5 -0.5 0 0 0]';
Impose_Strain_MCC;
Update_State_MCC;
pt = pswell + q / 3; % Total stress
u = pt - p;
Save_State;
end
%%
% Plot the results
figure(901);
semilogx(ResultsStruct.p, ResultsStruct.v);
xlabel('p (kPa, log scale)');
ylabel('v');
figure(902);
plot(ResultsStruct.pc);
title('pc');
figure(903);
hold all;
plot(ResultsStruct.p, ResultsStruct.q);
plot(ResultsStruct.pt, ResultsStruct.q);
xlabel('p (kPa)');
ylabel('q (kPa)');
legend('Effective stress', 'Total stress');
strcat('Peak excess pore pressure: ', num2str(max(ResultsStruct.pt - ResultsStruct.p)))
strcat('Ultimate excess pore pressure: ', num2str(ResultsStruct.pt(end) - ResultsStruct.p(end)))
strcat('Peak shear strength: ', num2str(max(ResultsStruct.q)/2))
strcat('Residual shear strength: ', num2str(ResultsStruct.q(end)/2))
figure(913);
plot(ResultsStruct.e(:,1) - ResultsStruct.e(tiStartShear,1), ResultsStruct.u);
xlabel('Axial strain');
ylabel('Excess pore pressure (kPa)');
figure(904);
plot(ResultsStruct.e(:,1) - ResultsStruct.e(tiStartShear,1), ResultsStruct.q);
xlabel('Axial strain');
ylabel('q (kPa)');
figure(905);
plot(ResultsStruct.s(:,2));
title('s2');
figure(906);
hold on;
plot(ResultsStruct.F);
plot(ResultsStruct.ndF);
legend('F', 'ndF');
figure(907);
hold on;
plot(ResultsStruct.v);
plot(ResultsStruct.vCSL);
legend('v', 'vCSL');