Waitwind'blog

随机纤维投放-abaqus二次开发

随机纤维投放-abaqus二次开发

整个投放过程就是不断的生成随机纤维，同时保证生成的纤维与之前生成的纤维不相交。相交的判断逻辑我用了两种方式实现，对应下面两段代码。第一种是从几何上判断两条线段是否相交，对已生成的所有线段判断。第二种是利用了abaqus自带的布尔操作，在进行cut布尔操作时，不相交会报错，使用try-except将错误捕捉后，就能对相交进行判断，理论上第二种方式对于所有的相交问题都能处理。

下面是两张生成的模型图

  • image-20240430153721091
  • image-20240430153830930

第一种

#-*- coding:utf-8 –*-
'''
随机投放纤维
copyright waitwind 2024年4月29日
'''
import numpy as np
from abaqus import *
from abaqusConstants import *
import random as rd
import datetime

# 记录开始时间
start_time = datetime.datetime.now()

#mdb=Mdb()
model=mdb.Model('mymodel')
part=model.Part('steel fibre')

a,b,h,length,area,pho_s,pho_c=100,100,800,30,0.5**2*np.pi,7800,2400
v_wire=area*length
m_wire=v_wire*pho_s
V=a*b*h
M=V*pho_c
rat=1.5e-2
n=int(V*rat/v_wire)
wires=[]
'''
参数：
a,b,c：长，宽，高
length：纤维长度 10-60mm
area：纤维面积
v_wire,m_ire:钢纤维体积，钢纤维质量
V,M:混凝土体积，混凝土质量
pho_s,pho_c:钢密度，混凝土密度
rat：体积分数，或者质量分数
n:钢纤维数量
'''
def isxiangjiao(A, B, C, D):
    """
    判断三维空间中的两条线段是否相交。

    参数:
    A, B: 线段1的起点和终点 (列表形式)
    C, D: 线段2的起点和终点 (列表形式)

    返回:
    True 线段不相交
    False 线段相交
    """
    # 将输入的点转换为numpy数组
    A, B, C, D = np.array(A), np.array(B), np.array(C), np.array(D)

    # 计算方向向量
    u = B - A
    v = D - C
    w = C - A
    x = D - A

    # 检查四个点是否共面
    if np.isclose(np.dot(np.cross(u, v), w), 0):

        # 平行且共线
        if np.allclose(np.cross(u,v), 0) and np.allclose(np.cross(u,w),0):
            #计算线段两点到线段另一点的距离，线段长度均小于距离，不相交，返回True
            if np.dot(u,u)<np.dot(w,w) and np.dot(u,u)<np.dot(x,x):
                return True
            else:
            #相交返回False
                return False
          
        # 仅平行不共线，必不相交，返回True。
        elif np.allclose(np.cross(u,v),0):
            return True
      
        #不平行，必有交点，交点不一定在两条线段内。
        else:
            #计算交点位置
            if np.allclose(np.cross([u[0],-v[0]],[u[1],-v[1]]),0):#是否线性相关
                p=np.array([[u[2],-v[2]],[u[1],-v[1]]])
                q=np.array([C[2]-A[2],C[1]-A[1]])
            else:#线性无关
                p=np.array([[u[0],-v[0]],[u[1],-v[1]]])
                q=np.array([C[0]-A[0],C[1]-A[1]])
            #解直线方程
            solution=np.linalg.solve(p,q)
          
            if 0<=solution[0]<=1 and 0<=solution[1]<=1:
                #交点位置在线段内，相交
                return False
            else:
                #交点位置在线段外，不相交
                return True

    # 不共面必不相交，返回True
    else:
        return True


def createWire():
    x_start,y_start,z_start=rd.uniform(0,a),rd.uniform(0,b),rd.uniform(0,h)#起点坐标
    while True:
        theta,phi=rd.uniform(0,2*np.pi),rd.uniform(0,np.pi)
        x_end,y_end,z_end=x_start+length*np.sin(phi)*np.cos(theta),y_start+length*np.sin(phi)*np.sin(theta),z_start+length*np.cos(phi)
        wirePoint = ((x_start, y_start, z_start), (x_end, y_end, z_end))
        if all((0<=x_end<=a,0<=y_end<=b,0<=z_end<=h)):
            if not wires:
                break
            values=(isxiangjiao(wires[i][0],wires[i][1],wirePoint[0],wirePoint[1]) for i in range(len(wires)))
            if all(values):
                break
    part.WirePolyLine(points=(wirePoint, ),
        mergeType=IMPRINT, meshable=ON)
    return wirePoint
def createBeam():
    s = model.ConstrainedSketch(name='__profile__',sheetSize=200.0)
    g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
    s.setPrimaryObject(option=STANDALONE)
    s.rectangle(point1=(0.0, 0.0), point2=(a, b))
    p = mdb.models['mymodel'].Part(name='Beam', dimensionality=THREE_D,
                                   type=DEFORMABLE_BODY)
    p.BaseSolidExtrude(sketch=s, depth=h)
    s.unsetPrimaryObject()
    del mdb.models['mymodel'].sketches['__profile__']
def assembly():
    a1 = mdb.models['mymodel'].rootAssembly
    a1.DatumCsysByDefault(CARTESIAN)
    p = mdb.models['mymodel'].parts['Beam']
    a1.Instance(name='Beam-1', part=p, dependent=ON)
    p = mdb.models['mymodel'].parts['steel fibre']
    a1.Instance(name='steel fibre-1', part=p, dependent=ON)
