厚壁圆筒问题径向、切向应力求解及Extra云图显示

# case-01import itasca as itimport mathit.command("python-reset-state false")# 模型参数r1 = 50 # 内半径/mmr2 = 100 # 外半径/mmth = 5 # 轴向厚度/mms_r = 30 # 径向网格数量s_th = 1 # 轴向网格数量s_t = 30 # 环向网格数量E = 2e5 # 弹性模量/MPav = 0.3 # 泊松比p = -10 # 内孔压力/MPa# 建模计算it.command("""           model new           model largestrain off           ;建模           zone create cylindrical-shell point 0(0,0,0) point 1({0},0,0) point 2(0,{1},0) point 3(0,0,{0}) point 4({0},{1},0)...                                         point 5(0,{1},{0}) point 8({2},0,0) point 9(0,0,{2}) point 10({2},{1},0) point 11(0,{1},{2}) size {3} {4} {5}           zone reflect origin (0,0,0) normal (-1,0,0)           zone reflect origin (0,0,0) normal (0,0,-1)           ;赋予本构、参数           zone cmodel assign elastic           zone property young {6} poisson {7}           ;边界条件           zone face skin           zone face apply velocity-normal 0 range group 'south' or 'north'           zone face apply stress-normal {8} range group 'west'           ;求解           model solve           """.format(r2,th,r1,s_r,s_th,s_t,E,v,p))# 圆心坐标xc = 0.0zc = 0.0# 遍历求解各单元径向、切向应力for z in it.zone.list():    ;计算夹角正弦、余弦    xz = z.pos_x()    zz = z.pos_z()    r = math.sqrt((xz-xc)**2+(zz-zc)**2)    sin = zz/r    cos = xz/r    ;获取应力分量    sxx = z.stress()[0][0]    szz = z.stress()[2][2]    sxz = z.stress()[2][0]    ;根据图2-1计算各单元径向、切向应力    sig_r = sxx*cos**2+szz*sin**2+2*sxz*cos*sin    sig_t = sxx*sin**2+szz*cos**2-2*sxz*cos*sin    ;设置径向、切向应力Extra云图    z.set_extra(1,sig_r)    z.set_extra(2,sig_t)    # 解析解sig_r_max = str(-1.0*r1**2*p/(r2**2-r1**2)*(1-r2**2/r1**2))sig_t_min = str(-1.0*r1**2*p/(r2**2-r1**2)*(1+r2**2/r2**2))sig_t_max = str(-1.0*r1**2*p/(r2**2-r1**2)*(1+r2**2/r1**2))print(sig_r_max,sig_t_min,sig_t_max)