ABAQUS中椭圆形移动载荷DLOAD和UTRACLOAD子程序详解：从定义到实现

   在有限元模拟中，重复移动载荷（Repeated moving pressure）是结构受力分析中用于等效模拟接触载荷的一个重要手段，尤其在轮轨接触、滚珠接触、焊接热源移动等问题研究中极为常见。本文主要介绍ABAQUS中椭圆形移动载荷定义、法向和切向载荷模拟、子程序DLOAD和UTRACLOAD编程实现，实现建议与注意事项。

1、椭圆形移动载荷定义

  移动载荷指的是随时间或空间位置变化而不断变化施加位置的载荷，其典型例子包括：1）行驶车辆对桥梁的作用力；2）火车车轮与轨道之间的接触力；3）滚动体在接触面上滑移产生的局部接触载荷；4）焊接过程中热源的沿路径移动。这些载荷不是固定不动的，而是随时间在接触体上“移动”，从而引发结构响应的动态变化。在应力应变分析、疲劳寿命评估等方面，考虑载荷的移动性尤为关键。

  在滚动体的接触中，Hertz型椭圆形接触斑较为常见，其形状可根据Hertz接触理论表示为：

其中，P为总法向力，a和b分别为横向x和纵向z上的接触斑半宽，p0为最大接触压力。

2、法向和切向移动载荷模拟

  在ABAQUS中，模拟移动载荷的两种典型方法分别对应法向载荷和切向载荷。

2.1 法向移动载荷

  法向载荷定义见式(1)所示。在给定总法向力P或者轴重，以及接触斑长半轴和短半轴大小后，即可确定出来p(x,z)空间分布。其中，P、a和b可以通过Hertz接触理论或者有限元法计算得到，也可以通过一些网站去快速计算，比如：https://www.tribology-abc.com/sub10.htm以及https://www.pecms.cn/hz/hzb2p。

图1 法向接触压力

2.2 切向移动载荷

  在滚动接触过程中，除了接触表面的法向接触压力外，接触体还存在局部滑动或者蠕滑，导致接触斑区域被划分为黏着区和滑动区。其中，沿着滚动方向的后沿为滑动区，前沿则为黏着区。在滑动区内，切向力大小为摩擦系数乘以法向力，而在黏着区内切向力大小与蠕滑率有关。

  为了计算接触斑内的黏滑分布，首先采用Haines和Ollerton条带理论将接触斑划分为若干个纵向条带，对每个条带使用Carter接触理论（图3）划分黏着区和滑动区大小，进而确定切向力大小分布。

图2 条带理论

图3 利用Crter接触理论划分黏滑分布

    1）在滑动区内，切向力大小可表示为：

其中，f为摩擦系数。

    2）在黏着区内，切向力分布由q1和q2两部分组成，合起来可表示为：

其中，d为黏着区与滑动区中心之间的距离，可表示为d=a−c，c为黏着区半宽。这里，c与黏着区黏着因子K具有如下关系：

其中，K为黏着区黏着因子，而K与蠕滑率ξ具有如下关系：

  K与蠕滑率ξ之间的关系曲线如图4：

图4  K与ξ关系曲线

  因此，切向接触载荷大小与最大接触压力p0、接触斑形状a和b、蠕滑率ξ、摩擦系数f密切相关。如果给定一个蠕滑率ξ，那么就可以依次确定出K、黏着区大小c、切向力大小q分布（如图5）。

图5 切向载荷分布

3、子程序DLOAD和UTRACLOAD编程实现

  确定好法向和切向移动载荷分布后，利用FORTRAN编写DLOAD和UTRACLOAD用户子程序，实现法向和切向移动载荷施加，载荷施加在滚动接触体表面（以图6所示的钢轨踏面施加移动载荷为例）。

图6 钢轨表面施加法向和切向移动载荷

3.1 法向移动载荷实现DLOAD

  DLOAD子程序是ABAQUS中定义体载荷、面载荷、线载荷等的一种接口，通过Fortran代码自定义每个积分点上的载荷值。DLOAD适用于定义在单元上的载荷（如压力、密度效应等）；可以利用时间（TIME(1)）、空间坐标（COORD）、元素编号等信息，来定义移动的载荷区域或强度，主要用于法向载荷的模拟。

   部分代码示意如下：

