利用matlab画出修正剑桥模型在三种三轴试验下的曲线图

1 写在前面

本人是在读研二研究生一枚，这篇文章的内容只是高等土力学课程作业，写下来只是希望能够帮助有需要的人，如果有帮助的话希望可以给个赞，不是计算机专业出身，确实代码水平有限，欢迎评论交友。

2 正文

2.1 理论知识

所要计算的三种试验包括，排水试验，不排水试验以及等P试验；

根据修正剑桥模型，塑性体应变和剪应变按照下面的公式进行计算：

修正剑桥模型的总塑性体应变和剪应变按照下式计算：

                                                  （1）

其中，

分别计算三种试验条件下的：

排水试验下：

                                                                 

等p试验下：

不排水试验下：

首先，总体积应变为0，可以根据式（1）解出

然后根据

计算

，

太多字母了，不打了

2.2 使用matlab进行编程

使用matlab中的ode23()函数求解上面三种试验条件下的三组常微分方程组。

2.2.1先编写函数，命名为：

排水试验：myfun1

等p试验：myfun2

不排水试验：myfun3

其中，t为偏应力q，y(1)为应力比，y(2)为平均应力p'，y(3)为塑性体应变，y(4)为塑性剪应变

程序中的函数方程没有化到最简，而是按照上面的公式直接算式求解的，可能看着比较繁琐，见谅

2.2.2使用ode23()函数求解

2.2.3画图

画图过程中，代码延长了一些线的长度，应为属性剪切应变算到后面就会趋于无穷大，不延长一些线图不好看

2.2.4代码

clc;clear;close all
[t3,r3]=ode23(@myfun3,[0:0.01:155.28],[0;196;0;0]);
[t1,r1]=ode23(@myfun1,[0:0.01:487.60],[0;196;0;0]);
[t2,r2]=ode23(@myfun2,[0:0.01:266.56],[0;196;0;0]);

subplot(2,2,2);
A2=plot(100*r2(:,3),r2(:,1));
xtickformat('percentage');hold on
A1=plot(100*r1(:,3),r1(:,1));
xtickformat('percentage');hold on
A3=plot(100*[r3(:,3);0.2],[r3(:,1);r3(end,1)]);
xtickformat('percentage');axis([0 20 0 1.5]);hold on
xlabel('塑性剪应变','color','k','fontsize',10);
ylabel('应力比','color','k','fontsize',10);
legend([A1 A2 A3],'排水试验','等P试验','不排水试验','Location','SouthEast','FontSize',5);

subplot(2,2,1);
y=@(x)1.36*x;
B2=plot(r2(:,2),t2);hold on
B1=plot(r1(:,2),t1);hold on
B3=plot(r3(:,2),t3);hold on
fplot(y);hold on
xlabel('p(kPa)','color','k','fontsize',10);axis([0 400 0 600])
ylabel('q(kPa)','color','k','fontsize',10);
legend([B1 B2 B3],'排水试验','等P试验','不排水试验','Location','NorthWest','FontSize',5);

subplot(2,2,4);
C2=plot(100*r2(:,3),100*r2(:,4));axis equal
hold on
C1=plot(100*r1(:,3),100*r1(:,4));
hold on
C3=plot(100*[r3(:,3);0.2],100*[r3(:,4);r3(end,4)]);
xtickformat('percentage');ytickformat('percentage');axis([0 20 0 10]);hold on
xlabel('塑性剪应变','color','k','fontsize',10);
ylabel('塑性体应变','color','k','fontsize',10);
legend([C1 C2 C3],'排水试验','等P试验','不排水试验','Location','NorthWest','FontSize',5);
function z=myfun3(t3,y)
    z=[(y(2)+t3*(0.0665*2*y(1)/(1.68*y(2)*(1.36^2+y(1)^2)))/((0.0188+0.0665*(1.8496-y(1)^2)/(1.8496+y(1)^2))/(1.68*y(2))))/y(2)^2;...
        -(0.0665*2*y(1)/(1.68*y(2)*(1.36^2+y(1)^2)))/((0.0188+0.0665*(1.8496-y(1)^2)/(1.8496+y(1)^2))/(1.68*y(2)));...
        (0.0665*4*y(1)^2)/((1.68*y(2))*(1.36^4-y(1)^4))-(0.0665*2*y(1)/(1.68*y(2)*(1.36^2+y(1)^2)))^2/((0.0188+0.0665*(1.8496-y(1)^2)/(1.8496+y(1)^2))/(1.68*y(2)));...
        0];
end
function z=myfun1(t,y)
    z=[(y(2)-t/3)/(y(2)^2);...
        1/3;...
        (0.0665*2*y(1)/(1.68*y(2)*(1.36^2+y(1)^2)))/3+(0.0665*4*y(1)^2)/((1.68*y(2))*(1.36^4-y(1)^4));...
        ((0.0188+0.0665*(1.8496-y(1)^2)/(1.8496+y(1)^2))/(1.68*y(2)))/3+(0.0665*2*y(1)/(1.68*y(2)*(1.36^2+y(1)^2)))];
end
function z=myfun2(t,y)
    z=[1/y(2);...
        0;...
        (0.0665*2*y(1)/(1.65*y(2)*(1.36^2+y(1)^2)))*0/3+(0.0665*4*y(1)^2)/((1.68*y(2))*(1.36^4-y(1)^4));...
        ((0.0188+0.0665*(1.8496-y(1)^2)/(1.8496+y(1)^2))/(1.68*y(2)))*0/3+(0.0665*2*y(1)/(1.68*y(2)*(1.36^2+y(1)^2)))];
end

2.2.5结果

3写在最后

感谢阅读，希望文章对你有所帮助，只是简单的编程和推导，欢迎有更好相关程序的大佬与我交流。