DONG_all1.for(228): error #5082: Syntax error, found REAL_CONSTANT '.843750' when expecting one of: , /
*80639.843750,80601.390625,80562.492188,80523.187500,80483.507812,80443.468750,80403.109375,80362.445312,80321.500000,80280.304
-----------^
DONG_all1.for(229): error #5082: Syntax error, found REAL_CONSTANT '.875000' when expecting one of: , /
*80238.875000,80197.257812,80155.468750,80113.539062,80071.492188,80029.359375,79987.164062,79944.921875,79902.664062,79860.414
-----------^
import numpy as np
import os
def generate_final_fortran_subroutine(input_txt, output_for):
"""
生成最终优化的Fortran子程序(包含DLOAD和UTRACLOAD)
读取多圈数据(255°-285°范围内的所有数据)
"""
try:
print(f"开始处理文件: {input_txt}")
# 读取数据:第1列(时间)和第4列(法向力Fz)
print("正在读取数据文件...")
data = np.loadtxt(input_txt, usecols=(0, 3), dtype=np.float32)
print(f"读取到 {len(data)} 行数据")
# 基本参数
R = 0.457497 # 接触斑中点所在半径(m)
velocity = 400 / 3.6 # 速度(m/s)
omega = velocity / R # 角速度(rad/s)
# 计算每个数据点对应的角度(度) - 使用连续角度
print("计算每个数据点对应的角度(度)...")
angles_deg = np.degrees(data[:, 0] * omega) # 连续角度值
# 定义关键角度范围 (255°-285°)
start_angle = 255.0
end_angle = 285.0
angle_range = end_angle - start_angle
print(f"提取关键角度范围({start_angle}°-{end_angle}°)内的数据...")
# 计算每个点相对于255°的角度偏移
angles_offset = angles_deg - start_angle
# 找到所有在255°-285°范围内的点(包括多圈)
mask = (angles_offset % 360) <= angle_range
# 提取关键角度和载荷值
key_angles = angles_deg[mask].astype(np.float32)
key_times = data[mask, 0].astype(np.float32)
key_fz = data[mask, 1].astype(np.float32)
if len(key_angles) == 0:
print("错误: 未找到关键角度范围内的数据")
return False
nPoints = len(key_angles)
print(f"找到 {nPoints} 个关键角度范围内的数据点")
print("\n详细数据点对比 (时间, 角度(度), 载荷(N)):")
print("=" * 65)
print(f"{'时间(s)':<12}{'角度(度)':<12}{'载荷(N)':<12}")
print("-" * 65)
# 打印每个数据点的详细对比信息
for i in range(nPoints):
print(f"{key_times[i]:<12.6f}{key_angles[i]:<12.6f}{key_fz[i]:<12.6f}")
print("=" * 65)
# 打印关键数据预览
print("\n关键角度数据点预览:")
print(f" 起始角度: {key_angles[0]:.2f}°, Fz: {key_fz[0]:.6f}N")
print(f" 结束角度: {key_angles[-1]:.2f}°, Fz: {key_fz[-1]:.6}N")
# 生成Fortran代码(包含两个子程序)
fortran_code = f""" SUBROUTINE DLOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,
1 COORDS,JLTYP,SNAME)
C
INCLUDE 'ABA_PARAM.INC'
C
DIMENSION TIME(2),COORDS (3)
CHARACTER*80 SNAME
**************** 直接使用KINC作为索引 ****************
* 1. 直接使用增量步号KINC作为数组索引
* 2. 无需索引计算或循环
* 3. 数组从1开始与KINC匹配
***************************************************
! ==== 声明所有变量 ====
integer nPoints
parameter (nPoints = {nPoints})
real*8 angleData(nPoints) ! 存储角度(度)
real*8 FzData(nPoints) ! 法向力
real*8 Fz, Tforce, omega
real*8 x0,y0,z0,check,x,y,z,xt,yt,zt
real*8 R, b, a, velocity
real*8 pi
parameter (pi = 3.141592653589793d0)
integer idx ! 直接使用KINC作为索引
logical firstCall
save firstCall
real*8 t
real*8 current_angle_deg ! 当前数据点角度(度)
real*8 mapped_angle_rad ! 当前角度的弧度值
integer lastKINC ! 用于跟踪上一个增量步
save lastKINC ! 保存lastKINC值
! ==== 初始化数组 ====
data angleData /"""
# 添加角度数据 (度)
# 第一行特殊处理(不加&)
for i in range(min(10, nPoints)):
fortran_code += f"{key_angles[i]:.6f}"
if i < min(9, nPoints - 1): # 不是行内最后一个且不是数组最后一个
fortran_code += ","
elif i < nPoints - 1: # 行内最后一个但不是数组最后一个
fortran_code += ",\n &," # 行末加逗号,新行以&,开头
# 后续行
for i in range(10, nPoints):
if i % 10 == 0: # 每10个元素换行
fortran_code += f"{key_angles[i]:.6f}" # 新行第一个元素
else:
fortran_code += f",{key_angles[i]:.6f}" # 行内元素加逗号前缀
# 行末处理(不是数组最后一个)
if (i + 1) % 10 == 0 and i < nPoints - 1:
fortran_code += ",\n &," # 行末加逗号,新行以&,开头
fortran_code += " /"
fortran_code += "\n data FzData /"
# 添加法向力数据
# 第一行特殊处理(不加&)
for i in range(min(10, nPoints)):
fortran_code += f"{key_fz[i]:.6f}"
if i < min(9, nPoints - 1): # 不是行内最后一个且不是数组最后一个
fortran_code += ","
elif i < nPoints - 1: # 行内最后一个但不是数组最后一个
fortran_code += ",\n &," # 行末加逗号,新行以&,开头
# 后续行
for i in range(10, nPoints):
if i % 10 == 0: # 每10个元素换行
fortran_code += f"{key_fz[i]:.6f}" # 新行第一个元素
else:
fortran_code += f",{key_fz[i]:.6f}" # 行内元素加逗号前缀
# 行末处理(不是数组最后一个)
if (i + 1) % 10 == 0 and i < nPoints - 1:
fortran_code += ",\n &," # 行末加逗号,新行以&,开头
fortran_code += " /"
fortran_code += "\n ! ==== 结束嵌入数据 ====\n\n"
fortran_code += f""" ! ==== 初始化状态标志 ====
data firstCall /.true./
data lastKINC /0/ ! 初始化lastKINC为0
! ==== 参数设置 ====
x0 = -457.497d-3
y0 = 87.503d-3
z0 = 0.0d0
b = 8.7036d-3
a = 5.4657d-3
R = 457.497d-3
! ==== 蠕滑区尺寸 ====
a1 = 3.1851d-3 ! 蠕滑区y方向半轴长
b1 = 4.8659d-3 ! 蠕滑区z方向半轴长
! ==== 蠕滑区中心 ====
xn0 = x0
yn0 = y0
zn0 = z0 - (a - a1) ! 蠕滑区中心位置
! ==== 载荷速度 ====
velocity = 400.0d0 / 3.6d0
omega = velocity / R
! ==== 获取当前时间 ====
t = TIME(1)
! ==== 直接使用KINC作为索引 ====
idx = KINC
! ==== 检查索引是否在有效范围内 ====
if (idx < 1 .or. idx > nPoints) then
write(*,*) '错误: 索引超出范围! KINC=', idx, '有效范围:1-', nPoints
F = 0.0d0
return
endif
! ==== 获取当前载荷点角度和载荷值 ====
current_angle_deg = angleData(idx)
Fz = FzData(idx)
! ==== 输出调试信息 ====
if (firstCall) then
firstCall = .false.
write(*,*) '直接使用KINC索引载荷模式'
write(*,*) '数据点数: ', nPoints
write(*,*) '起始角度: ', angleData(1), ' Fz: ', FzData(1)
write(*,*) '结束角度: ', angleData(nPoints), ' Fz: ', FzData(nPoints)
endif
! ==== 输出当前信息(每个增量步只输出一次) ====
if (KINC .ne. lastKINC) then
write(*,*) '增量步号: ', KINC,
1 ' 角度: ', current_angle_deg, '°'
write(*,*) 'Fz = ', Fz, 'N'
lastKINC = KINC
endif
! ==== 计算接触应力 ====
Tforce = 3.0d0 * Fz / (2.0d0 * pi * a * b)
! ==== 检查点是否在接触斑内 ====
x = COORDS(1)
y = COORDS(2)
z = COORDS(3)
! 坐标转换(使用当前数据点角度 - 转换为弧度)
mapped_angle_rad = current_angle_deg * pi / 180.0d0
xt = x * dcos(mapped_angle_rad) + z * dsin(mapped_angle_rad)
yt = y
zt = -x * dsin(mapped_angle_rad) + z * dcos(mapped_angle_rad)
! 检查点是否在椭圆内
check = ((yt - y0) / a)**2 + ((zt - z0) / b)**2
! 施加载荷或零
if (check .le. 1.0d0 .and. dabs(xt - x0) .le. 0.1d0 * R) then
F = Tforce * dsqrt(1.0d0 - check)
else
F = 0.0d0
endif
RETURN
END
C
C ======================================================================
C
SUBROUTINE UTRACLOAD(ALPHA,T_USER,KSTEP,KINC,TIME,NOEL,NPT,
1 COORDS,DIRCOS,JLTYP,SNAME)
C
INCLUDE 'ABA_PARAM.INC'
C
DIMENSION T_USER(3), TIME(2), COORDS(3), DIRCOS(3,3)
CHARACTER*80 SNAME
*******************头部文件************************
! ==== 声明所有变量 ====
integer nPoints
parameter (nPoints = {nPoints})
real*8 angleData(nPoints) ! 存储角度(度)
real*8 FzData(nPoints) ! 法向力
real*8 Fz, Tforce, omega
real*8 x0,y0,z0,check,x,y,z,xt,yt,zt
real*8 R, b, a, velocity
real*8 pi
parameter (pi = 3.141592653589793d0)
integer idx ! 直接使用KINC作为索引
logical firstCall
save firstCall
real*8 t
real*8 current_angle_deg ! 当前数据点角度(度)
real*8 mapped_angle_rad ! 当前角度的弧度值
integer lastKINC ! 用于跟踪上一个增量步
save lastKINC ! 保存lastKINC值
real*8 alpha1, E, v, fv, m, n, xp, yp, zp, check1
real*8 xn0, yn0, zn0, xck, F_hua, temp1, temp2
real*8 a1, b1 ! 蠕滑区尺寸
! ==== 初始化数组 ====
data angleData /"""
# 添加角度数据 (度)
# 第一行特殊处理(不加&)
for i in range(min(10, nPoints)):
fortran_code += f"{key_angles[i]:.6f}"
if i < min(9, nPoints - 1): # 不是行内最后一个且不是数组最后一个
fortran_code += ","
elif i < nPoints - 1: # 行内最后一个但不是数组最后一个
fortran_code += ",\n &," # 行末加逗号,新行以&,开头
# 后续行
for i in range(10, nPoints):
if i % 10 == 0: # 每10个元素换行
fortran_code += f"{key_angles[i]:.6f}" # 新行第一个元素
else:
fortran_code += f",{key_angles[i]:.6f}" # 行内元素加逗号前缀
# 行末处理(不是数组最后一个)
if (i + 1) % 10 == 0 and i < nPoints - 1:
fortran_code += ",\n &," # 行末加逗号,新行以&,开头
fortran_code += " /"
fortran_code += "\n data FzData /"
# 添加法向力数据
# 第一行特殊处理(不加&)
for i in range(min(10, nPoints)):
fortran_code += f"{key_fz[i]:.6f}"
if i < min(9, nPoints - 1): # 不是行内最后一个且不是数组最后一个
fortran_code += ","
elif i < nPoints - 1: # 行内最后一个但不是数组最后一个
fortran_code += ",\n &," # 行末加逗号,新行以&,开头
# 后续行
for i in range(10, nPoints):
if i % 10 == 0: # 每10个元素换行
fortran_code += f"{key_fz[i]:.6f}" # 新行第一个元素
else:
fortran_code += f",{key_fz[i]:.6f}" # 行内元素加逗号前缀
# 行末处理(不是数组最后一个)
if (i + 1) % 10 == 0 and i < nPoints - 1:
fortran_code += ",\n &," # 行末加逗号,新行以&,开头
fortran_code += " /"
fortran_code += "\n ! ==== 结束嵌入数据 ====\n\n"
fortran_code += f""" ! ==== 初始化状态标志 ====
data firstCall /.true./
data lastKINC /0/ ! 初始化lastKINC为0
! ==== 参数设置 ====
x0 = -457.497d-3
y0 = 87.503d-3
z0 = 0.0d0
b = 8.7036d-3
a = 5.4657d-3
R = 457.497d-3
! ==== 蠕滑区尺寸 ====
a1 = 3.1851d-3 ! 蠕滑区y方向半轴长
b1 = 4.8659d-3 ! 蠕滑区z方向半轴长
! ==== 蠕滑区中心 ====
xn0 = x0
yn0 = y0
zn0 = z0 - (a - a1) ! 蠕滑区中心位置
! ==== 载荷速度 ====
velocity = 400.0d0 / 3.6d0
omega = velocity / R
! ==== 获取当前时间 ====
t = TIME(1)
! ==== 直接使用KINC作为索引 ====
idx = KINC
! ==== 检查索引是否在有效范围内 ====
if (idx < 1 .or. idx > nPoints) then
write(*,*) '错误: 索引超出范围! KINC=', idx, '有效范围:1-', nPoints
ALPHA = 0.0d0
return
endif
! ==== 获取当前载荷点角度和载荷值 ====
current_angle_deg = angleData(idx)
Fz = FzData(idx)
! ==== 输出调试信息 ====
if (firstCall) then
firstCall = .false.
write(*,*) 'UTRACLOAD: 直接使用KINC索引载荷模式'
write(*,*) '数据点数: ', nPoints
write(*,*) '起始角度: ', angleData(1), ' Fz: ', FzData(1)
write(*,*) '结束角度: ', angleData(nPoints), ' Fz: ', FzData(nPoints)
endif
! ==== 输出当前信息(每个增量步只输出一次) ====
if (KINC .ne. lastKINC) then
write(*,*) 'UTRACLOAD: 增量步号: ', KINC,
1 ' 角度: ', current_angle_deg, '°'
write(*,*) 'Fz = ', Fz, 'N'
lastKINC = KINC
endif
! ==== 计算接触应力 ====
Tforce = 3.0d0 * Fz / (2.0d0 * pi * a * b)
! ==== 蠕滑区参数 ====
E = 210e9 ! 弹性模量
v = 0.29 ! 泊松比
fv = 0.4 ! 摩擦系数
! ==== 获取当前点坐标 ====
x = COORDS(1)
y = COORDS(2)
z = COORDS(3)
! ==== 坐标转换 ====
mapped_angle_rad = current_angle_deg * pi / 180.0d0
alpha1 = mapped_angle_rad ! 使用当前角度作为alpha1
xt = x * dcos(alpha1) + z * dsin(alpha1)
yt = y
zt = -x * dsin(alpha1) + z * dcos(alpha1)
! ==== 检查点是否在接触斑内 ====
check = ((yt - y0) / a)**2 + ((zt - z0) / b)**2
check1 = ((yt - yn0) / a1)**2 + ((zt - zn0) / b1)**2
xck = dabs(xt - x0)
! ==== 计算切向力大小 ====
F_hua = fv * Tforce * dsqrt(1.0d0 - check) ! 滑动区切向力
! ==== 初始化切向力 ====
ALPHA = 0.0d0
! ==== 检查点位置并计算切向力 ====
if (xck .le. R / 9.0d0) then
! 滑动区切向力
if (check .le. 1.0d0 .and. check1 .gt. 1.0d0) then
ALPHA = F_hua
endif
! 蠕滑区切向力(衰减)
if (check1 .le. 1.0d0) then
temp1 = ((xt - a + a1) / a1)**2 + (yt / b1)**2
temp2 = a1 / a * dsqrt(1.0d0 - temp1)
ALPHA = F_hua - temp2
endif
endif
! ==== 设定切向力方向 ====
T_USER(1) = dsin(alpha1)
T_USER(2) = 0.0d0
T_USER(3) = -dcos(alpha1)
RETURN
END
"""
# 写入Fortran文件
with open(output_for, 'w') as f:
f.write(fortran_code)
print(f"\n成功生成包含DLOAD和UTRACLOAD的Fortran子程序: {output_for}")
print(f"关键角度数据点数: {nPoints}")
print(f"角度范围: {min(key_angles):.2f}° 到 {max(key_angles):.2f}°")
print(f"法向力范围: {min(key_fz):.6f}N 到 {max(key_fz):.6f}N")
return True
except Exception as e:
print(f"\n错误: {str(e)}")
import traceback
traceback.print_exc()
return False
if __name__ == "__main__":
# 输入文件路径
input_txt = r"D:\Desktop\FEIWENTAI\波磨结果\output_V400_u-0.05已修正_W0.08\Normalforce4-Longitudinalforce1-creepage右轮.txt"
# 输出文件路径 - 修改为 DONG_all1.for
output_for = r"D:\Desktop\DONG_all1.for"
print("=" * 50)
print("数据对比工具 + FORTRAN子程序生成器")
print("=" * 50)
print(f"输入文件: {input_txt}")
print(f"输出文件: {output_for}")
print("=" * 50)
# 执行转换
success = generate_final_fortran_subroutine(input_txt, output_for)
if success:
print("\n" + "=" * 50)
print("数据对比完成!FORTRAN代码已成功生成")
print("=" * 50)
else:
print("\n" + "=" * 50)
print("处理失败,请检查错误信息")
print("=" * 50)
最新发布