Matlab soil mechanics concrete pile number stress calculation program

I have nothing to do this week A fellow brother is going to graduate By the way, I wrote him a pile Stress procedure of ,

According to the formula written in the paper he referred to I always feel the atmosphere of making graduation design ,

The paper is Nonlinear analysis of pile load−settlement behavior in layered soil Go and have a look if you are interested Just a bar

In other words, how to calculate the beautiful curve inside ...

It is divided into 3 Three main functions of soil layer , Add two iterations

At present, it's strange to make up the paste The last program plus diagram for reference .

Because the iteration in the paper is too annoying, one layer of pile has to be calculated 65W Time ,matlab We all know the computational efficiency of ... So it saves a lot

function [ U,P ] = diedaihanshu( u,Z,P1)
%UNTITLED3  A summary of this function is shown here 
%    Detailed description here 
% Basic parameters of the third stage 
format long g
a=1.018e-5;
b=5.683e-8;
fsu=55.525e3;
emtsu=0;
Ai=1.0465e7;
if u<=4.626e-3
    fu=u/(b+a*u);
else
    fu=emtsu*u+fsu;
end
fu;
X=[];
fudao= 1/(b + a*u) - (a*u)/(b + a*u)^2;
fsi=fudao*u+(fu-fudao*u);
X=[X;fu,fudao,fsi];
B=-X(:,2);
emt1=sqrt(-B/Ai);
emt2=-sqrt(-B/Ai);
P1=2*1.196*207e3*u/(1-(0.3)^2);
Ci=-fu+fudao*u;
c1i=exp(-emt1*(Z))*(emt2*u +P1/Ai+emt2*Ci/B)/(emt2-emt1);
c2i=exp(-emt2*(Z))*(emt1*u +P1/Ai+emt1*Ci/B)/(emt1-emt2);
U=c1i*exp(emt1*(Z-0.1))+c2i*exp(emt2*(Z-0.1))-Ci/B;
syms x
y=c1i*exp(emt1*x)+c2i*exp(emt2*x)-Ci/B;
diff(y);
p=subs(diff(y),x,Z-0.1);
P=-35e6*(pi/4)*(1.196)^2*eval(p);
end

The following iterative procedure

p=[];
u=[];
X=[];
Y=[];
for i=2e-4:2e-4:4e-3
    x=i;
    P1=2*1.196*207e3*x/(1-(0.3)^2);
    for j=46.7:-0.1:30.3
        Z=j;
        [U,P]=diedaihanshu1(x,Z,P1);
        u=[u;U];
        p=[p;P];
        x=U;
        P1=P;
    end
    X=[X;u(length(u))]
    Y=[Y;p(length(p))]
end
u,p,X,Y;
plot(X,Y,'o')

graphics Abscissa is U The displacement of The iteration of each section of The ordinate is the stress. It can be seen that the curve increases before it reaches a certain value, and then it is the fixed value . 

After my analysis, the reason is  

if u<=4.626e-3
    fu=u/(b+a*u);
else
    fu=emtsu*u+fsu;
end

The cyclic discrimination here leads to , I don't care about the specific theory at all

But at first glance, this curve seems to be the same thing , Smooth basically no mutation .

iteration 2

X=[];
Y=[];
u1=1e-6
p1=2*1.196*207e3*u1/(1-(0.3)^2);
for Z=46.7:-0.1:0
    if Z>30.3
        [U,P]=diedaihanshu( u1,Z,p1);
        X=[X,Z];
        Y=[P,Y];
        u1=U;
        p1=P;
    elseif Z<=10.0
        [U,P]=diedaihanshu3(u1,Z,p1);
        X=[X,Z];
        Y=[P,Y];
        u1=U;
        p1=P;
    else 
        [U,P]=diedaihanshu2(u1,Z,p1);
        X=[X,Z];
        Y=[P,Y];
        u1=U;
        p1=P;
    end
end

plot(X,Y,'o')

Produced a strange three segment line

I can only think that the initial data it gives is more toxic , The soil layer is relatively poor ？

In fact, I always think soil mechanics is very fucking , Is a bunch of empirical formulas really good .

Ah Tomorrow, the elder brother will take this broken achievement to see his mentor , Let's look forward to it .