Geddes应力系数的matlab代码

zz from www.slope.com.cn

function main

%*********************参数说明
%L为桩长，r为计算点与桩轴线的水平距离，z为计算点距离地面深度。都取正值
%f1,f2,f3依次是点荷载、均布荷载以及三角形荷载下的Geddes应力解函数
%Ip,Ir,It依次是对应三个荷载下的应力系数，负表示拉，正表示压
%参考文献:J.D.Geddes Stresses in Foundation Soils Due to Vertical Subsurface Loading, 1966
%计算结果于Geddes的附加应力表是一致的
clc
clear
L=10;
r=2;
v=0.5;
z=3;

disp(’Ip, It, Ir’)

[Ip,Ir,It]=Geddes(L,r,v,z)

disp(’***************’)

function [Ip,Ir,It]=Geddes(L,r,v,z)

m=z/L;
n=r/L;
disp(’m and n’)
[m n]

if n==0
disp(’计算点位于桩轴线上，采用n=0.002进行修正，并模拟刺入沉降!’)
n=0.002;
end

Ip=f1(m,n,v);
Ir=f2(m,n,v);
It=f3(m,n,v);

function Ip_sub=f1(m,n,v) %点荷载

F=sqrt(m^2+n^2);
A=sqrt(n^2+(m-1)^2);
B=sqrt(n^2+(m+1)^2);
temp0=1/(8*pi*(1-v));
temp1=(1-2*v)*(m-1)/A^3-(1-2*v)*(m-1)/B^3+3*(m-1)^3/A^5;
temp2=(3*(3-4*v)*m*(m+1)^2-3*(m+1)*(5*m-1))/B^5;
temp3=30*m*(m+1)^3/B^7;
Ip_sub=-temp0*(temp1+temp2+temp3);

function Ir_sub=f2(m,n,v) %矩形荷载
F=sqrt(m^2+n^2);
A=sqrt(n^2+(m-1)^2);
B=sqrt(n^2+(m+1)^2);

temp0=1/(8*pi*(1-v))
temp1=-2*(2-v)/A;
temp2=(2*(2-v)+2*(1-2*v)*(1+m)*m/n^2)/B;
temp3=-2*(1-2*v)*m^2/n^2/F+n^2/A^3;
temp4=(4*m^2-4*(1+v)*(m/n)^2*m^2)/F^3;
temp5=(4*m*(1+v)*(m+1)*(m/n+1/n)^2-(4*m^2+n^2))/B^3;
temp6=6*m^2*((m^4-n^4)/n^2)/F^5;
temp7=6*m*(m*n^2-(m+1)^5/n^2)/B^5;
Ir_sub=temp0*(temp1+temp2+temp3+temp4+temp5+temp6+temp7);

function It_sub=f3(m,n,v) %三角形荷载
F=sqrt(m^2+n^2);
A=sqrt(n^2+(m-1)^2);
B=sqrt(n^2+(m+1)^2);
temp0=1/(4*pi*(1-v));
temp1=-2*(2-v)/A;
temp2=(2*(2-v)*(4*m+1)-2*(1-2*v)*(1+m)*m^2/n^2)/B;
temp3=(2*(1-2*v)*m^3/n^2-8*(2-v)*m)/F+(m*n^2+(m-1)^3)/A^3;
temp4=(4*v*m*n^2+4*m^3-15*n^2*m-2*(5+2*v)*(m/n)^2*(m+1)^3+(m+1)^3)/B^3;
temp5=(2*(7-2*v)*m*n^2-6*m^3+2*(5+2*v)*(m/n)^2*m^3)/F^3;
temp6=(6*m*n^2*(n^2-m^2)+12*(m/n)^2*(m+1)^5)/B^5;
temp7=-(12*(m/n)^2*m^5+6*m*n^2*(n^2-m^2))/F^5;
temp100=((A+m-1)/(F+m))*((B+m+1)/(F+m));
temp8=-2*(2-v)*log(temp100);
It_sub=temp0*(temp1+temp2+temp3+temp4+temp5+temp6+temp7+temp8);