Parallel network models

[CF]

Hello all
I am currently working with a constitutive model that consists of two parallel networks acting in unison- something like the Bergstrom-Boyce model. The total stress is equal to the sum of the stresses in each network i.e. S_tot = S_A + S_B.
I am running simple one-element uniaxial compression analyses (displacement control) and I've noticed something, which I find peculiar. When I change the material parameters governing Network B the stress in Network A (S_A) is different. When I examine it further I see the deformation gradient components F(2,2) and F(3,3) are different each time. Should this be the case? When I solve the problem numerically in Matlab this does not happen.
I've tried it with a modified version of the vumat_nh.for that Jorgen posted a while back and I get the same thing happening. I've pasted it below. Perhaps I'm missing something obvious :? but if anyone can point me the right way I'd be much obliged.

Best regards,
Cormac

c
c Copyright 2004 Jorgen Bergstrom
c
subroutine vumat(nblock, ndi, nshr, nstatev, nfieldv, nprops,
. lanneal, stepTime, totTime, dt, cmname, coordMp,
. charLen, props, density, Dstrain, rSpinInc, temp0,
. U0, F0, field0, stressVec0, state0,
. intEne0, inelaEn0, temp1, U1,
. F1, field1, stressVec1, state1, intEne1, inelaEn1)
implicit none
integer nblock, ndi, nshr, nstatev, nfieldv, nprops, lanneal
real stepTime, totTime, dt
character*80 cmname
real coordMp(nblock,*)
integer charLen
real props(nprops), density(nblock),
. Dstrain(nblock,ndi+nshr), rSpinInc(nblock,nshr),
. temp0(nblock), U0(nblock,ndi+nshr),
. F0(nblock,ndi+nshr+nshr), field0(nblock,nfieldv),
. stressVec0(nblock,ndi+nshr), state0(nblock,nstatev),
. intEne0(nblock), inelaEn0(nblock), temp1(nblock),
. U1(nblock,ndi+nshr), F1(nblock,ndi+nshr+nshr),
. field1(nblock,nfieldv), stressVec1(nblock,ndi+nshr),
. state1(nblock,nstatev), intEne1(nblock), inelaEn1(nblock),
. stressVecA(nblock,ndi+nshr),stressVecB(nblock,ndi+ nshr)
c local variables
real F(3,3), B(3,3), J, t1, t2, t3, mu, mu2, kappa,dfgrd1(3,3),
&Finverse(3,3),rotation(3,3)
integer i, x

mu = 5.0
kappa = 100.0
mu2 = 10.0
c loop through all blocks
do i = 1, nblock
c setup F (upper diagonal part)
F(1,1) = U1(i,1)
F(2,2) = U1(i,2)
F(3,3) = U1(i,3)
F(1,2) = U1(i,4)
if (nshr .eq. 1) then
F(2,3) = 0.0
F(1,3) = 0.0
end if
!
c calculate J
t1 = F(1,1) * (F(2,2)*F(3,3) - F(2,3)**2)
t2 = F(1,2) * (F(2,3)*F(1,3) - F(1,2)*F(3,3))
t3 = F(1,3) * (F(1,2)*F(2,3) - F(2,2)*F(1,3))
J = t1 + t2 + t3
t1 = J**(-2.0/3.0)
!
c Bstar = J^(-2/3) F Ft (upper diagonal part)
B(1,1) = t1*(F(1,1)**2 + F(1,2)**2 + F(1,3)**2)
B(2,2) = t1*(F(1,2)**2 + F(2,2)**2 + F(2,3)**2)
B(3,3) = t1*(F(1,3)**2 + F(2,3)**2 + F(3,3)**2)
B(1,2) = t1*(F(1,1)*F(1,2)+F(1,2)*F(2,2)+F(1,3)*F(2,3))
B(2,3) = t1*(F(1,2)*F(1,3)+F(2,2)*F(2,3)+F(2,3)*F(3,3))
B(1,3) = t1*(F(1,1)*F(1,3)+F(1,2)*F(2,3)+F(1,3)*F(3,3))
!
c StressA = mu/J * Dev(Bstar) + kappa*log(J)/J * Eye
t1 = (B(1,1) + B(2,2) + B(3,3)) / 3.0
t2 = mu/J
t3 = kappa * log(J) / J
StressVecA(i,1) = t2 * (B(1,1) - t1) + t3
StressVecA(i,2) = t2 * (B(2,2) - t1) + t3
StressVecA(i,3) = t2 * (B(3,3) - t1) + t3
StressVecA(i,4) = t2 * B(1,2)
if (nshr .eq. 3) then
StressVecA(i,5) = t2 * B(2,3)
StressVecA(i,6) = t2 * B(1,3)
end if

c StressB = mu/J * Dev(Bstar) + kappa*log(J)/J * Eye
t1 = (B(1,1) + B(2,2) + B(3,3)) / 3.0
t2 = mu2/J
t3 = kappa * log(J) / J
StressVecB(i,1) = t2 * (B(1,1) - t1) + t3
StressVecB(i,2) = t2 * (B(2,2) - t1) + t3
StressVecB(i,3) = t2 * (B(3,3) - t1) + t3
StressVecB(i,4) = t2 * B(1,2)
if (nshr .eq. 3) then
StressVecB(i,5) = t2 * B(2,3)
StressVecB(i,6) = t2 * B(1,3)
end if

! Total stress
StressVec1 = StressVecA + StressVecB

! Spit out to sta file
write(6,*) stressVec1(1,1), stressVecA(1,1), stressVecB(1,1)

end do
return
end

[Jorgen]

Hi Cormac,

I have not yet tested your subroutine, but I have a few comments to your questions.

:arrow: I am not suprised that the stress in Network A changes when you change the properties of Network B. The reason for that is that in uniaxial compression the transverse stresses (S_tot_22 and S_tot_33) are always zero. But the stress components S_tot_22 and S_tot_33 are obtained by the sum of S_A and S_B, and if you change the properties of Network B then the stress S_B will change and as a results the stress S_A will also need to change in order to fulfill the zero stress condition in the transverse direction.

How does that sound?
- Jorgen

[CF]

Thanks for the quick response Jorgen. I guess that sounds ok to me. My initial intuition made me think that S_A(2,2), S_A(3,3), S_B(2,2), S_B(3,3) should all be zero in uniaxial compression.
I've updated my own solver and all results now agree.

Thanks
Cormac