橡胶Mooney-Rivlin材料模型编程

woshimaidi

最近在往自己开发的软件上添加橡胶Mooney-Rivlin材料模型，但是结果不知道为什么老是不对。我用fortran95写的代码。把代码贴出来，有熟悉的同志点拨一下！我个人觉得在计算得到的变形梯度应该是没问题的（deformation gradients in undeformed configuration），应该是计算柯西-格林张量，第二类PK应力，柯西应力的时候少乘了一些量。谢谢了！

subroutine constitution_MR_material(eid,mid,ihq,pkxj,timestep,sig,vol)    !hyperelastic mooney-rivlin material  
!w = a*(i-3) + b*(ii-3)
!            + c*[iii**(-2)-1] + d*(iii-1)**2
!
!      where:   c = a/2 + b
!               d = [a*(5*pr-2) + b*(11*pr-5)] / (2-4*pr)
!               a , b , pr  are input variables
!               i , ii , iii  are invariants of right c-g tensor   
  implicit none
  integer,intent(in)::eid,mid,ihq
  real(kind=RK),intent(in)::pkxj(3,8),timestep
  real(kind=RK),intent(out)::sig(6)
  real(kind=RK),intent(in),optional::vol
  real(kind=RK)::tab(2),tc,dd,scb
  real(kind=RK)::pr
  integer::nodes(8)
  integer::size,i
  real(kind=RK)::xyz(3,8)
  real(kind=RK)::xyz17(3),xyz28(3),xyz35(3),xyz46(3)
  real(kind=RK)::eps(9),eg(6),trace,egi(6),cdfac,sp(6),q(9)
  real(kind=RK)::cdet,df,rx,cdeti
  character(len=MNL)::msg

    call obtain_MR_material( mid, tab)                                                                       !获得a,b的值
    pr=get_poisson_ratio( mid )                                                                               !获得泊松比   
    tc=0.5*tab(1) + tab(2)                                                                                      !c 的值
    dd=(tab(1)*(5*pr-2.) + tab(2)*(11.*pr-5.))/(2.-4.*pr)                                      !d 的值
    scb=2.*(tab(1)+tab(2))*(1.+pr)/(3.-6.*pr)                                                        !声速

  !compute deformation gradients in undeformed configuration
    size=get_fem_nid(eid,nodes)
    do i=1,8
       call obtain_node_kinematic(nodes(i),pos=xyz(:,i))                                              !获得节点坐标。
    end do
    if (ihq.ne.3.and.ihq.lt.5) then
       xyz17=xyz(:,1)-xyz(:,7)
       xyz28=xyz(:,2)-xyz(:,8)
       xyz35=xyz(:,3)-xyz(:,5)
       xyz46=xyz(:,4)-xyz(:,6)
       eps(1)=pkxj(1,1)*xyz17(1)+pkxj(1,2)*xyz28(1)+pkxj(1,3)*xyz35(1)+pkxj(1,4)*xyz46(1)
       eps(2)=pkxj(2,1)*xyz17(1)+pkxj(2,2)*xyz28(1)+pkxj(2,3)*xyz35(1)+pkxj(2,4)*xyz46(1)
       eps(3)=pkxj(3,1)*xyz17(1)+pkxj(3,2)*xyz28(1)+pkxj(3,3)*xyz35(1)+pkxj(3,4)*xyz46(1)
       eps(4)=pkxj(1,1)*xyz17(2)+pkxj(1,2)*xyz28(2)+pkxj(1,3)*xyz35(2)+pkxj(1,4)*xyz46(2)
       eps(5)=pkxj(2,1)*xyz17(2)+pkxj(2,2)*xyz28(2)+pkxj(2,3)*xyz35(2)+pkxj(2,4)*xyz46(2)
       eps(6)=pkxj(3,1)*xyz17(2)+pkxj(3,2)*xyz28(2)+pkxj(3,3)*xyz35(2)+pkxj(3,4)*xyz46(2)
       eps(7)=pkxj(1,1)*xyz17(3)+pkxj(1,2)*xyz28(3)+pkxj(1,3)*xyz35(3)+pkxj(1,4)*xyz46(3)
       eps(8)=pkxj(2,1)*xyz17(3)+pkxj(2,2)*xyz28(3)+pkxj(2,3)*xyz35(3)+pkxj(2,4)*xyz46(3)
       eps(9)=pkxj(3,1)*xyz17(3)+pkxj(3,2)*xyz28(3)+pkxj(3,3)*xyz35(3)+pkxj(3,4)*xyz46(3)
     else
       eps(1)=pkxj(1,1)*xyz(1,1)+pkxj(1,2)*xyz(1,2)+pkxj(1,3)*xyz(1,3)+pkxj(1,4)*xyz(1,4)+&
              &pkxj(1,5)*xyz(1,5)+pkxj(1,6)*xyz(1,6)+pkxj(1,7)*xyz(1,7)+pkxj(1,8)*xyz(1,8)
       eps(5)=pkxj(2,1)*xyz(2,1)+pkxj(2,2)*xyz(2,2)+pkxj(2,3)*xyz(2,3)+pkxj(2,4)*xyz(2,4)+&
              &pkxj(2,5)*xyz(2,5)+pkxj(2,6)*xyz(2,6)+pkxj(2,7)*xyz(2,7)+pkxj(2,8)*xyz(2,8)        
       eps(9)=pkxj(3,1)*xyz(3,1)+pkxj(3,2)*xyz(3,2)+pkxj(3,3)*xyz(3,3)+pkxj(3,4)*xyz(3,4)+&
              &pkxj(3,5)*xyz(3,5)+pkxj(3,6)*xyz(3,6)+pkxj(3,7)*xyz(3,7)+pkxj(3,8)*xyz(3,8)  
       eps(2)=pkxj(2,1)*xyz(1,1)+pkxj(2,2)*xyz(1,2)+pkxj(2,3)*xyz(1,3)+pkxj(2,4)*xyz(1,4)+&
              &pkxj(2,5)*xyz(1,5)+pkxj(2,6)*xyz(1,6)+pkxj(2,7)*xyz(1,7)+pkxj(2,8)*xyz(1,8)  
       eps(3)=pkxj(3,1)*xyz(1,1)+pkxj(3,2)*xyz(1,2)+pkxj(3,3)*xyz(1,3)+pkxj(3,4)*xyz(1,4)+&
              &pkxj(3,5)*xyz(1,5)+pkxj(3,6)*xyz(1,6)+pkxj(3,7)*xyz(1,7)+pkxj(3,8)*xyz(1,8)  
       eps(6)=pkxj(3,1)*xyz(2,1)+pkxj(3,2)*xyz(2,2)+pkxj(3,3)*xyz(2,3)+pkxj(3,4)*xyz(2,4)+&
              &pkxj(3,5)*xyz(2,5)+pkxj(3,6)*xyz(2,6)+pkxj(3,7)*xyz(2,7)+pkxj(3,8)*xyz(2,8)  
       eps(4)=pkxj(1,1)*xyz(2,1)+pkxj(1,2)*xyz(2,2)+pkxj(1,3)*xyz(2,3)+pkxj(1,4)*xyz(2,4)+&
              &pkxj(1,5)*xyz(2,5)+pkxj(1,6)*xyz(2,6)+pkxj(1,7)*xyz(2,7)+pkxj(1,8)*xyz(2,8)  
       eps(7)=pkxj(1,1)*xyz(3,1)+pkxj(1,2)*xyz(3,2)+pkxj(1,3)*xyz(3,3)+pkxj(1,4)*xyz(3,4)+&
              &pkxj(1,5)*xyz(3,5)+pkxj(1,6)*xyz(3,6)+pkxj(1,7)*xyz(3,7)+pkxj(1,8)*xyz(3,8)  
       eps(8)=pkxj(2,1)*xyz(3,1)+pkxj(2,2)*xyz(3,2)+pkxj(2,3)*xyz(3,3)+pkxj(2,4)*xyz(3,4)+&
              &pkxj(2,5)*xyz(3,5)+pkxj(2,6)*xyz(3,6)+pkxj(2,7)*xyz(3,7)+pkxj(2,8)*xyz(3,8)  
       if(.not.present(vol)) goto 111
       eps=eps/vol
     end if

    !right cauchy-green tensor
     eg(1)=eps(1)*eps(1)+eps(4)*eps(4)+eps(7)*eps(7)
     eg(2)=eps(2)*eps(2)+eps(5)*eps(5)+eps(8)*eps(8)
     eg(3)=eps(3)*eps(3)+eps(6)*eps(6)+eps(9)*eps(9)
     eg(4)=eps(1)*eps(2)+eps(4)*eps(5)+eps(7)*eps(8)
     eg(5)=eps(2)*eps(3)+eps(5)*eps(6)+eps(8)*eps(9)
     eg(6)=eps(1)*eps(3)+eps(4)*eps(6)+eps(7)*eps(9)

