在好例子网,分享、交流、成长!
您当前所在位置:首页Python 开发实例Python语言基础 → 二维裂纹扩展J积分以及应力强度因子计算

二维裂纹扩展J积分以及应力强度因子计算

Python语言基础

下载此实例
  • 开发语言:Python
  • 实例大小:4.62KB
  • 下载次数:5
  • 浏览次数:42
  • 发布时间:2024-11-23
  • 实例类别:Python语言基础
  • 发 布 人:游大侠
  • 文件格式:.py
  • 所需积分:2
 相关标签: 二维 扩展 积分 计算

实例介绍

【实例简介】abaqus的二维(2D)的扩展有限元(XFEM)算J积分和应力强度因子(SIF)。从原理上想,应该不难实现,于是乎自己想写一个简单的小程序
【实例截图】from clipboard
【核心代码】
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)

实例下载地址

二维裂纹扩展J积分以及应力强度因子计算

不能下载?内容有错? 点击这里报错 + 投诉 + 提问

好例子网口号:伸出你的我的手 — 分享

网友评论

发表评论

(您的评论需要经过审核才能显示)

查看所有0条评论>>

小贴士

感谢您为本站写下的评论,您的评论对其它用户来说具有重要的参考价值,所以请认真填写。

  • 类似“顶”、“沙发”之类没有营养的文字,对勤劳贡献的楼主来说是令人沮丧的反馈信息。
  • 相信您也不想看到一排文字/表情墙,所以请不要反馈意义不大的重复字符,也请尽量不要纯表情的回复。
  • 提问之前请再仔细看一遍楼主的说明,或许是您遗漏了。
  • 请勿到处挖坑绊人、招贴广告。既占空间让人厌烦,又没人会搭理,于人于己都无利。

关于好例子网

本站旨在为广大IT学习爱好者提供一个非营利性互相学习交流分享平台。本站所有资源都可以被免费获取学习研究。本站资源来自网友分享,对搜索内容的合法性不具有预见性、识别性、控制性,仅供学习研究,请务必在下载后24小时内给予删除,不得用于其他任何用途,否则后果自负。基于互联网的特殊性,平台无法对用户传输的作品、信息、内容的权属或合法性、安全性、合规性、真实性、科学性、完整权、有效性等进行实质审查;无论平台是否已进行审查,用户均应自行承担因其传输的作品、信息、内容而可能或已经产生的侵权或权属纠纷等法律责任。本站所有资源不代表本站的观点或立场,基于网友分享,根据中国法律《信息网络传播权保护条例》第二十二与二十三条之规定,若资源存在侵权或相关问题请联系本站客服人员,点此联系我们。关于更多版权及免责申明参见 版权及免责申明

;
报警