许可优化
软件
产品
本文所阐述的过程为在给定区域内生成随机纤维的方法，随机变量包括纤维的长度，随机区域的大小等，算法不拘泥于某一特定的编程语言或者计算软件，只是在文末给出了便于在Abaqus内应用的Python源码，采用该Python源码可以实现在Abaqus内生成任意需求的纤维分布，源码涉及到本人的知识产权，仅供科研交流和学习使用（包括论文写作），请勿对其进行二次贩售，本文转载时请注明出处：【知乎-蓝莓黑加仑】。
纤维混凝土建立过程中包含的两部分为：
1.外层混凝土的尺寸
2.内部纤维的参数
外层混凝土一般为规则的矩形，其建立中涉及到长（L）、宽（B）、高（H）三个尺寸参数，这里在代码中定义好相应的参数变量即可。
比较棘手的是内部纤维的参数，纤维建立中涉及到纤维长度、纤维朝向、以及多个纤维之间的空间关系三个参数。
纤维长度较为简单，也是在代码中定义好相应的参数变量即可，纤维长度本文代码中采用参数fiberlength表示。
纤维朝向一般认为具有随机性，即每个纤维默认其在空间360°内自由旋转，因此朝向这里需要引入均匀分布的随机变量，对此本文的算法是先确定纤维的一个端点，记为ipoint，接着以ipoint为起点，给定绝对长度在空间中进行360°的随机确定终点epoint。为了便于整理和重复调用，本文代码对此过程直接编了一个函数来实现，函数名为：fiberlocation。
在随机完每根纤维后，需要判定该纤维是否超出混凝土边界，对于超出边界的还需要删除重新生成。这个判定条件很简单，因为随机起点的时候肯定是在给定区域内随机获得的，所以起点一定没有超出混凝土边界，故仅需要判定纤维终点是否在边界之外即可，对比终点坐标与混凝土边界坐标的关系来得出。
多个纤维之间的空间关系主要是为了保证多个纤维之间不发生相交，这也符合实际中纤维的摆放要求，在这里需要引入的是空间直线不相交的判断方法，即如何判断任意空间直线是否相交。我们都知道判断两个空间直线是否相交最常规的方法是判断两个直线方程所构成的方程组是否有解，但这个方式未免过于低效，这意味着建立数百根纤维时候将会有成倍的方程组需要求解。
现在发挥你的想象力，究竟有什么快捷的方法可以直接判断空间直线是否相交，这种方法最好不要牵扯到多大的数值或解方程问题，因为这会使得求解过程过于漫长从而无法应用于实战。
最好可以将直线是否相交转化为正负号问题，这样仅需判断数值的正负就可以来判断相交问题了，能有这个想法真的太赞了！如果你也想到这个点子，不妨奖励自己一款3A大作犒劳一下大脑，写这篇文章时候马上就是Steam夏促了，E胖让我剁手的时间又到了，对尼尔：机械纪元的续作尼尔：伪装者，我本人很感兴趣，在机械纪元里我花费了超过400小时来达成全成就，这是一款很棒的游戏，我也把它推荐给你。说的有点跑题了，我们回到相交的判别上来，本文对此采用了矢量叉乘与点乘相结合的方法来进行直线是否相交的快速判别，我们来具体看看怎么进行数学建模吧。
如果两个直线相交时会有什么情况发生？我们做一个简图来进行说明。
相交时，对于其中一条直线而言，另一直线的两端点必在其两侧，两侧说明了什么呢？
让我们把图画得更直观一点，看看引入向量以后会产生什么不同。
天哪！结果显而易见，针对不同端点进行叉乘所得出的向量方向相反，对于上部小图而言，叉乘后向量指向下，而对于下部则是叉乘后指向上，叉乘所得结果表明所得向量的方向恰好相反！
所以我们理所当然可以对叉乘所得的向量进行点乘，因为如果他们方向相反，那么点乘中的cos值必为-1，嗯，就是这样，依据cos值是否为-1就可以判断两直线是否相交。
那么，有没有特殊情况呢？会使得这种判别方式失效呢？答案是肯定的，请看下图两纤维的摆放情况，它的计算结果是不是和前述的一致呢？但是很明显这两个纤维并没有产生相交。
如何将这种特殊情况排除呢？要知道如果把本不相交的纤维判定为相交，我们就需要重新再计算一个纤维，直到他们不相交为止，对此我们应尽可能地避免将不相交归到相交这一类中，这对大量纤维建模时会很有帮助。
解决方法很简单，只需以另一条直线为基准重复一次上述操作即可，这样如果两次计算的cos值都为-1时，那理所当然是相交关系，如果仅一次计算的cos值为-1，即图中那个特殊情况，那应判断为不相交，这个纤维可以使用！
以下是上述过程应用到实战的Python源码，将下述代码复制到记事本，保存为.py结尾的文件，在Abaqus中点击左上角的File-Run script-选择保存好后的.py文件运行，模型即可自动生成。
# coding=UTF-8
from abaqusConstants import *
from caeModules import *
import numpy as np
import math
def intersection (a,b):
a1 = list(a[0])
a2 = list(a[1])
b1 = list(b[0])
b2 = list(b[1])
v1 = np.cross(np.array(a1) - np.array(b1), np.array(b2) - np.array(b1))
v2 = np.cross(np.array(a2) - np.array(b1), np.array(b2) - np.array(b1))
dot_value = np.dot(v1, v2)/(np.linalg.norm(v1) * np.linalg.norm(v2))
v11 = np.cross(np.array(b1) - np.array(a1), np.array(a2) - np.array(a1))
v21 = np.cross(np.array(b2) - np.array(a1), np.array(a2) - np.array(a1))
dot_value1 = np.dot(v11, v21)/(np.linalg.norm(v11) * np.linalg.norm(v21))
if dot_value == -1 and dot_value1 == -1:
return 1
else:
return 0
def fiberlocation (b):
Psi = np.random.uniform(0,2 * math.pi)
Theta = np.random.uniform(0,2 * math.pi)
x = b[0] + fiberlength * math.sin(Psi) * math.sin(Theta)
y = b[1] + fiberlength * math.sin(Psi) * math.cos(Theta)
z = b[2] + fiberlength * math.cos(Psi)
endpoint = (x,y,z)
return endpoint
# 输入投放区域的长、高、宽，投放区域纤维个数与纤维长度值(这5个参数根据自身需求选取)
L = 150.0
H = 150.0
B = 550.0
fiberlength = 12
fibernumber = 800
fiber = [] * fibernumber
initialfiber = (np.random.uniform(-(H/2-fiberlength),(H/2-fiberlength)),np.random.uniform(-(B/2-fiberlength),(B/2-fiberlength)),np.random.uniform(-(L/2-fiberlength),(L/2-fiberlength)))
fiber.append((initialfiber,fiberlocation(initialfiber)))
i = 1
while i < fibernumber:
ipoint = (np.random.uniform(-H/2,H/2),np.random.uniform(-B/2,B/2),np.random.uniform(-L/2,L/2))
epoint = fiberlocation(ipoint)
if abs(epoint[0])<H/2 and abs(epoint[1])<B/2 and abs(epoint[2])<L/2:
for j in range(0, i):
signal = intersection(fiber[j], (ipoint, epoint))
if signal is True:
break
fiber.append((ipoint, epoint))
i = i + 1
else:
continue
myModel = mdb.Model(name='Model-1')
p = myModel.Part(name='fiber', dimensionality=THREE_D, type=DEFORMABLE_BODY)
p.ReferencePoint(point=(0.0, 0.0, 0.0))
for j in range (0,fibernumber):
p.WirePolyLine(points=(fiber[j][0], fiber[j][1]), mergeType=IMPRINT, meshable=ON)
del p.features['RP']
#构建混凝土部分
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
s.setPrimaryObject(option=STANDALONE)
s.Line(point1=(-H/2, -B/2), point2=(-H/2, B/2))
s.VerticalConstraint(entity=g[2], addUndoState=False)
s.Line(point1=(-H/2, B/2), point2=(H/2, B/2))
s.HorizontalConstraint(entity=g[3], addUndoState=False)
s.PerpendicularConstraint(entity1=g[2], entity2=g[3], addUndoState=False)
s.Line(point1=(H/2, B/2), point2=(H/2, -B/2))
s.VerticalConstraint(entity=g[4], addUndoState=False)
s.PerpendicularConstraint(entity1=g[3], entity2=g[4], addUndoState=False)
s.Line(point1=(H/2, -B/2), point2=(-H/2, -B/2))
s.HorizontalConstraint(entity=g[5], addUndoState=False)
s.PerpendicularConstraint(entity1=g[4], entity2=g[5], addUndoState=False)
p = mdb.models['Model-1'].Part(name='Concrete', dimensionality=THREE_D, type=DEFORMABLE_BODY)
p = mdb.models['Model-1'].parts['Concrete']
p.BaseSolidExtrude(sketch=s, depth=L)
s.unsetPrimaryObject()
session.viewports['Viewport: 1'].setValues(displayedObject=p)
del mdb.models['Model-1'].sketches['__profile__']
#组合拼装
a = mdb.models['Model-1'].rootAssembly
a.DatumCsysByDefault(CARTESIAN)
p = mdb.models['Model-1'].parts['Concrete']
a.Instance(name='Concrete-1', part=p, dependent=ON)
a.translate(instanceList=('Concrete-1', ), vector=(0.0, 0.0, -L/2))
p = mdb.models['Model-1'].parts['fiber']
a.Instance(name='fiber-1', part=p, dependent=ON)
此脚本仅仅是完成几何建模，具体计算的参数需要单独输入，另外纤维与混凝土之间的粘结关系其实并不好确定，按照以往的文献来看可以定义为embed约束 。这里算例就取的比较随意了，取比较常规的立方体，各个参数不去单独说明，自己做材性试验或者部分文献里都有，仅展示一下最终计算结果。
关注我，学习更多Abaqus的知识
技术文档
热门文章
155-2731-8020