******************************* Top **********************************       	
************************* Dload 法向移动载荷 *************************		  
	SUBROUTINE DLOAD(F, KSTEP, KINC, TIME, NOEL, NPT, LAYER, KSPT,
     &           COORDS, JLTYP, SNAME)
	INCLUDE 'ABA_PARAM.INC'
	DIMENSION TIME(2), COORDS(3)
	CHARACTER*80 SNAME
	real*8 x, y, z, Velocity, Move, centerx, para, P0, A, B, P0Y, AA
*     接触斑参数
   P0 = 1258.0D0 ! 最大接触压力
	A = 7.3D0   ! 长半轴
	B = 5.5D0   ! 短半轴
*     移动载荷参数
	num_cycles = 100  !循环100次
	StepTime = 0.0072  !每个循环时间  t=800mm/27777.78mm/s=0.0288
	FixedMove = 200.0D0
	Velocity = FixedMove / StepTime   ! v=100*1000/3.6=27777.78mm/s
	TotalCycles = 100      ! 总循环次数
*      坐标
	x0 = -100.0D0
	y0 = 2.78D0
	x = COORDS(1)
	y = COORDS(2)
	z = COORDS(3)
    
*      确定载荷中心位置
	IF (TIME(1) .LE. StepTime * TotalCycles) THEN
		tmod = MOD(TIME(1), StepTime)   ! 当前时间在单次循环中的时间
		Move = Velocity * tmod          ! 移动距离
		centerx = x0 + Move
	ELSE
		centerx = x0    ! 超过 100 次循环后，载荷静止
	END IF	


*      判断椭圆接触斑内
	para = (x - centerx)**2 / A**2 + y**2 / B**2   

3.2 切向移动载荷实现UTRACLOAD

  UTRACLOAD子程序用于定义接触面上的切向载荷（traction load），尤其适合模拟轮轨滚动、摩擦拖曳、履带驱动等问题。与接触对配合使用；可通过位置（COORD）、法向和切向方向、时间等参数动态指定接触面上移动的切向载荷；适合模拟摩擦力、驱动力等切向移动载荷。

   部分代码示意如下：

************************* UtracLoad 切向移动载荷 ***************	
	SUBROUTINE UTRACLOAD(ALPHA,T_USER,KSTEP,KINC,TIME,NOEL,NPT,
     &		   COORDS,DIRCOS,JLTYP,SNAME)
	INCLUDE 'ABA_PARAM.INC'
	DIMENSION T_USER(3), TIME(2), COORDS(3), DIRCOS(3,3)
	CHARACTER*80 SNAME
	real*8 x,y,z,centerx,para,KI,CC,StepTime,TotalCycles,Velocity
	real*8 P0,A,B,FRI,AA,P0X,tmod,ALP1


*      参数初始化
	P0 = 1258.0D0
	A  = 7.3D0
	B  = 5.5D0
	FRI = 0.5D0
	KI = 0.8D0  !蠕滑分布系数，并非牵引系数，K=1，纯滚；K=0，全滑
*      移动载荷参数
	num_cycles = 100  !循环100次
	StepTime = 0.0072  !每个循环时间
	FixedMove = 200.0D0
	Velocity = FixedMove / StepTime 
	TotalCycles = 100      ! 总循环次数
    
	x=COORDS(1)
	y=COORDS(2)
	z=COORDS(3)
*      移动载荷的坐标
	x0 = -100.0D0
	y0 = 2.78D0
	T_USER(1)=1.0
	T_USER(2)=0.0
	T_USER(3)=0.0

4、实现建议与注意事项

  路径控制：若载荷沿非直线路径移动，可通过函数表或外部数据点定义轨迹；

  速度控制：建议从较低速度开始测试，确保模拟稳定；

  边界条件耦合：对于周期性或滚动加载问题，可结合对称/循环边界实现持续滚动效果；

  输出检查：使用可视化工具（如 contour 图或 field 输出）确认载荷移动过程是否正确。

如果你正在处理类似问题，欢迎留言交流，后文也分享了一个成功测试过的子程序代码！