subroutine microsphere(ud,fb,theta,sig,dd) c-----[--.----+----.----+----.-----------------------------------------] c Purpose: User Constitutive Model: Microsphere model of rubber c elasticity (C. Miehe, S. Goektepe, F. Lulei, 2003) c Input: c fb(3,3) - Fbar: unimodular part of Deformation gradient c at point (finite deformation) c theta - Determinant of deformation gradient c ud(5) - User material parameters (nud) c Output: c sig(6) - cauchy stress c dd(6,6) - Current stiffness c Variables: c p ud(1) non-affine stretch parameter c N ud(2) number of chain segments c mu ud(3) shear modulus c q ud(4) non-affine tube parameter c U ud(5) tube geometry parameter c bb ud(6) bulk modulus c-----[--.----+----.----+----.-----------------------------------------] implicit none integer i,j,k,a,b integer co(6,2) real*8 theta,ud(6),fb(3,3) real*8 sig(6),dd(6,6) real*8 bb,r(21,4),t(3,21),lab(21),la real*8 h(6),hh(6,6),tauf,cf real*8 taufb(6),cfb(6,6),const,det,fbit(3,3) real*8 n(3,21),vb(21),kk(6),kkk(6,6),G(6,6) real*8 ij(3,3),taucb(6),ccb(6,6) real*8 taub(6),cb(6,6) real*8 p,bulk real*8 tmp1,tmp2,tmp3,tmp4,tmp5(6),tmp6(6),tmp7 save c conversion matrix (conversion to voigt) data co /1,2,3,1,2,1 ,1,2,3,2,3,3/ c bulk modulus bb=ud(6) call pzero(r,84) c integration points and weights on unit sphere r(3,1) = 1 r(6,1) = 0.707106781187 r(7,1) = -0.707106781187 r(8,1) = 0.707106781187 r(9,1) = -0.707106781187 r(10,1) = 0.836095596749 r(11,1) = -0.836095596749 r(12,1) = 0.836095596749 r(13,1) = -0.836095596749 r(14,1) = 0.387907304067 r(15,1) = -0.387907304067 r(16,1) = 0.387907304067 r(17,1) = -0.387907304067 r(18,1) = 0.387907304067 r(19,1) = -0.387907304067 r(20,1) = 0.387907304067 r(21,1) = -0.387907304067 r(2,2) = 1 r(4,2) = 0.707106781187 r(5,2) = -0.707106781187 r(8,2) = 0.707106781187 r(9,2) = 0.707106781187 r(10,2) = 0.387907304067 r(11,2) = 0.387907304067 r(12,2) = -0.387907304067 r(13,2) = -0.387907304067 r(14,2) = 0.836095596749 r(15,2) = 0.836095596749 r(16,2) = -0.836095596749 r(17,2) = -0.836095596749 r(18,2) = 0.387907304067 r(19,2) = 0.387907304067 r(20,2) = -0.387907304067 r(21,2) = -0.387907304067 r(1,3) = 1 r(4,3) = 0.707106781187 r(5,3) = 0.707106781187 r(6,3) = 0.707106781187 r(7,3) = 0.707106781187 r(10,3) = 0.387907304067 r(11,3) = 0.387907304067 r(12,3) = 0.387907304067 r(13,3) = 0.387907304067 r(14,3) = 0.387907304067 r(15,3) = 0.387907304067 r(16,3) = 0.387907304067 r(17,3) = 0.387907304067 r(18,3) = 0.836095596749 r(19,3) = 0.836095596749 r(20,3) = 0.836095596749 r(21,3) = 0.836095596749 r(1,4) = 0.0265214244093 r(2,4) = 0.0265214244093 r(3,4) = 0.0265214244093 r(4,4) = 0.0199301476312 r(5,4) = 0.0199301476312 r(6,4) = 0.0199301476312 r(7,4) = 0.0199301476312 r(8,4) = 0.0199301476312 r(9,4) = 0.0199301476312 r(10,4) = 0.0250712367487 r(11,4) = 0.0250712367487 r(12,4) = 0.0250712367487 r(13,4) = 0.0250712367487 r(14,4) = 0.0250712367487 r(15,4) = 0.0250712367487 r(16,4) = 0.0250712367487 r(17,4) = 0.0250712367487 r(18,4) = 0.0250712367487 r(19,4) = 0.0250712367487 r(20,4) = 0.0250712367487 r(21,4) = 0.0250712367487 c non-affine stretch part call pzero(t,63) do j = 1,21 do i = 1,3 do k = 1,3 t(i,j) = t(i,j) + fb(i,k)*r(j,k) end do end do end do call pzero(lab,21) do j = 1,21 do i = 1,3 lab(j) = lab(j)+t(i,j)**(2) end do lab(j) = sqrt(lab(j)) end do la=0 do i = 1,21 la = la + (lab(i)**(ud(1)))*r(i,4)*2.d0 end do la = la**(1./ud(1)) call pzero(h,6) call pzero(hh,36) do i = 1,21 do j = 1,6 h(j)=h(j)+t(co(j,1),i)*t(co(j,2),i) 1 *2.d0*r(i,4)*lab(i)**(ud(1)-2) do k= j,6 hh(j,k)=hh(j,k)+t(co(j,1),i)*t(co(j,2),i) 1 *t(co(k,1),i)*t(co(k,2),i) 2 *2.d0*r(i,4)*(ud(1)-2)*lab(i)**(ud(1)-4) end do end do end do c micro-stresses & micro-moduli tauf=ud(3)*la*(3.d0*ud(2)-la**(2))/(ud(2)-la**(2)) cf=ud(3)*(la**(4)+3*ud(2)**(2))/(ud(2)-la**(2))**(2) c macro-stresses & macro moduli call pzero(cfb,36) const=cf*la**(2.d0-2.d0*ud(1))-(ud(1)-1)*tauf*la**(1-2.d0*ud(1)) do i = 1,6 taufb(i)=h(i)*tauf*la**(1-ud(1)) do j = i,6 cfb(i,j)=hh(i,j)*tauf*la**(1-ud(1)) 1 +h(i)*h(j)*const end do end do c non-affine tube part det = fb(1,1)*(fb(2,2)*fb(3,3)-fb(2,3)*fb(3,2))- 1 fb(1,2)*(fb(2,1)*fb(3,3)-fb(2,3)*fb(3,1))+ 2 fb(1,3)*(fb(2,1)*fb(3,2)-fb(2,2)*fb(3,1)) C inverse transpose of fb fbit(1,1) = (fb(2,2)*fb(3,3)-fb(2,3)*fb(3,2))/det fbit(1,2) = (fb(2,3)*fb(3,1)-fb(2,1)*fb(3,3))/det fbit(1,3) = (fb(2,1)*fb(3,2)-fb(2,2)*fb(3,1))/det fbit(2,1) = (fb(1,3)*fb(3,2)-fb(1,2)*fb(3,3))/det fbit(2,2) = (fb(1,1)*fb(3,3)-fb(1,3)*fb(3,1))/det fbit(2,3) = (fb(1,2)*fb(3,1)-fb(1,1)*fb(3,2))/det fbit(3,1) = (fb(1,2)*fb(2,3)-fb(1,3)*fb(2,2))/det fbit(3,2) = (fb(1,3)*fb(2,1)-fb(1,1)*fb(2,3))/det fbit(3,3) = (fb(1,1)*fb(2,2)-fb(1,2)*fb(2,1))/det call pzero(n,63) do j = 1,21 do i = 1,3 do k = 1,3 n(i,j) = n(i,j)+fbit(i,k)*r(j,k) end do end do end do call pzero(vb,21) do i = 1,21 do j = 1,3 vb(i) = vb(i)+n(j,i)**(2) end do vb(i) = sqrt(vb(i)) end do call pzero(kk,6) call pzero(kkk,36) call pzero(G,36) call pzero(ij,9) do i = 1,3 ij(i,i) = 1.d0 end do do i = 1,21 do a = 1,6 kk(a)=kk(a)+n(co(a,1),i) 1 *n(co(a,2),i) 2 *2.d0*ud(4)*r(i,4)*vb(i)**(ud(4)-2) do b = a,6 kkk(a,b)=kkk(a,b) 1 +n(co(a,1),i)*n(co(a,2),i) 2 *n(co(b,1),i)*n(co(b,2),i) 3 *(ud(4)-2)*ud(4) 4 *2.d0*r(i,4)*vb(i)**(ud(4)-4) tmp1=1.d0/2.d0*(ij(co(a,1),co(b,1))*n(co(a,2),i)*n(co(b,2),i) 1 + ij(co(a,2),co(b,1))*n(co(a,1),i)*n(co(b,2),i)) tmp2=1.d0/2.d0*(ij(co(a,2),co(b,2))*n(co(a,1),i)*n(co(b,1),i) 1 + ij(co(a,1),co(b,2))*n(co(a,2),i)*n(co(b,1),i)) G(a,b)=G(a,b)+ud(4)*(tmp1+tmp2) 1 *4.d0*r(i,4)*vb(i)**(ud(4)-2) end do end do end do c macro-stresses and macro-moduli tmp7 = ud(3)*ud(2)*ud(5) do i = 1,6 taucb(i) = -tmp7*kk(i) do j = i,6 ccb(i,j) = tmp7*(kkk(i,j)+G(i,j)) end do end do c superimposed stress response do i = 1,6 taub(i) = taufb(i)+taucb(i) do j = i,6 cb(i,j) = cfb(i,j) + ccb(i,j) end do end do c pressure and bulk modulus p = bb*(theta*theta-1)/2.d0 bulk = bb*(theta*theta+1)/2.d0 c cauchy stress (sig) tmp3=(taub(1)+taub(2)+taub(3))/3.d0 sig(1)=(p-tmp3+taub(1))/theta sig(2)=(p-tmp3+taub(2))/theta sig(3)=(p-tmp3+taub(3))/theta sig(4)=( taub(4))/theta sig(5)=( taub(5))/theta sig(6)=( taub(6))/theta c calculate final stiffness tmp4=cb(1,1)+2*cb(1,2)+2*cb(1,3) 1 +cb(2,2)+2*cb(2,3) 2 +cb(3,3) tmp5(1)=cb(1,1)+cb(1,2)+cb(1,3) tmp5(2)=cb(1,2)+cb(2,2)+cb(2,3) tmp5(3)=cb(1,3)+cb(2,3)+cb(3,3) tmp5(4)=cb(1,4)+cb(2,4)+cb(3,4) tmp5(5)=cb(1,5)+cb(2,5)+cb(3,5) tmp5(6)=cb(1,6)+cb(2,6)+cb(3,6) c projection of taub tmp6(1)=taub(1)-tmp3 tmp6(2)=taub(2)-tmp3 tmp6(3)=taub(3)-tmp3 tmp6(4)=taub(4) tmp6(5)=taub(5) tmp6(6)=taub(6) do i=1,6 do j=i,6 dd(i,j)=((p+bulk-2.d0*tmp3/3.d0+(1.d0/9.d0)*tmp4) 1 *ij(co(i,1),co(i,2))*ij(co(j,1),co(j,2)) 2 +(ij(co(i,1),co(j,1))*ij(co(i,2),co(j,2)) 3 +ij(co(i,1),co(j,2))*ij(co(i,2),co(j,1)))*(tmp3-p) 4 +cb(i,j)-(1.d0/3.d0)*tmp5(i)*ij(co(j,1),co(j,2)) 5 -(1.d0/3.d0)*tmp5(j)*ij(co(i,1),co(i,2)) 6 -(2.d0/3.d0)*(ij(co(i,1),co(i,2))*tmp6(j) 7 +ij(co(j,1),co(j,2))*tmp6(i)))/theta end do end do do i=1,5 do j=i+1,6 dd(j,i)=dd(i,j) end do end do end