!     determinant, trace, and inverse
      df=eps(1)*(eps(5)*eps(9)-eps(6)*eps(8))+eps(4)*(eps(3)*eps(8)-eps(2)*eps(9))+eps(7)*(eps(2)*eps(6)-eps(3)*eps(5))
      cdet=df*df
      rx=1./df
      cdeti=rx*rx
      trace=eg(1)+eg(2)+eg(3)
      egi(1)=cdeti*(eg(2)*eg(3)-eg(5)*eg(5))
      egi(2)=cdeti*(eg(1)*eg(3)-eg(6)*eg(6))
      egi(3)=cdeti*(eg(1)*eg(2)-eg(4)*eg(4))
      egi(4)=cdeti*(eg(5)*eg(6)-eg(4)*eg(3))
      egi(5)=cdeti*(eg(4)*eg(6)-eg(1)*eg(5))
      egi(6)=cdeti*(eg(4)*eg(5)-eg(2)*eg(6))

!     second p-k stress
      cdfac=2.*dd*cdet*(cdet-1.)-2.*tc*cdeti*cdeti
      sp(1)=tab(1)+tab(2)*(trace-eg(1))+cdfac*egi(1)
      sp(2)=tab(1)+tab(2)*(trace-eg(2))+cdfac*egi(2)
      sp(3)=tab(1)+tab(2)*(trace-eg(3))+cdfac*egi(3)
      sp(4)=-eg(4)*tab(2)+cdfac*egi(4)
      sp(5)=-eg(5)*tab(2)+cdfac*egi(5)
      sp(6)=-eg(6)*tab(2)+cdfac*egi(6)

