# -*- coding: utf-8 -*-
from abaqus import *
from abaqusConstants import *
from caeModules import *
from odbAccess import *
import math
import csv
import os
import numpy as np
from scipy.stats import qmc
def run_fem_simulation(params):
try:
tool_position = params['tool_position']
tool_step = params['tool_step']
ae = params['ae']
R = params['R']
# ---------------几何参数定义----------------
ap = 10.0 # 轴向切深
ae = 3.0 # 径向切深
ft = 0.5 # 每齿进给量
min_mesh = 0.5 # 切削区网格尺寸
max_mesh = 3.0 # 非切削区网格尺寸
# ---------------材料参数定义----------------
length = 120.0 # 工件长度(X方向)
width = 60.0 # 工件宽度(Z方向)
depth = 5.0 # 工件厚度(Y方向)
material_name = 'TC4'
E = 113800.0 # 弹性模量 (MPa)
nu = 0.33 # 泊松比
density = 4.43e-6 # 密度 (kg/mm³)
# ---------------刀具参数定义----------------
R = 6.0 # 刀具半径
helix_angle1 = 38 # 螺旋角1,单位:度
helix_angle2 = 41
E_tool = 2.1e5 # MPa
I = math.pi * R ** 4 / 4 # 圆截面惯性矩
# ---------------切削力系数------------------
kt1, kte1 = 2026.48, 7.16 # 切向力系数
kn1, kne1 = 972.23, 16.99 # 法向力系数
ka1, kae1 = 764.08, 0.51 # 轴向力系数
kt2, kte2 = 2294.84, 10.69 # 切向力系数
kn2, kne2 = 988.74, 18.20 # 法向力系数
ka2, kae2 = 1015.45, 1.39 # 轴向力系数
N = 4 # 刀具齿数
omega = 5000 # 主轴转速 (rpm)
steps = 64 # 旋转步数
phis = math.pi / 2 - math.acos(1 - ae / R)
phie = math.pi / 2
convergence_threshold = 0.1 # 收敛阈值
model_name = 'Model-%d-%d' % (tool_position, tool_step)
if model_name in mdb.models.keys():
del mdb.models[model_name]
myModel = mdb.Model(name=model_name)
# ---------------几何建模---------------------
# Part创建
mySketch = myModel.ConstrainedSketch(name='sketch', sheetSize=200)
mySketch.rectangle((0, 0), (length, depth))
myPart = myModel.Part(name='WorkPiece', dimensionality=THREE_D,
type=DEFORMABLE_BODY)
myPart.BaseSolidExtrude(sketch=mySketch, depth=width)
# 分割工件(用于不同区域网格划分)
datum = myPart.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=width / 2)
myPart.PartitionCellByDatumPlane(
datumPlane=myPart.datums[datum.id],
cells=myPart.cells.getSequenceFromMask(mask=('[#1 ]',), )
)
# ---------------网格划分---------------------
# 切削区域(上半部分)精细网格
region_edges = myPart.edges.getByBoundingBox(
xMin=-1e-5, xMax=length + 1e-5,
yMin=-1e-5, yMax=depth + 1e-5,
zMin=width / 2 - 1e-5, zMax=width + 1e-5
)
myPart.seedEdgeBySize(edges=region_edges, size=min_mesh, constraint=FINER)
# 非切削区域(下半部分)粗网格
other_edges = myPart.edges.getByBoundingBox(
xMin=-1e-5, xMax=length + 1e-5,
yMin=-1e-5, yMax=depth + 1e-5,
zMin=-1e-5, zMax=width / 2 + 1e-5
)
myPart.seedEdgeBySize(edges=other_edges, size=max_mesh, constraint=FINER)
myPart.generateMesh()
# ---------------材料定义---------------------
myMaterial = myModel.Material(name=material_name)
myMaterial.Density(table=((density,),))
myMaterial.Elastic(table=((E, nu),))
# 截面属性
mySection = myModel.HomogeneousSolidSection(name='Section',
material=material_name)
region = myPart.Set(cells=myPart.cells[:], name='AllCells')
myPart.SectionAssignment(region=region, sectionName='Section')
# ---------------装配体操作-------------------
myAssembly = myModel.rootAssembly
instance_name = 'WorkPiece-Instance'
myInstance = myAssembly.Instance(name=instance_name, part=myPart,
dependent=ON)
# ---------------边界条件设置-----------------
# 固定底面(Z=0)
# 选择在 xz 平面上的所有节点
nodes = myInstance.nodes
fixed_nodes = nodes.getByBoundingBox(xMin=-1e-05, xMax=length,
yMin=-1e-05, yMax=depth,
zMin=-1e-05, zMax=1e-05)
fixed_node_set = myAssembly.Set(nodes=fixed_nodes,
name='FixedNodeSet'
)
myModel.DisplacementBC(name='BC-FIX',
createStepName='Initial',
region=fixed_node_set,
u1=0.0, u2=0.0, u3=0.0
)
# ---------------分析步设置-------------------
if tool_position > 1:
myModel.StaticStep(name='Step-del',
previous='Initial',
maxNumInc=200,
timePeriod=1e-05,
initialInc=1e-10,
minInc=1e-10,
maxInc=1e-05
)
myModel.StaticStep(name='Step-1',
previous='Step-del',
nlgeom=ON,
maxNumInc=200,
timePeriod=0.1,
initialInc=1e-05,
minInc=1e-05,
maxInc=1e-02
)
else:
myModel.StaticStep(name='Step-1',
previous='Initial',
nlgeom=ON,
maxNumInc=200,
timePeriod=0.1,
initialInc=1e-05,
minInc=1e-05,
maxInc=1e-02
)
###---------------------------------------------创建被切削区域集合-----------------------------------------------------
# 切削区节点坐标
xmin = -1e-05
xmax = length
ymin = -1e-05
ymax = ae + 1e-05
zmin = width - ap + 1e-05
zmax = width + 1e-05
# 创建被切削节点集合
nodes = myInstance.nodes
milling_nodes = nodes.getByBoundingBox(xMin=xmin, xMax=xmax,
yMin=ymin, yMax=ymax,
zMin=zmin, zMax=zmax)
node_set_name = 'MillingNodes'
myAssembly.Set(nodes=milling_nodes, name=node_set_name)
print(f'Created node set [{node_set_name}] containing {len(milling_nodes)} nodes')
# 创建被切削单元集合
elements = myInstance.elements
milling_elements = elements.getByBoundingBox(xMin=xmin, xMax=xmax + min_mesh,
yMin=ymin, yMax=ymax,
zMin=zmin - min_mesh, zMax=zmax)
element_set_name = 'MillingElements'
myAssembly.Set(elements=milling_elements, name=element_set_name)
print(f'Created element set [{element_set_name}] containing {len(milling_elements)} elements')
#生死单元单元格集合
milling_cells = milling_elements.getByBoundingBox(xMin=- 1e-5, xMax=(tool_position - 1) * R + 1e-5,
yMin=ymin, yMax=ymax,
zMin=zmin - min_mesh, zMax=zmax)
removed_name = f'remove-P{tool_position}'
myAssembly.Set(elements=milling_cells, name=removed_name)
print(f'Created removed feature set [{removed_name}] containing {len(milling_cells)} elements')
#变形节点集合
deformation_nodes = nodes.getByBoundingBox(xMin=(tool_position - 1) * R - min_mesh/2,
xMax=(tool_position - 1) * R + min_mesh/2,
yMin=ae - min_mesh + min_mesh/2,
yMax=ae + min_mesh - min_mesh/2,
zMin=width - ap + tool_step * min_mesh - min_mesh/2,
zMax=width - ap + tool_step * min_mesh + min_mesh/2)
if deformation_nodes:
deformation_set_name = f'DEFORMATION-P{tool_position}-M{tool_step}'
myAssembly.Set(nodes=deformation_nodes, name=deformation_set_name)
print(f'Created deformation set [{deformation_set_name}] with {len(deformation_nodes)} nodes')
else:
print(f'[Warning] No deformation nodes found for tool position {tool_position}')
###----------------------------------------------计算并施加载荷---------------------------------------------------------
convergence_threshold = 0.2
displacement = 0
previous_ae = 0
iteration_number = 1
if tool_position % 2 == 1:
helix_angle = helix_angle1
else:
helix_angle = helix_angle2
dphi = math.pi * 2 / steps
if helix_angle != 0:
db = R * dphi / math.tan(math.radians(helix_angle))
else:
db = ap / steps
steps_axial = int(ap / db + 1)
milling_set = {}
for node in milling_nodes:
node_id = node.label
node_x, node_y, node_z = node.coordinates
milling_set[node_id] = (node_x, node_y, node_z)
loaded_node_coords = []
# === 计算所有齿中底部微元(cnt4=0)的最小切出步 ===
dphi = math.pi * 2 / steps # 旋转步长
min_exit_step = steps # 初始化为最大值
for cnt3 in range(N):
# 计算当前齿底部微元切出所需步数
phia_needed = math.pi / 2 - phis - cnt3 * 2 * math.pi / N
# 调整phia_needed到0-2π范围内
while phia_needed < 0:
phia_needed += 2 * math.pi
while phia_needed >= 2 * math.pi:
phia_needed -= 2 * math.pi
# 计算切出步数
cnt1_exit = math.ceil(phia_needed / dphi)
# 更新最小切出步
if cnt1_exit < min_exit_step:
min_exit_step = cnt1_exit
while abs(ae - previous_ae) > convergence_threshold:
Fy_total = [0.0] * steps_axial
zj_list = [width - ap + j * db for j in range(steps_axial)] # 每个cnt4对应的Z轴坐标
delta_y_current = [0.0] * steps_axial # 当前轮的刀具变形,初始化为0
loaded_node_coords = [] # 清空加载节点坐标
for cnt1 in range(steps):
if cnt1 == (min_exit_step - 1) + tool_step:
for cnt3 in range(N):
if cnt3 in [0, 2]: # 齿编号0和2使用第一组系数
kt, kte = kt1, kte1
kn, kne = kn1, kne1
ka, kae = ka1, kae1
else: # 齿编号1和3使用第二组系数
kt, kte = kt2, kte2
kn, kne = kn2, kne2
ka, kae = ka2, kae2
for cnt4 in range(steps_axial):
phia = phis + cnt1 * dphi + cnt3 * 2 * math.pi / N - cnt4 * db * math.tan(
math.radians(helix_angle)) / R
phia = phia % (2 * math.pi)
if phis <= phia <= phie:
h = ft * math.sin(phia)
Ft = kt * h * db + kte * db
Fn = kn * h * db + kne * db
Fa = ka * h * db + kae * db
else:
Ft, Fn, Fa = 0, 0, 0
# 将切削力分解到直角坐标系
Fx = -Ft * math.cos(phia) - Fn * math.sin(phia)
Fy = Ft * math.sin(phia) - Fn * math.cos(phia)
Fz = Fa
# 累加Fy
Fy_total[cnt4] += Fy
for i in range(int(length / R) - 1):
# 计算刀具坐标系切削微元位置
micro_x = R * math.cos(phia) + (tool_position - 1) * R
micro_y = R * math.sin(phia) - (depth - ae) - delta_y_current[cnt4]
micro_z = width - ap + cnt4 * db
# 找到最近的节点
min_distance = float('inf')
nearest_node = None
for node_id, (node_x, node_y, node_z) in milling_set.items():
distance = math.sqrt(
(micro_x - node_x) ** 2 + (micro_y - node_y) ** 2 + (micro_z - node_z) ** 2)
if distance < min_distance:
min_distance = distance
nearest_node = node_id
# 创建载荷
if nearest_node is not None and (Fx != 0 or Fy != 0 or Fz != 0):
node = milling_nodes.sequenceFromLabels((nearest_node,))
node_name = 'Load-' + str(cnt1) + '-' + str(cnt3) + '-' + str(cnt4)
load_name = 'Force-' + str(cnt1) + '-' + str(cnt3) + '-' + str(cnt4)
node_set = myAssembly.Set(nodes=node, name=node_name)
myModel.ConcentratedForce(name=load_name, createStepName='Step-1',
region=node_set, cf1=Fx, cf2=Fy, cf3=Fz)
node_x, node_y, node_z = milling_set[nearest_node]
loaded_node_coords.append((node_x, node_y, node_z))
# 创建参考点
adjusted_point = (micro_x, micro_y, micro_z)
reference_point_name = 'RP-' + str(cnt1) + '-' + str(cnt3) + '-' + str(cnt4)
reference_point = myAssembly.ReferencePoint(point=adjusted_point)
# 刀具变形计算
delta_y_new = [] # 创建新列表存储本轮的变形计算结果
for j in range(steps_axial):
zj = zj_list[j]
delta = 0.0
for jp in range(steps_axial):
zjp = zj_list[jp]
if zj <= zjp:
K = zjp ** 2 * (3 * zj - zjp) / (6 * E_tool * I)
else:
K = zj ** 2 * (3 * zjp - zj) / (6 * E_tool * I)
delta += K * Fy_total[jp]
delta_y_new.append(delta)
# 更新delta_y_current用于下一次迭代(只在有计算值时更新)
if delta_y_new: # 确保列表不为空
delta_y_current = delta_y_new
###------------------------------------------------------生死单元------------------------------------------------
if tool_position > 1:
intname = 'Int-' + str(tool_position)
stepname = 'Step-del'
element_name = removed_name
int_region = myAssembly.sets[element_name]
myModel.ModelChange(name=intname,
createStepName=stepname,
region=int_region,
regionType=ELEMENTS,
activeInStep=False,
includeStrain=False
)
###-------------------------------------------------------创建作业---------------------------------------------------
job_name = f'Model-{tool_position}-{tool_step}-{iteration_number}'
myJob = mdb.Job(name=job_name, model=model_name, type=ANALYSIS, numCpus=24, numDomains=24)
myJob.submit()
mdb.jobs[job_name].waitForCompletion()
###-----------------------------------------------------打开ODB文件--------------------------------------------------
try:
odb = openOdb(path=job_name + '.odb')
step_name = 'Step-1'
displacement_field = odb.steps[step_name].frames[-1].fieldOutputs['U']
instance_name = 'WORKPIECE-INSTANCE'
instance = odb.rootAssembly.instances[instance_name]
deformation_name = f'DEFORMATION-P{tool_position}-M{tool_step}'
if deformation_name in odb.rootAssembly.nodeSets:
region = odb.rootAssembly.nodeSets[deformation_name]
displacement_subset = displacement_field.getSubset(region=region)
displacement = displacement_subset.values[0].data[2]
else:
print(f"[Critical] Node set {deformation_name} missing in ODB!")
displacement = 0 # 默认值或终止
print(f'ae={ae},displacement={displacement}')
except Exception as e:
print(f"Error reading ODB: {str(e)}")
displacement = 0
finally:
odb.close()
pre_ae = ae
ae -= abs(displacement)
ae = (ae + previous_ae)/2
previous_ae = pre_ae
if iteration_number > 10:
break
else:
iteration_number += 1
final_displacement = displacement # 从ODB获取最终位移
final_ae = ae # 最终径向切深
print('final ae is:', final_ae)
print('final displacement is:', final_displacement)
return {
'displacement': final_displacement,
'final_ae': final_ae,
'status': 'success',
'message': 'Simulation completed successfully'
}
except Exception as e:
return {
'status': 'error',
'message': str(e),
'displacement': None,
'final_ae': None
}
def generate_parameters(samples):
"""生成拉丁超立方采样参数组合"""
param_ranges = {
'tool_position': (1, 19),
'tool_step': (1, 11)
}
# 拉丁超立方采样
sampler = qmc.LatinHypercube(d=len(param_ranges))
sample_data = sampler.random(n=samples)
# 构造参数列表
params_list = []
for sample in qmc.scale(sample_data,
[v[0] for v in param_ranges.values()],
[v[1] for v in param_ranges.values()]):
params = {
'tool_position': int(round(sample[0])),
'tool_step': int(round(sample[1])),
'ae': 3.0, # 固定值
'R': 5.0 # 固定值
}
params_list.append(params)
return params_list
# 数据保存模块
class DataRecorder:
def __init__(self, filename="simulation_data.csv"):
self.filename = filename
self._init_file()
def _init_file(self):
"""初始化数据文件"""
if not os.path.exists(self.filename):
with open(self.filename, 'w', newline='') as f:
writer = csv.writer(f)
writer.writerow([
'tool_position', 'tool_step',
'ae_input', 'R_input',
'displacement', 'final_ae',
'status'
])
def save_record(self, params, result):
"""保存单条记录"""
record = {
'tool_position': params['tool_position'],
'tool_step': params['tool_step'],
'ae_input': params['ae'],
'R_input': params['R'],
'displacement': result.get('displacement', None),
'final_ae': result.get('final_ae', None),
'status': result['status']
}
with open(self.filename, 'a', newline='') as f:
writer = csv.DictWriter(f, fieldnames=record.keys())
writer.writerow(record)
if __name__ == '__main__':
# 初始化数据记录器
recorder = DataRecorder()
# 生成参数组合
param_combinations = generate_parameters(samples=200)
# 执行参数扫描
for idx, params in enumerate(param_combinations,1):
print(f"\n 正在执行参数组合 {idx}/{len(param_combinations)}")
print(f" 刀具位置: {params['tool_position']}")
print(f" 刀具旋转步: {params['tool_step']}")
print(f" 径向切深: {params['ae']}mm")
print(f" 刀具半径: {params['R']}mm")
# 执行仿真
result = run_fem_simulation(params)
# 处理结果
if result['status'] == 'success':
print(f"成功完成: 位移={result['displacement']}mm, 最终ae={result['final_ae']}mm")
else:
print(f"执行失败: {result['message']}")
# 保存数据
recorder.save_record(params, result)
print("\n所有参数组合执行完成!")
print(f"数据文件已保存至: {os.path.abspath(recorder.filename)}")
以上是我的py代码,在abaqus中运行速度特别缓慢
最新发布