【技术分享】Drucker-Prager 准则屈服面绘制

6天前浏览21

引言

在岩土数值计算中大本构模型是一个绕不开的问题。岩土材料用的最多的屈服准则是Mohr Coulomb和Drucker Prager准则。Mohr Coulomb虽然用的很广泛，但是屈服面在空间形状是一个六棱锥，意味着屈服函数对主应力的导数不连续。因此我们看到大量的论文基于Drucker Prager进行本构的开发研究。笔者也在这方面开展了一些研究，但是对于屈服面的绘制一直没有具体画过。正好感觉有几个小时空闲，尝试用Python简单绘制一个Drucker-Prager 准则的屈服面。未来不忙的时候再更新其他准则如Mohr Coulomb和Hoek Brown的屈服面绘制。可能存在错误，如果发现了请指正，提前感谢！

本构的形式

Drucker-Prager准则是一种用于描述材料屈服行为的模型，广泛应用于岩土工程和金属塑性分析。其表达式通常基于应力张量的不变量，考虑了材料的内摩擦角和黏聚力。以下是Drucker-Prager准则的常见表达式：

+ -k≤0

其中：

：第一应力不变量(平均正应力的三倍)，即

：第二偏应力不变量，表示剪应力强度，定义为：

{ + + }

或偏应力表示为 = ,其中为偏应力张量。

：与材料的内摩擦角相关的参数，表征材料对静水压力的敏感性。

k：与材料的粘聚力c相关的参数，表征材料的抗剪强度。

与k的表达式与具体应用场景有关，通常根据Mohr-Coulomb内摩擦角与凝聚力c确定。常见形式为：

这些参数假设Drucker-Prager圆锥与Mohr-Coulomb六边形在外圆周处相切(外接圆假设)。若采用内接圆或中间圆假设，与k的表达式会有所不同。

Haigh-Westergaard应力空间坐标转换

Haigh-Westergaard坐标(也称为应力不变量空间，有时也叫做Lode坐标)与主应力空间之间的坐标转换是岩土工程、材料力学和有限元分析中常用的方法，特别是在研究屈服准则(如 Drucker-Prager 或 Mohr-Coulomb)时。Haigh-Westergaard 空间采用柱坐标表示，使用三个应力不变量( ,s, )(或类似的表示)来描述应力状态，如下图。而主应力空间使用三个主应力( , , )，两者之间的转换关系如下：

图1 Haigh-Westergaard应力空间

Haigh-Westergaard应力空间坐标转换对应的代码如下：

def haigh_westergaard_to_cartesian(rou_t, s_t, theta_t):
    sig1 = rou_t / np.sqrt(3) + np.sqrt(2 / 3) * s_t * np.cos(theta_t)
    sig2 = rou_t / np.sqrt(3) + np.sqrt(2 / 3) * s_t * np.cos(theta_t - np.pi / 3 * 2)
    sig3 = rou_t / np.sqrt(3) + np.sqrt(2 / 3) * s_t * np.cos(theta_t + np.pi / 3 * 2)
    return sig1, sig2, sig3

Drucker-Prager 准则屈服面绘制

利用上述公式变换，Drucker-Prager准则在Haigh-Westergaard 应力空间下的表达式为：

+s- =0

我们可以在先计算出s为0时对应的，然后假定计算区域为[- ,10 ]，计算屈服面上的一系列s。假定的变化范围为0~360°，则可以得到屈服面上所有的点，进而采用matplotlib绘制屈服面。

图2 Drucker-Prager 屈服面

来源：数字孪生与工程计算

【仿真技术】使用Fish实现FLAC3D计算结果快速导入TECPLOT