!     cauchy stress
      q(1)=eps(1)*sp(1)+eps(2)*sp(4)+eps(3)*sp(6)
      q(2)=eps(2)*sp(2)+eps(1)*sp(4)+eps(3)*sp(5)
      q(3)=eps(3)*sp(3)+eps(1)*sp(6)+eps(2)*sp(5)
      q(4)=eps(4)*sp(1)+eps(5)*sp(4)+eps(6)*sp(6)
      q(5)=eps(5)*sp(2)+eps(4)*sp(4)+eps(6)*sp(5)
      q(6)=eps(6)*sp(3)+eps(4)*sp(6)+eps(5)*sp(5)
      q(7)=eps(7)*sp(1)+eps(8)*sp(4)+eps(9)*sp(6)
      q(8)=eps(8)*sp(2)+eps(7)*sp(4)+eps(9)*sp(5)
      q(9)=eps(9)*sp(3)+eps(7)*sp(6)+eps(8)*sp(5)

!    calculate stress
     sig(1)=rx*(eps(1)*q(1)+eps(2)*q(2)+eps(3)*q(3))            
      sig(2)=rx*(eps(4)*q(4)+eps(5)*q(5)+eps(6)*q(6))
      sig(3)=rx*(eps(7)*q(7)+eps(8)*q(8)+eps(9)*q(9))
      sig(4)=rx*(eps(4)*q(1)+eps(5)*q(2)+eps(6)*q(3))
      sig(5)=rx*(eps(7)*q(4)+eps(8)*q(5)+eps(9)*q(6))
      sig(6)=rx*(eps(7)*q(1)+eps(8)*q(2)+eps(9)*q(3))
      return
end subroutine

cbgtc

楼主自己编感觉比较强大啊哈哈。不过目前很多有限元软件里面已经有Mooney-Rivlin本构了。

有个问题请教一下楼主，我看到有文献说没有充足实验数据的时候可以假设C01=0.25C10，请问一下楼主这个合适么？有没有出处？