首页/文章/ 详情

Drucker-Prager准则：基础理论、MATLAB代码、优点与不足

现代石油人

7月前浏览652

符号列表

：Drucker-Prager材料常数；

：应力偏张量第二不变量；

：有效应力张量第一不变量；

：第一有效主应力；

：第二有效主应力；

：第三有效主应力；

：剪应力（octahedral:八面体的）；

：有效正应力（octahedral:八面体的）；

：单轴抗压强度；

：单轴抗拉强度；

：洛德角；

： (Mises– Schleicher and Drucker–Prager unified）参数，该参数定义了Drucker-Parger准则中平面的形状；

： (Mises– Schleicher and Drucker–Prager unified）参数;

：内摩擦角；

：内聚力；

1. 描述

Drucker-Prager破坏准则是一个三维的考虑压力的数学模型，用于估算岩石达到极限强度时的应力状态。该准则基于以下假设：破坏时的剪应力（三维情况下为八面体形状）通过材料常数与法向应力（三维情况下为八面体形状）成线性关系。

2. 背景

Drucker-Prager破坏准则被视作土体材料的广义的Mohr-Coulomb破坏准则，其数学表达为：(1)

其中：

(2)

(3)

因此Drucker-Prager准则认为是的函数。

3. 数学表达

参数和可以通过三轴测试获得内摩擦角和内聚力的大小得到：

(4)

4. MATLAB源代码（Drucker-Prager准则）

4.1 Drucker-Prager参数和数学表达

function f1 = fail(fai,I1,R,theta,i,j,m)

   %%Drucker-Prager参数和数学表达
  
  
    rad23=sqrt(2/3);
    
    T(1)=I1(m)/3+rad23*R(i)*sin(theta(j)+2/3*pi);
    T(2)=I1(m)/3+rad23*R(i)*sin(theta(j));
    T(3)=I1(m)/3+rad23*R(i)*sin(theta(j)-2/3*pi);
    
     T1=max(T);
     T3=min(T);
    
    trT=T(1)+T(2)+T(3);  
    Ts=T-1/3*trT;
    
    I1= T(1)+T(2)+T(3);
    I2=T(1)*T(2)+T(2)*T(3)+T(3)*T(1);
    I3=T(1)*T(2)*T(3)       ;   
    J2=-(Ts(1)*Ts(2)+Ts(2)*Ts(3)+Ts(3)*Ts(1));
    
    f1= sqrt(J2) - 2*sind(fai)/(sqrt(3)*(3-sind(fai)))*I1;
      
end

4.2 Drucker-Prager破坏准则主代码

%% Drucker-Prager破坏准则主代码

clear all;
clc;

fai=20;
     
for m=1:10:101;
    I1(m)=m;
%     I1=10;
    i=1;
    R(i)=0;
%Increase radius to get failure 
   for j=1:1:62;
   theta(j,m)=-pi/6+pi/30*(j-1);
       while fail(fai,I1,R,theta,i,j,m)<0
            i = i+1;                
            R(i)=R(i-1)+0.01;
        end
        Rf(j,m)=R(i);
   
        sig1(j,m)=I1(m)/3+sqrt(2/3)*Rf(j,m)*sin(theta(j)+2/3*pi);
        sig2(j,m)=I1(m)/3+sqrt(2/3)*Rf(j,m)*sin(theta(j));
        sig3(j,m)=I1(m)/3+sqrt(2/3)*Rf(j,m)*sin(theta(j)-2/3*pi);   

       i=1;
       R(i)=0;
   end
end   

   
  
%%
% close all
% subplot(1,2,1);

figure(1)
mesh(sig1,sig2,sig3,'EdgeColor','k'); hold on
view([20,5]);
axis on;
axis equal;
set(gca,'Fontsize',15,'FontName','Times new Roman');
set(get(gca,'XLabel'),'FontSize',18,'FontName','Times new Roman');
set(get(gca,'YLabel'),'FontSize',18,'FontName','Times new Roman');

% subplot(1,2,2);
% polar(theta,Rf,'b');
% axis equal;