实例介绍
【实例简介】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)
好例子网口号:伸出你的我的手 — 分享!
小贴士
感谢您为本站写下的评论,您的评论对其它用户来说具有重要的参考价值,所以请认真填写。
- 类似“顶”、“沙发”之类没有营养的文字,对勤劳贡献的楼主来说是令人沮丧的反馈信息。
- 相信您也不想看到一排文字/表情墙,所以请不要反馈意义不大的重复字符,也请尽量不要纯表情的回复。
- 提问之前请再仔细看一遍楼主的说明,或许是您遗漏了。
- 请勿到处挖坑绊人、招贴广告。既占空间让人厌烦,又没人会搭理,于人于己都无利。
关于好例子网
本站旨在为广大IT学习爱好者提供一个非营利性互相学习交流分享平台。本站所有资源都可以被免费获取学习研究。本站资源来自网友分享,对搜索内容的合法性不具有预见性、识别性、控制性,仅供学习研究,请务必在下载后24小时内给予删除,不得用于其他任何用途,否则后果自负。基于互联网的特殊性,平台无法对用户传输的作品、信息、内容的权属或合法性、安全性、合规性、真实性、科学性、完整权、有效性等进行实质审查;无论平台是否已进行审查,用户均应自行承担因其传输的作品、信息、内容而可能或已经产生的侵权或权属纠纷等法律责任。本站所有资源不代表本站的观点或立场,基于网友分享,根据中国法律《信息网络传播权保护条例》第二十二与二十三条之规定,若资源存在侵权或相关问题请联系本站客服人员,点此联系我们。关于更多版权及免责申明参见 版权及免责申明
网友评论
我要评论