如何确保Abaqus中随机生成的纤维不与修剪圆相交？

我使用Python编写代码在Abaqus中生成随机纤维，希望确保生成的纤维不会与左侧的修剪圆相交，但当前的相交检测逻辑并未生效。

我已尝试修改代码，添加了纤维与修剪圆的相交检测逻辑：生成新纤维后，同时检查其与现有纤维及修剪圆是否相交，若检测到相交则标记intersection_exists为True，本应排除这类纤维，但问题仍存在。

代码如下 #

from abaqus import *
from abaqusConstants import *
from caeModules import *
from driverUtils import executeOnCaeStartup
from random import *
from math import *

executeOnCaeStartup()

# Function to check the intersection of a circle with a list of circles
def intersection(R1, Cen1, R2, Cen2, Tol):
    Distan = ((Cen1[0] - Cen2[0]) ** 2.0 + (Cen1[1] - Cen2[1]) ** 2.0) ** 0.5
    Min_Dis = R1 + R2 + Tol

    if Distan > Min_Dis:
        Intersection = 'No'
    else:
        Intersection = 'Yes'

    return Intersection

def ran(X1, X2):
    value = (X2 - X1) * random() + X1
    return value

Wcomp = 2.0  # Adjust the matrix width
Hcomp = 6.0  # Adjust the matrix height
Rmin = 0.1
Rmax = 0.1
Tol = 0.001

Vf_target = 0.3  # Target fiber volume fraction (30%)
Vff = 0.0

CentList = []

while Vff < Vf_target:
    CentList = []

    Vff = 0.0

    while Vff < Vf_target:
        Newcircle = [ran(Rmax, Wcomp - Rmax), ran(Rmax, Hcomp - Rmax), ran(Rmin, Rmax)]

        # Check for intersection with existing fibers
        intersection_exists = False
        for Cent in CentList:
            R = Cent[2]
            Check = intersection(R, Cent, Newcircle[2], Newcircle, Tol)
            if Check == 'Yes':
                intersection_exists = True
                break

        # Check for intersection with matrix boundaries
        if (Newcircle[0] - Newcircle[2] < 0 + Tol or
                Newcircle[0] + Newcircle[2] > Wcomp - Tol or
                Newcircle[1] - Newcircle[2] < 0 + Tol or
                Newcircle[1] + Newcircle[2] > Hcomp - Tol):
            intersection_exists = True

        if not intersection_exists:
            CentList.append(Newcircle)
            Vff += pi * Newcircle[2] ** 2.0 / (Wcomp * Hcomp)

        # Check for intersection with trim circle
        trim_circle_center = (0.0, Hcomp / 2)
        trim_circle_radius = 3.375
        trim_check = intersection(Newcircle[2], Newcircle, trim_circle_radius, trim_circle_center, Tol)
        if trim_check == 'Yes':
            intersection_exists = True

        if not intersection_exists:
            CentList.append(Newcircle)
            Vff += pi * Newcircle[2] ** 2.0 / (Wcomp * Hcomp)


    print(Vff)

print(CentList)

# Create Part (Matrix)
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=20.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=(Wcomp, Hcomp))
s.CircleByCenterPerimeter(center=(0.0, Hcomp/2), point1=(0.0, 3.375))
s.CoincidentConstraint(entity1=v[4], entity2=g[2], addUndoState=False)
s.EqualDistanceConstraint(entity1=v[0], entity2=v[1], midpoint=v[4], 
    addUndoState=False)
s.autoTrimCurve(curve1=g[6], point1=(-0.232322692871094, 2.57432436943054))
s.autoTrimCurve(curve1=g[2], point1=(-0.0296030044555664, 3.15962839126587))
p = mdb.models['Model-1'].Part(name='Plate', dimensionality=TWO_D_PLANAR, 
    type=DEFORMABLE_BODY)
p = mdb.models['Model-1'].parts['Plate']
p.BaseShell(sketch=s)
s.unsetPrimaryObject()
session.viewports['Viewport: 1'].setValues(displayedObject=p)
del mdb.models['Model-1'].sketches['__profile__']

# Partition
f, e, d1 = p.faces, p.edges, p.datums
t = p.MakeSketchTransform(sketchPlane=f[0], sketchPlaneSide=SIDE1, origin=(0.0, 0.0, 0.0))
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=1.0, gridSpacing=0.01, transform=t)
g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
s.setPrimaryObject(option=SUPERIMPOSE)
p = mdb.models['Model-1'].parts['Plate']

# Generate random fibers
for Cent in CentList:
    X = Cent[0]
    Y = Cent[1]
    R = Cent[2]
    p.projectReferencesOntoSketch(sketch=s, filter=COPLANAR_EDGES)
    s.CircleByCenterPerimeter(center=(X, Y), point1=(X + R, Y))

p = mdb.models['Model-1'].parts['Plate']
f = p.faces
pickedFaces = f.getSequenceFromMask(mask=('[#1 ]',), )
pickedFaces = f[:]
e1, d2 = p.edges, p.datums
p.PartitionFaceBySketch(faces=pickedFaces, sketch=s)
s.unsetPrimaryObject()
del mdb.models['Model-1'].sketches['__profile__']



#Property

#Matrix
mdb.models['Model-1'].Material(name='Matrix')
mdb.models['Model-1'].materials['Matrix'].Elastic(table=((3000.0, 0.3), ))
mdb.models['Model-1'].HomogeneousSolidSection(name='Section-2', 
    material='Matrix', thickness=None)
p = mdb.models['Model-1'].parts['Plate']
f = p.faces
faces = f.findAt((0,0,0))
q=faces.index
Fac = f[q:q+1]
region = p.Set(faces=Fac, name='Set-Matrix')
p = mdb.models['Model-1'].parts['Plate']
p.SectionAssignment(region=region, sectionName='Section-2', offset=0.0, 
    offsetType=MIDDLE_SURFACE, offsetField='', 
    thicknessAssignment=FROM_SECTION)

#fibers
mdb.models['Model-1'].Material(name='Fiber')
mdb.models['Model-1'].materials['Fiber'].Elastic(table=((70000.0, 0.3), ))
mdb.models['Model-1'].HomogeneousSolidSection(name='Section-1', 
    material='Fiber', thickness=None)
f = p.faces

i = 1
for Cent in CentList:
    X = Cent[0]
    Y = Cent[1]
    R = Cent[2]
    
# Find faces and assign section to each fiber
for i, Cent in enumerate(CentList):
    X = Cent[0]
    Y = Cent[1]
    R = Cent[2]
    
    faces = f.findAt((X, Y, 0))
    if faces:
        q = faces.index
        Fac = f[q:q+1]
        region = p.Set(faces=Fac, name='Set-1'+str(i+1))
        p.SectionAssignment(region=region, sectionName='Section-1', offset=0.0,
                            offsetType=MIDDLE_SURFACE, offsetField='',
                            thicknessAssignment=FROM_SECTION)
    else:
        print("No faces found at coordinates ({X}, {Y})".format(X=X, Y=Y))

内容的提问来源于stack exchange，提问作者awjaf