|
|
朋友们，我在计算子程序时发现，直接运行显示NAN，但是Cnf是有计算出来的值的；如果将Cnf直接设为计算出来的值，计算没有问题，下面是我的核心程序，UMAT开头格式已省略，全文源代码网址Fcode Pastebin
integer i, j, k, l
real(KIND=8) ZERO, ONE, TWO, THREE, A0, n, q1, q2, v, ABFF,ABFF1,
1 pivot, Cnf, k0, q4, Smid, S0, DTIME1, k1, A, E1, Ac, Bc, kapa0,
2 kapa, AA, choose1, choose2, choose3, AAA, dtime2, AAAA, Smid1,
3 SS, SSS, SSSS, SSSSS, ssssss, SA, Smid2,sssssss,ACA,ABFF2,
4 Cnv, Taomiu, Amiu, Cnv1,abff4,abff5,timenext,alphal,alphar,
5 alphamid,na,Cnv2, VV, aging,elastic, GC, Cnf1
parameter(ZERO=0.0d0, TWO=2.0d0, ONE=1.0d0, ACA=1.0d1, n=1.0d-1,
1 A0=0.0d0, THREE=3.0d0, na=1.9d0)
real(kind=8), dimension(ntens) :: dstrainv
real(kind=8), dimension(ntens) :: dstrainf
real(kind=8), dimension(ntens) :: gamamiu
real(kind=8), dimension(ntens,ntens) :: Iddsdde
real(kind=8), dimension(ntens,ntens) :: IIddsdde
real(kind=8), dimension(ntens,2*ntens) :: aug
real(kind=8), dimension(ntens) :: DSTRESS1
real(kind=8), dimension(ntens) :: DSTRESS2
real(kind=8), dimension(ntens) :: STRESS2
real(kind=8), dimension(ntens) :: STRESS1
real(kind=8), dimension(2*ntens) :: ABC
real(kind=8), dimension(2*ntens) :: ABD
real(kind=8), dimension(2*ntens) :: ABE
real(kind=8), dimension(2*ntens) :: ABF
real(kind=8), dimension(ntens,ntens) :: Dve1
real(kind=8), dimension(ntens,ntens) :: DD
real(kind=8), dimension(ntens,ntens) :: II
real(kind=8), dimension(ntens,2*ntens) :: aug1
real(kind=8), dimension(ntens,ntens) :: Dve
real(kind=8), dimension(ntens) :: PS
real(kind=8), dimension(ntens) :: SB
LOGICAL flok
C material properties
q1 = props(1)
v = props(2)
q2 = props(3)
Ac = props(4)
Bc = props(5)
k0 = props(6)
q4 = props(7)
S0 = props(8) * 1.0d0
kapa0 = props(9)
E1 = props(10)
C
if (KINC.EQ.ONE) then
DTIME1 = ZERO
do i = 1, ntens
DSTRESS1(i) = ZERO
STRESS1(i) = ZERO
end do
else
DTIME1 = STATEV(1) / (2.4d1 * 3.6d3)
Smid = STATEV(14)
do i = 1, ntens
k1 = i + 1.0d0
STRESS1(i) = STATEV(k1)
end do
do i = 1, ntens
k1 = i + 7.0d0
DSTRESS1(i) = STATEV(k1)
end do
end if
C G matrix
do i = 1, ntens
do j = 1, ntens
Iddsdde(i,j) = 0.0d0
end do
end do
do i = 1, ndi
do j = 1, ndi
Iddsdde(i,j) = -v
end do
Iddsdde(i,i) = 1.0d0
end do
do i = ndi+1, ntens
Iddsdde(i,i) = TWO * (ONE + v)
end do
dtime2 = dtime / (2.4d1 * 3.6d3)
timenext = time(1) / (2.4d1 * 3.6d3) + 5.0d-1 * dtime2
do i = 1,ntens
STRESS2(i) = stress(i)
DSTRESS2(i) = DSTRESS1(i)
end do
C viscous
choose2 = 1
if (choose2.EQ.1) then
do i = 1, ntens
dstrainf(i) = 0.0d0
end do
Cnf = ZERO
if (KINC.EQ.ONE) then
Smid = S0
else
Smid = Smid
end if
Smid1 = ((1.0d0 + 2.0d0 * k0 * dtime / 2.4d1 / 3.6d3 * Smid)
1 ** 5d-1 - 1.0d0) / (k0 * dtime / 2.4d1 / 3.6d3)
Smid = Smid1
Cnf1 = 5.0d-1 * q4 * k0 * Smid * dtime2
do i =1, ntens
A = 0.0d0
do j = 1, ntens
SB(i) = STRESS(i)
A = A + Iddsdde(i,j) * stress(j)
end do
dstrainf(i) = q4 * k0 * Smid * A * dtime2
end do
else
Cnf = 0
do i = 1,ntens
dstrainf(i) = 0
end do
end if
Cnf = Cnf1
C Damage
elastic = 0
choose3 = 0
if (choose3.EQ.1) then
AA = 0.0d0
CALL SPRINC(STRAN,PS,2,NDI,NSHR)
do i = 1, ndi
AA = AA + 1.0d-4
end do
kapa = kapa0 + AA ** 5.0d-1
d = 1.0d0 - kapa0 * (1.0d0 - Ac) / kapa - Ac /
1 (exp(Bc * (kapa - kapa0)))
else
d = 0
end if
C E
if (choose3.EQ.1) then
AAA = 1.0d0
else
AAA = 0.0d0
end if
if (elastic.EQ.1) then
q1 = q1
else
q1 = 0.0d0
end if
do i = 1, ntens
do j = 1, ntens
GC = Iddsdde(i,j) * (Cnv + Cnf + q1 + AAA / E1)
IIddsdde(i,j) = GC
end do
end do
do i = 1, ntens
do j = 1, ntens
Iddsdde(i,j) = IIddsdde(i,j)
end do
end do
C Dve1
if (choose3.EQ.1) then
do i = 1, ntens
do j = 1,ntens
II(i,j) = 0
end do
end do
do i = 1,ntens
II(i,i) = 1
end do
do i = 1,ntens
do j =1,ntens
Dve1(i,j) = ((1 - d) * E1 * II(i,j) + d *
1 Iddsdde(i,j)) / ((1 - d) * E1)
end do
end do
else
do i = 1,ntens
do j = 1,ntens
Dve1(i,j) = 0.0d0
end do
end do
end if
C Inversion of the matrix Iddsdde
aug(:,1:TWO*ntens) = 0.0d0
aug(:,1:ntens) = Iddsdde(:,1:ntens)
aug(:,ntens+1:TWO*ntens) = 0.0d0
do k = 1, ntens
aug(k,k+ntens) = 1.0d0
end do
do i = 1, ntens
pivot = aug(i,i)
do l = 1, 2*ntens
aug(i,l) = aug(i,l) / pivot
end do
do j = 1, ntens
if (j /= i) then
pivot = aug(j,i)
do k = 1, 2*ntens
aug(j,k) = aug(j,k) - pivot * aug(i,k)
end do
end if
end do
end do
DD(:,1:ntens) = aug(:,ntens+1:TWO*ntens)
C Inversion of the matrix Dve1
if (choose3.EQ.1) then
aug1(:,1:TWO*ntens) = 0.0d0
aug1(:,1:ntens) = Dve1(:,1:ntens)
aug1(:,ntens+1:TWO*ntens) = 0.0d0
do i = 1, ntens
aug1(i,i+ntens) = 1.0d0
end do
do i = 1, ntens
pivot = aug1(i,i)
do j = 1, 2*ntens
aug1(i,j) = aug1(i,j) / pivot
end do
do j = 1, ntens
if (j /= i) then
pivot = aug1(j,i)
do k = 1, 2*ntens
aug1(j,k) = aug1(j,k) - pivot * aug1(i,k)
end do
end if
end do
end do
Dve(:,1:ntens) = aug1(:,ntens+1:TWO*ntens)
else
end if
C Reall Dve
do l = 1,ntens
do k = 1,ntens
ddsdde(l,k) = 0
end do
end do
if (choose3.EQ.1) then
do l = 1,ntens
do j = 1,ntens
do k = 1,ntens
ddsdde(l,j) = ddsdde(l,j) + Dve(l,k) * DD(k,j)
end do
end do
end do
else
do l = 1,ntens
do j = 1,ntens
ddsdde(l,j) = ddsdde(l,j) + DD(l,j)
end do
end do
end if
C STATEV
k1 = 1
STATEV(k1) = dtime
do i = 1, ntens
k1 = i + 1
STATEV(k1) = stress(i)
end do
do i = 1, ntens
k1 = i + 7
do j = 1, ntens
STATEV(k1) = ddsdde(i,j) * (dstran(j) -
1 dstrainv(j) - dstrainf(j))
end do
end do
STATEV(14) = Smid
C Stress increment evaluation
do l = 1, ntens
do k = 1, ntens
stress(l) = stress(l) + ddsdde(l,k) * (dstran(k) -
1 dstrainv(k) - dstrainf(k))
end do
end do
|
|
F 币
贡献
权杖
窥视卡
发表于 2024-7-11 16:50:48
微信
收藏
点赞
点踩
提升卡
置顶卡
变色卡
千斤顶
显身卡
楼主