威布尔分布可详见A three-dimensional numerical meso-approach to modeling time- independent deformation and fracturing of brittle rocks(10.1016/j.compgeo.2019.103274)一文,岩体参数计算公式如下:

对于某个具体的zone,u0为力学参数初值,ωi是0~1范围内的随机数,m由现场试验获得,本次取m=47,ui为满足威布尔分布的参数值。以上为三群群友x.口述,如有问题,请和他battle,然后告诉我正确的理解,谢谢。
2.代码
以下代码省略了应力、应变监测,相关内容可查阅彭文斌老师所著教材。
import random as randimport mathimport itasca as itit.command("python-reset-state off")# 岩体参数初值C0 = 1e3 #initial value of cohesionD0 = 1800.0 #initial value of densityF0 = 30.0 #initial value of frictionT0 = 1e10 #initial value of tensionE0 = 6e5 #initial value of Young's ModulesV0 = 0.35 #initial value of poisson ratio#威布尔分布参数值m = 47.0# 用于存储岩体参数listC=[]D=[]F=[]T=[]E=[]V=[]# 设置模型it.command("""model newmodel largestrain offzone import 'liexi.f3grid'zone face skinzone cmodel assign mohr-coulombzone property dilation 10.0zone face apply velocity (0,0,9e-5) range group 'bottom2'zone face apply velocity (0,0,-9e-5) range group 'top'"""# 计算岩体参数并存储于对应list中count = it.zone.count()for i in range(count):w = rand.uniform(0,1)C.append(C0*(-1*math.log(1-w))**(1/m))D.append(D0*(-1*math.log(1-w))**(1/m))F.append(F0*(-1*math.log(1-w))**(1/m))T.append(T0*(-1*math.log(1-w))**(1/m))E.append(E0*(-1*math.log(1-w))**(1/m))V.append(V0*(-1*math.log(1-w))**(1/m))# 赋参count = 0for z in it.zone.list():z.set_prop('cohesion',C[count])z.set_prop('density',D[count])z.set_prop('friction',F[count])z.set_prop('tension',T[count])z.set_prop('young',E[count])z.set_prop('poisson',V[count])count += 1# 求解it.command("model step 8000")



图3-2 模型塑性区
本文只保证代码能跑,不保证结果合理。