n=1000
for i in range(n):
    wires.append(createWire())
    print('total:{0},{1}'.format(n,i+1))

createBeam()
assembly()
mdb.saveAs('Ecc.cae')

# 记录结束时间
end_time = datetime.datetime.now()

print(end_time-start_time)

第二种

#-*- coding:utf-8 –*-
'''
随机投放纤维
copyright waitwind 2024年4月29日
'''
import numpy as np
from abaqus import *
from abaqusConstants import *
import random as rd
import datetime

# 记录开始时间
start_time = datetime.datetime.now()

#mdb=Mdb()
model=mdb.Model('mymodel')
a,b,h,length,area,pho_s,pho_c=100.0,100.0,800.0,30.0,0.5**2*np.pi,7800,2400
v_wire=area*length
m_wire=v_wire*pho_s
V=a*b*h
M=V*pho_c
rat=0.5e-2
n=int(V*rat/v_wire)
a1= model.rootAssembly
'''
参数：
a,b,c：长，宽，高
length：纤维长度 10-60mm
area：纤维面积
v_wire,m_ire:钢纤维体积，钢纤维质量
V,M:混凝土体积，混凝土质量
pho_s,pho_c:钢密度，混凝土密度
rat：体积分数，或者质量分数
n:钢纤维数量
'''

def createWire():
    x_start,y_start,z_start=rd.uniform(0,a),rd.uniform(0,b),rd.uniform(0,h)#起点坐标
    while True:
        theta,phi=rd.uniform(0,2*np.pi),rd.uniform(0,np.pi)
        x_end,y_end,z_end=x_start+length*np.sin(phi)*np.cos(theta),y_start+length*np.sin(phi)*np.sin(theta),z_start+length*np.cos(phi)
        wirePoint = ((x_start, y_start, z_start), (x_end, y_end, z_end))
        if all((0<=x_end<=a,0<=y_end<=b,0<=z_end<=h)):
            break
    part=model.Part('temp')
    part.WirePolyLine(points=(wirePoint, ),mergeType=IMPRINT, meshable=ON)

def createBeam():
    s = model.ConstrainedSketch(name='__profile__',sheetSize=200.0)
    g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
    s.setPrimaryObject(option=STANDALONE)
    s.rectangle(point1=(0.0, 0.0), point2=(a, b))
    p = mdb.models['mymodel'].Part(name='Beam', dimensionality=THREE_D,
                                   type=DEFORMABLE_BODY)
    p.BaseSolidExtrude(sketch=s, depth=h)
    s.unsetPrimaryObject()
    del mdb.models['mymodel'].sketches['__profile__']
  
def steel_fibre():
    try:
        p=model.parts['temp']
        a1.Instance(name='temp-1',part=p,dependent=ON)
        p=model.parts['steel fibre']
        a1.Instance(name='steel fibre-1',part=p,dependent=ON)
        a1.InstanceFromBooleanCut(name='temp2', instanceToBeCut=a1.instances['steel fibre-1'], cuttingInstances=(a1.instances['temp-1'], ), originalInstances=DELETE)
        del a1.features['temp2-1']
        del model.parts['temp2']
        del model.parts['temp']
    except:
        p=model.parts['steel fibre']
        a1.Instance(name='steel fibre-1',part=p,dependent=ON)
        p=model.parts['temp']
        a1.Instance(name='temp-1',part=p,dependent=ON)
        a1.InstanceFromBooleanMerge(name='steel fibre2', instances=(a1.instances['steel fibre-1'], a1.instances['temp-1'], ), originalInstances=DELETE, domain=GEOMETRY)
        del a1.features['steel fibre2-1']
        del model.parts['steel fibre']
        del model.parts['temp']
        model.parts.changeKey(fromName='steel fibre2',toName='steel fibre')

def assembly():
    a1.DatumCsysByDefault(CARTESIAN)
    p = mdb.models['mymodel'].parts['Beam']
    a1.Instance(name='Beam-1', part=p, dependent=ON)
    p = mdb.models['mymodel'].parts['steel fibre']
    a1.Instance(name='steel fibre-1', part=p, dependent=ON)
n=1000
for i in range(n):
    createWire()
    if i==0:
        model.Part(name='steel fibre', objectToCopy=model.parts['temp'])
    steel_fibre()
    print('total:{0},{1}'.format(n,i+1))

createBeam()
assembly()
mdb.saveAs('Ecc.cae')

# 记录结束时间
end_time = datetime.datetime.now()

print(end_time-start_time)

最后

  1. 随机投放1000个纤维，第一种用时 0:01:24.416000 ，第二种用时 0:02:58.252000。就效率而言，似乎第一种方式更好，但第二种方式通用性更好。
  2. 对代码进行改动，可以生成多种纤维，将生成纤维改成骨料，可以实现随机骨料的投放。甚至骨料与纤维的混杂。
  3. 代码主要难点在相交的逻辑的处理，布尔操作提供了一个通用的方式来处理，想要提高代码效率需要一个合适的处理逻辑。
  4. 纤维生成区域为长方体，可以根据需要改动。
  5. 如果有更好的方式，欢迎交流。

Hi,Friend

© 2025 Waitwind'blog