Vdload子程序的理解
vdload子程序关于坐标问题
Colin 发表于: 五月7th, 2011 | 阅读: 893 次阅读| 评论数: (2)
vdload是abaqus显式分析（explicit）的一个子程序。
官方的帮助文档：（博主译）
用户子程序VDLOAD:
可以用来定义一组点关于位置, 时间, 速度等的函数的分布载荷变量，每一个点出现在an element-based or
surface-based nonuniform load definition;
will be called for load integration points associated with each nonuniform load definition including PENU and PINU loads applicable for pipe elements;
does not make available the current value of the nonuniform distributed loads for file output purposes; and recognizes an amplitude reference (“Amplitude curves,” Section 30.1.2 of the Abaqus Analysis User’s Manual) if it appears with the associated nonuniform load definition.
User subroutine interface（用户子程序接口）
subroutine vdload (
C Read only (unmodifiable)variables -
1 nblock, ndim, stepTime, totalTime,
2 amplitude, curCoords, velocity, dirCos, jltyp, sname,
C Write only (modifiable) variable -
1 value )
C
include ‘vaba_param.inc’
C
dimension curCoords(nblock,ndim), velocity(nblock,ndim),
1 dirCos(nblock,ndim,ndim), value(nblock)
character*80 sname
C
do 100 km = 1, nblock
user coding to define value
100 continue
return
end
需要定义的变量：
value (nblock)
均布载荷的大小。面载荷的单位为FL–2 ，体积力的单位为FL–3。
传入信息的变量
nblock
Number of points to be processed in this call to VDLOAD.
ndim
坐标方向的数值：二维模型值为2，三维模型值为3。若有任意一个三维单元则模型被认为时三维的。(包括SPRINGA 单元).
stepTime
Value of time since the step began.
totalTime
Value of total time. The time at the beginning of the step is given by totalTime – stepTime.
amplitude
Current value of the amplitude referenced for this load (set to unity if no amplitude is referenced). You must multiply the load by the current amplitude value within the user subroutine if the amplitude is required.
curCoords (nblock, ndim)
Current coordinates of each point for which the load is to be calculated.
velocity (nblock, ndim)
Current velocity of each point for which the load is to be calculated.
dirCos (nblock, ndim, ndim)
Current orientation of the face, edge, pipe, or beam for pressure type loads (not applicable for body force type loads). The second dimension indicates the vector, and the third dimension indicates the components of that vector. For faces (pressures on three-dimensional continuum, shell, and membrane elements), the first and second vectors are the local directions in the plane of the surface and the third vector is the normal to the face, as defined in “Conventions,” Section 1.2.2 of the Abaqus Analysis User’s Manual. For solid elements the normal points inward, which is the opposite of what is defined in the conventions; for shell elements the normal definition is consistent with the defined conventions. For edges (pressures on two-dimensional continuum elements and two-dimensional beams and pipes), the first vector is the normal to the edge, the second vector is the tangent to the edge, and, if ndim=3, the third vector will be a unit normal in the out-of-plane direction. For three-dimensional beam and pipe elements, the first and second vectors are the local axes (, ) and the third vector is the tangent vector (), as defined in “Beam element cross-section orientation,” Section 26.3.4 of the Abaqus Analysis User’s Manual.
jltyp
Key that identifies the distributed load type. The load type may be a body force, a surface-based load, or an element-based surface load. For element-based surface loads, this variable identifies the element face for which this call to VDLOAD is being made. See Part VI, “Elements,” of the Abaqus Analysis User’s Manual for element load type identification. This information is useful when several different nonuniform distributed loads are being imposed on an element at the same time.
sname
Surface name for a surface-based load definition (JLTYP=0). For a body force or an element-based load the surface name is passed in as a blank.
以前一直纠结于curCoords (nblock, ndim)的参数，据文档值nblock代表调用VDLOAD子程序的点的个数，ndim代表着维数。2维模型ndim=2,三维模型ndim=3。我一直在想函数的返回值是什么呢，是一个数组？还是一个具体数，含义是什么呢？
最近终于知道我的理解大错特错了，ndim在3维情况下取值为1~3，当ndim=1是，返回值为1方向坐标；ndim=2，返回值为
2方向作品。ndim=3，返回值为3方向作品。而ndim在2维情况下取值为1~2。这样，我们就可以很轻松的调用这个函数来表示坐标了。
2个例子
C User subroutine VDLOAD
subroutine vdload (
C Read only (unmodifiable) variables -
* nblock, ndim, stepTime, totalTime,
* amplitude, curCoords, velocity, dircos,
* jltyp, sname,
C Write only (modifiable) variable -
* value )
C
include 'vaba_param.inc'
parameter (
* d0 = 120.00d0,
* pa = 0.0d0,
* ps = 364.0d0,
* theta = 0.124651928275d-3,
c * theta = 0.1569262d-3,
* timea = 0.01875d-3,
* timed = 0.195d-3,
* zero = 0.0d0,
* one = 1.0d0
* )
C
dimension curCoords(nblock,ndim),
* velocity(nblock,ndim),
* dircos(nblock,ndim,ndim),
* value(nblock)
character*80 sname
c
c local variables
c
real xx,yy,zz,dd,dd0,xc,yc,zc,tt,Tint,press
c
xc = zero
yc = zero
zc = zero
time = totalTime
dd0 = d0 * d0
Tint = timed-timea
tt = time - timea
do k = 1, nblock
press = zero
if ((time .gt. timea) .and. (time .le. timed)) then
xx = curCoords(k,1) - xc
yy = curCoords(k,2) - yc
zz = curCoords(k,3) - zc
xx = xx * xx
yy = yy * yy
zz = zz * zz
dd = xx + yy + zz
press = (ps-pa) * (one - (tt/Tint))
press = press * (exp(-tt/theta))
press = press * (exp(-dd/dd0))
end if
value(k) = amplitude * press
c write(*,*) time, value(k)
end do
*
return
end
c
c=================================================== subroutine vdload (
C Read only (unmodifiable)variables -
1 nblock, ndim, stepTime, totalTime,
2 amplitude, curCoords, velocity, dirCos, jltyp, sname,
C Write only (modifiable) variable -
1 value )
C
include 'vaba_param.inc'
C
dimension curCoords(nblock,ndim), velocity(nblock,ndim),
1 dirCos(nblock,ndim,ndim), value(nblock)
character*80 sname
C
parameter(vel=5.0d0,dlen=0.1d0,pressure=500000.0d0)
integer k
disp=vel*steptime
zc=disp
zmax=zc
zmin=zmax-dlen
zmax1=zmax+0.5
zmin1=zmin+0.5
zmax2=zmax+1.5
zmin2=zmin+1.5
zmax3=zmax+2
zmin3=zmin+2
IF (SNAME(1:15).EQ.'ASSEMBL Y_SURF-1') THEN
do k=1, nblock
if(curCoords(k,3).lt.zmax.and.curCoords(k,3).ge.zmin) then value(k)=amplitude * pressure
else if(curCoords(k,3).lt.zmax1.and.curCoords(k,3).ge.zmin1) then value(k)=amplitude * pressure
else if(curCoords(k,3).lt.zmax2.and.curCoords(k,3).ge.zmin2) then value(k)=amplitude * pressure
else if(curCoords(k,3).lt.zmax3.and.curCoords(k,3).ge.zmin3) then value(k)=amplitude * pressure
else
value(k)=0.0
end if
enddo
ENDIF
return
end