基本信息

---

下载

   友情提示：（无需注册或充值，赞助后即可获取资源下载链接）

---

     嘿，亲！知识可是无价之宝呢，但咱这精心整理的资料也耗费了不少心血呀。小小地破费一下，绝对物超所值哦！如有下载和支付问题，请联系我们QQ(微信同号)：813200300

本次赞助数额为： 2 元　

支付宝支付

微信支付

   源码介绍

---

abaqus的二维（2D）的扩展有限元（XFEM）算J积分和应力强度因子（SIF）。从原理上想，应该不难实现，于是乎自己想写一个简单的小程序

import numpy as np def J_integral(ue,ve,xe,ye,qe,xcenter,ycenter,E,mu):
   ue=np.matrix(ue); # element 1 direction displacement (1x4)  ve=np.matrix(ve); # element 2 direction displacement (1x4)  xe=np.matrix(xe); # element 1 direction coordinates  (1x4)  ye=np.matrix(ye); # element 2 direction coordinates  (1x4)  qe=np.matrix(qe); # element q function value at element node (1x4)  g3 = 1/np.sqrt(3) 
   rs=np.matrix([[-g3,g3,g3,-g3],[-g3,-g3,g3,g3]]) # integral point's position   kapa = (3-mu)/(1 mu) # material constants  nu = E/2/(1 mu) #  const = E/(1-mu**2)
   Ce = const*np.matrix([[1,mu,0],[mu,1,0],[0,0,(1-mu)/2]]) # elastic matrix  Je = 0.0  Je2 = 0.0  Je4 = 0.0  for ip in range(0,4): # loop over 4 integral point  r = rs[0,ip]
      s = rs[1,ip]
      N = 0.25*np.matrix([(1-r)*(1-s),(1 r)*(1-s),(1 r)*(1 s),(1-r)*(1 s)]) # element shape function  Nr = 0.25*np.matrix([s-1,1-s,1 s,-1-s]) # the 1-direction derivative of shape function  Ns = 0.25*np.matrix([r-1,-1-r,1 r,1-r]) # the 2-direction derivative of shape funciton  xr = Nr*xe.H;xs=Ns*xe.H;yr=Nr*ye.H;ys=Ns*ye.H; # d(x,y)/d(ksi,eta)  detJe = xr*ys-xs*yr # element Jaucobi det at integral point  Nx = (ys*Nr-yr*Ns)/detJe # dN/dx  Ny = (-xs*Nr xr*Ns)/detJe # dN/dy  ux = Nx*ue.H # du/dx  uy = Ny*ue.H # du/dy  vx = Nx*ve.H # dv/dx  vy = Ny*ve.H # dv/dy  qx = Nx*qe.H # dq/dx  qy = Ny*qe.H # dq/dy  exx = ux[0,0] # exx  eyy = vy[0,0]
      exy = uy[0,0] vx[0,0]
      ee1 = np.matrix([exx,eyy,exy])
      ss1 = Ce*ee1.H
      ss1 = ss1.H
      xq = N*xe.H;yq=N*ye.H # integral point coordinates  rr = np.sqrt((xq- xcenter)**2 (yq- ycenter )**2)[0,0]
      theta = np.arctan2((yq- ycenter),(xq- xcenter))[0,0]
      ss2 = 1/np.sqrt(2*np.pi*rr)*np.matrix([1-np.sin(theta/2)*np.sin(theta*3/2),1 np.sin(theta/2)*np.sin(theta*3/2),np.sin(theta/2)*np.cos(theta*3/2)])*np.cos(theta/2)
      ux2 = 1/np.sqrt(2*np.pi*rr)/4/nu*\
      np.matrix([(8*np.cos(theta/2)**4-10*np.cos(theta/2)**2 kapa 1)*np.cos(theta/2),(8*np.cos(theta/2)**4-6*np.cos(theta/2)**2-kapa-1)*np.sin(theta/2)])
      ss3 = 1/np.sqrt(2*np.pi*rr)*np.matrix([-np.sin(theta/2)*(2 np.cos(theta/2)*np.cos(theta*3/2)),np.sin(theta/2)*np.cos(theta/2)*np.cos(theta*3/2),np.cos(theta/2)*(1-np.sin(theta/2)*np.sin(theta*3/2))])
      ux3 = 1/np.sqrt(2*np.pi*rr)/4/nu*\
      np.matrix([-(8*np.cos(theta/2)**4-6*np.cos(theta/2)**2 kapa 1)*np.sin(theta/2),(8*np.cos(theta/2)**4-10*np.cos(theta/2)**2-kapa 3)*np.cos(theta/2)])
      w = ss2*ee1.H
      Je=Je ((ss1[0,0]*ux2[0,0] ss1[0,2]*ux2[0,1] ss2[0,0]*ux[0,0] ss2[0,2]*vx[0,0]-w[0,0])*qx[0,0]   (ss1[0,1]*ux2[0,1] ss1[0,2]*ux2[0,0] ss2[0,1]*vx[0,0] ss2[0,2]*ux[0,0])*qy[0,0])*detJe
      w2 = ss3*ee1.H
      Je2=Je2 ((ss1[0,0]*ux3[0,0] ss1[0,2]*ux3[0,1] ss3[0,0]*ux[0,0] ss3[0,2]*vx[0,0]-w2[0,0])*qx[0,0]   (ss1[0,1]*ux3[0,1] ss1[0,2]*ux3[0,0] ss3[0,1]*vx[0,0] ss3[0,2]*ux[0,0])*qy[0,0])*detJe
      
      w4 = 0.5*ss1*ee1.H
      Je4 = Je4   ((ss1[0,0]*ux ss1[0,2]*vx-w4)*qx (ss1[0,2]*ux ss1[0,1]*vx)*qy)*detJe return [Je[0,0],Je2[0,0],Je4[0,0]] ## 前半段为参数设定及理论过程过程最后得到结果，后半段为数据提取引用 from odbAccess import * def post(Jr,xcenter,ycenter,problemtype,enter,fileName,centernode):
   odb = openOdb(fileName) try: if not enter:
         xcenter = centernode.coordinates[0]
         ycenter = centernode.coordinates[1] except: pass   step1 = odb.steps.values()[0]
   lastFrame = step1.frames[-1]
   displacement = lastFrame.fieldOutputs['U'].values
   totalnode=odb.rootAssembly.instances['PART-1-1'].nodes
   totalele = odb.rootAssembly.instances['PART-1-1'].elements
   mat = odb.materials
   E = mat['MATERIAL-1'].elastic.table[0][0]
   mu = mat['MATERIAL-1'].elastic.table[0][1] if problemtype=='plane strain':
      E = E/(1-mu**2)
      mu = mu/(1-mu)
   
   xe = [0,0,0,0]
   ye = [0,0,0,0]
   dis = [0,0,0,0]
   ue = [0,0,0,0]
   ve = [0,0,0,0]
   qe = [0,0,0,0]
   J1 = 0  J2 = 0  J4 = 0  eleee=[] for ele in totalele: for i in range(0,4):
         xe[i] = totalnode[ ele.connectivity[i]-1 ].coordinates[0]
         ye[i] = totalnode[ ele.connectivity[i]-1 ].coordinates[1]
         dis[i] = np.sqrt((xe[i]- xcenter)**2 (ye[i]- ycenter)**2)-Jr pass  if dis[0]*dis[1]<0 or dis[0]*dis[2]<0 or dis[0]*dis[3]<0: for i in range(0,4): for v in displacement: if v.nodeLabel==ele.connectivity[i]:
                  ue[i] = v.data[0]
                  ve[i] = v.data[1]
                  qe[i] = 0.5*(1-np.sign(dis[i]))
         [Je,Je2,Je4] = J_integral(ue,ve,xe,ye,qe,xcenter,ycenter,E,mu)
         J1 = J1   Je
         J2 = J2   Je2
         J4 = J4   Je4
         eleee.append(ele.label) print('Integral radius : %f'%Jr) print('KI = %f , KII = %f ,  J = %f '%(J1*E/2,J2*E/2,J4)) print(eleee) def main_Jintegral(Jr,xcenter,ycenter,problemtype,enter,fileName,centernode=None):
   post(Jr,xcenter,ycenter,problemtype,enter,fileName,centernode)