在之前的一篇文章中，我们介绍了如何使用Python编程语言将FLAC3D计算结果快速导入TECPLOT软件。公众号链接如下：【仿真技术】使用Python实现FLAC3D计算结果快速导入TECPLOTFish语言简介今天，我们将继续探讨将FLAC3D计算结果快速导入TECPLOT软件，但这次我们将使用一种专为FLAC3D而设计的编程语言——Fish。Fish是一种脚本编程语言，它针对FLAC3D软件进行了优化，使得处理和操作FLAC3D模型和数据变得更加便捷。如果您之前已经阅读了我们的文章，您会发现本质上代码的功能是相同的，只是编程语言不同。编程语言：Fish适用软件版本：FLAC3D7.0和TECPLOT2017注：本文代码是在前人基础上进一步编写的使用方法本节基于以FLAC3D使用手册中 “Simple Slope in Hoek-Brown Material”为例，展示使用Fish编程语言来实现将FLAC3D计算结果导入TECPLOT的功能。根据官方提供的代码，模型计算结果如下：打开FLAC3D7.0_TO_Tecplot2017.dat并直接运行，获得“0502-KL1-1-XSGK+DZ.dat”文件。打开TECPLOT，点击File-Load Data，载入“0502-KL1-1-XSGK+DZ.dat”文件，即可实现将FLAC3D计算结果导入TECPLOT。代码详解代码逻辑与PYTHON篇的代码基本一致，因此只介绍关键步骤。第一步：统计单元数量、创建TECPLOT文件并写入头部信息首先，我们需要创建一个TECPLOT文件，并设置文件的基本信息。以下是使用Fish语言编写代码：fish define ini_mesh2tec global IO_READ = 0 global IO_WRITE = 1 global IO_FISH = 0 global IO_ASCII = 1 global N_RECORD = 10 ;;;;;;;;;;;;;;;;The figure count of each line global ZONE_NGP = zone.gp.num(zone.head) array buf(1) global tec_file = '0502-KL1-1-XSGK+DZ.dat';; Edit the tec_range to set plot range command model range create 'tec_range' endcommandend@ini_mesh2tec;; Get number of zones to plotfish define get_nzoneglobal p_zglobal n_zoneglobal plotitglobal numzones n_zone = 0 numzones = 0 p_z = zone.head loop while p_z # null numzones = numzones + 1 plot_test if plotit = 1 then n_zone = n_zone + 1 endif p_z = zone.next(p_z) endloopend@get_nzone;; Write Tecplot File Headfish define write_head global n_zone buf(1) = 'TITLE = "FLAC3D to Tecplot 10"\n' buf(1) = buf(1) + 'VARIABLES = "X(m)" \n"Y(m)" \n"Z(m)" \n' buf(1) = buf(1) + '"DISP(mm)" \n"XDISP(mm)" \n"YDISP(mm)" \n"ZDISP(mm)" \n' buf(1) = buf(1) + '"SIG1(MPa)" \n"SIG2(MPa)" \n"SIG3(MPa)" \n' buf(1) = buf(1) + '"SXX(MPa)" \n"SYY(MPa)" \n"SZZ(MPa)" \n' buf(1) = buf(1) + '"SXY(MPa)" \n"SYZ(MPa)" \n"SZX(MPa)" \n' buf(1) = buf(1) + 'ZONE T="GLOBAL" \n' buf(1) = buf(1) + ' N=' + string(gp.num) + ',' buf(1) = buf(1) + ' E=' + string(n_zone) + ',' buf(1) = buf(1) + ' ZONETYPE=FEBrick \n' buf(1) = buf(1) + ' DATAPACKING=BLOCK \n' buf(1) = buf(1) + ' VARLOCATION=([8-16]=CELLCENTERED) \n' buf(1) = buf(1) + ' DT=(SINGLE SINGLE SINGLE' buf(1) = buf(1) + ' SINGLE SINGLE SINGLE SINGLE' buf(1) = buf(1) + ' SINGLE SINGLE SINGLE' buf(1) = buf(1) + ' SINGLE SINGLE SINGLE' buf(1) = buf(1) + ' SINGLE SINGLE SINGLE )' status = file.write(buf,1)end 第二步：写入节点坐标信息FLAC3D模型的节点坐标信息是构成模型的基础数据之一。我们需要将这些节点信息（例如节点坐标和位移）写入TECPLOT文件。相关代码如下：fish define write_gp_info global p_gp global info_flag global gp_disp1 p_gp = gp.head loop while p_gp # null buf(1) = '' loop i(1,N_RECORD) if p_gp # null then caseof info_flag case 0 buf(1) = buf(1) + string(gp.pos.x(p_gp)) + ' ' case 1 buf(1) = buf(1) + string(gp.pos.y(p_gp)) + ' ' case 2 buf(1) = buf(1) + string(gp.pos.z(p_gp)) + ' ' case 4 get_gp_disp_Y buf(1) = buf(1) + string(gp_disp1*1000) + ' ' case 8 buf(1) = buf(1) + string(gp.disp.x(p_gp)*1000) + ' ' case 16 buf(1) = buf(1) + string(gp.disp.y(p_gp)*1000) + ' ' case 32 buf(1) = buf(1) + string(gp.disp.z(p_gp)*1000) + ' ' endcase p_gp = gp.next(p_gp) endif endloop status = file.write(buf,1) endloopend 第三步：写入单元信息FLAC3D模型由各种类型的单元组成。我们需要将这些单元的信息（例如应力和单元的连接性）写入TECPLOT文件。要注意的是代码中使用了zone.code（）函数，在FLAC3D中（0 = 砖形，1 = 楔形，2 = 金字塔形，3 = dbrick，4 = 四面体）。获取主应力的代码如下：fish define write_zone_info global p_z global info_flag global plotit p_z = zone.head loop while p_z # null buf(1) = '' loop i(1,N_RECORD) if p_z # null then plot_test if plotit = 1 then caseof info_flag case 0 buf(1) = buf(1) + string(zone.stress.prin.x(p_z)*0.000001) + ' ' case 1 buf(1) = buf(1) + string(zone.stress.prin.y(p_z)*0.000001) + ' ' case 2 buf(1) = buf(1) + string(zone.stress.prin.z(p_z)*0.000001) + ' ' case 4 buf(1) = buf(1) + string(zone.stress.xx(p_z)*0.000001) + ' ' case 8 buf(1) = buf(1) + string(zone.stress.yy(p_z)*0.000001) + ' ' case 16 buf(1) = buf(1) + string(zone.stress.zz(p_z)*0.000001) + ' ' case 32 buf(1) = buf(1) + string(zone.stress.xy(p_z)*0.000001) + ' ' case 64 buf(1) = buf(1) + string(zone.stress.yz(p_z)*0.000001) + ' ' case 128 buf(1) = buf(1) + string(zone.stress.xz(p_z)*0.000001) + ' ' endcase endif p_z = zone.next(p_z) endif endloop status = file.write(buf,1) endloopend 通过上述三步，我们能够将FLAC3D计算结果中的模型数据以Tecplot格式写入文件。这样一来，我们就可以方便地在TECPLOT软件中加载和分析FLAC3D的计算结果，进行进一步的可视化和后处理。END来源：数字孪生与工程计算