请教怎么求 四点薄膜或薄壳单元的变形梯度

chenqian1982

怎么求四点薄膜或薄壳单元的变形梯度？以及在厚度方向的积分点得值是多少？（在整体坐标的）

chenqian1982

终于搞定了

molen

那就给大伙介绍一下吧，谢谢

chenqian1982

好啊，发一个子程序U 是位移，ETA KSI是高师点,N型函数,DNKSI ,DNETA 型函数求导，COORDS坐标
因为我的子程序是给ABAQUS的，所以含有KBLOCK,和NBLOCK，其实可以去调，NNODE点个数,NDOFEL自由度个数
C**********************************************************************
      SUBROUTINE CalculerF(eta,ksi,g,N,dNksi,dNeta,coords,u,nblock,
        1                                         nnode,ncrd,ndofel,kblock,epais,F,G1,G2)
c      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*4        dNksi(4),dNeta(4),N(4),u(nblock,ndofel),m(3),s(3),v(3),
        1                coords(nblock,nnode,ncrd),J(3,3),dNk(3,12),dNe(3,12),
        2                dNkv(3,1),dNev(3,1),node(12,1),dNkc(3,1),dNec(3,1),
     3                ev(12,1),EN(3,12),NV(3,1),INVJ(3,3),nodenew(12,1),
     4        dNkvn(3,1),dNevn(3,1),dF(3,3),F(3,3),G1(3,1),G2(3,1),
        5                dNksieta(3),dveta(3),dvksi(3),du(12,1),mn(3),sn(3),vn(3),
        7                dvetan(3),dvksin(3),dNksietan(3),df1(3,3)

C**********************************************************************       

        do k1=1,12
        node(k1,1)=0
        nodenew(k1,1)=0
        enddo

        do k1=1,3
         do k2=1,12
         dNe(k1,k2)        =0
         dNk(k1,k2)        =0
         EN(k1,k2) =0
         enddo
        enddo
               
        do k=1,4
        do i=1,3
         node((k-1)*3+i,1)=coords(k,i)
         nodenew((k-1)*3+i,1)=coords(k,i)+u(（k-1)*3+i)
        enddo
        enddo

        do i=1,3
         do k1=1,3
        J(i,k1)=0
        df(i,k1)=0
        f(i,k1)=0
        enddo
        enddo

        m(1)=dNksi(1)*coords(kblock,1,1)+dNksi(2)*coords(kblock,2,1)+
        1         dNksi(3)*coords(kblock,3,1)+dNksi(4)*coords(kblock,4,1)
        m(2)=dNksi(1)*coords(kblock,1,2)+dNksi(2)*coords(kblock,2,2)+
        1         dNksi(3)*coords(kblock,3,2)+dNksi(4)*coords(kblock,4,2)
        m(3)=dNksi(1)*coords(kblock,1,3)+dNksi(2)*coords(kblock,2,3)+
        1         dNksi(3)*coords(kblock,3,3)+dNksi(4)*coords(kblock,4,3)

        s(1)=dNeta(1)*coords(kblock,1,1)+dNeta(2)*coords(kblock,2,1)+
        1         dNeta(3)*coords(kblock,3,1)+dNeta(4)*coords(kblock,4,1)
        s(2)=dNeta(1)*coords(kblock,1,2)+dNeta(2)*coords(kblock,2,2)+
        1         dNeta(3)*coords(kblock,3,2)+dNeta(4)*coords(kblock,4,2)
        s(3)=dNeta(1)*coords(kblock,1,3)+dNeta(2)*coords(kblock,2,3)+
        1         dNeta(3)*coords(kblock,3,3)+dNeta(4)*coords(kblock,4,3)

        dNksieta(1)=(coords(kblock,1,1)-coords(kblock,2,1)+
        1                        coords(kblock,3,1)-coords(kblock,4,1))/4

        dNksieta(2)=(coords(kblock,1,2)-coords(kblock,2,2)+
        1                        coords(kblock,3,2)-coords(kblock,4,2))/4

        dNksieta(3)=(coords(kblock,1,3)-coords(kblock,2,3)+
        1                        coords(kblock,3,3)-coords(kblock,4,3))/4

        dvksi(1)=m(2)*dNksieta(3)-m(3)*dNksieta(2)
        dvksi(2)=m(3)*dNksieta(1)-m(1)*dNksieta(3)
        dvksi(3)=m(1)*dNksieta(2)-m(2)*dNksieta(1)

        dveta(1)=s(3)*dNksieta(2)-s(2)*dNksieta(3)
        dveta(2)=s(1)*dNksieta(3)-s(3)*dNksieta(1)
        dveta(3)=s(2)*dNksieta(1)-s(1)*dNksieta(2)

        v(1)=m(2)*s(3)-m(3)*s(2)
        v(2)=m(3)*s(1)-m(1)*s(3)
        v(3)=m(1)*s(2)-m(2)*s(1)

        dbarv=sqrt(v(1)**2+v(2)**2+v(3)**2)

        v(1)=v(1)/dbarv
        v(2)=v(2)/dbarv
        v(3)=v(3)/dbarv

        do k1=1,3
        dNe(k1,k1)        =dNeta(1)
        dNe(k1,k1+3)=dNeta(2)
        dNe(k1,k1+6)=dNeta(3)
        dNe(k1,k1+9)=dNeta(4)

        dNk(k1,k1)        =dNksi(1)
        dNk(k1,k1+3)=dNksi(2)
        dNk(k1,k1+6)=dNksi(3)
        dNk(k1,k1+9)=dNksi(4)
        enddo       

        call BRMUL(dNk,node,3,12,1,dNkv)
        call BRMUL(dNe,node,3,12,1,dNev)


        J(1,1)=dNkv(1,1)+(dNksi(1)+dNksi(3)+dNksi(2)+dNksi(4))*v(1)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dvksi(1)*g

        J(1,2)=dNkv(2,1)+(dNksi(1)+dNksi(3)+dNksi(2)+dNksi(4))*v(2)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dvksi(2)*g

        J(1,3)=dNkv(3,1)+(dNksi(1)+dNksi(3)+dNksi(2)+dNksi(4))*v(3)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dvksi(3)*g

        J(2,1)=dNev(1,1)+(dNeta(1)+dNeta(3)+dNeta(2)+dNeta(4))*v(1)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dveta(1)*g

        J(2,2)=dNev(2,1)+(dNeta(1)+dNeta(3)+dNeta(2)+dNeta(4))*v(2)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dveta(2)*g

        J(2,3)=dNev(3,1)+(dNeta(1)+dNeta(3)+dNeta(2)+dNeta(4))*v(3)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dveta(3)*g

        J(3,1)=(N(1)+N(3)+N(2)+N(4))*v(1)

        J(3,2)=(N(1)+N(3)+N(2)+N(4))*v(2)

        J(3,3)=(N(1)+N(3)+N(2)+N(4))*v(3)

c        calcul dxn/di

        mn(1)=dNksi(1)*nodenew(1,1)+dNksi(2)*nodenew(4,1)+
        1         dNksi(3)*nodenew(7,1)+dNksi(4)*nodenew(10,1)
        mn(2)=dNksi(1)*nodenew(2,1)+dNksi(2)*nodenew(5,1)+
        1         dNksi(3)*nodenew(8,1)+dNksi(4)*nodenew(11,1)
        mn(3)=dNksi(1)*nodenew(3,1)+dNksi(2)*nodenew(6,1)+
        1         dNksi(3)*nodenew(9,1)+dNksi(4)*nodenew(12,1)

        sn(1)=dNeta(1)*nodenew(1,1)+dNeta(2)*nodenew(4,1)+
        1         dNeta(3)*nodenew(7,1)+dNeta(4)*nodenew(10,1)
        sn(2)=dNeta(1)*nodenew(2,1)+dNeta(2)*nodenew(5,1)+
        1         dNeta(3)*nodenew(8,1)+dNeta(4)*nodenew(11,1)
        sn(3)=dNeta(1)*nodenew(3,1)+dNeta(2)*nodenew(6,1)+
        1         dNeta(3)*nodenew(9,1)+dNeta(4)*nodenew(12,1)

        dNksietan(1)=(nodenew(1,1)-nodenew(4,1)+
        1                         nodenew(7,1)-nodenew(10,1))/4
        dNksietan(2)=(nodenew(2,1)-nodenew(5,1)+
        1                         nodenew(8,1)-nodenew(11,1))/4
        dNksietan(3)=(nodenew(3,1)-nodenew(6,1)+
        1                         nodenew(9,1)-nodenew(12,1))/4

        dvksin(1)=mn(2)*dNksietan(3)-mn(3)*dNksietan(2)
        dvksin(2)=mn(3)*dNksietan(1)-mn(1)*dNksietan(3)
        dvksin(3)=mn(1)*dNksietan(2)-mn(2)*dNksietan(1)

        dvetan(1)=sn(2)*dNksietan(3)-sn(3)*dNksietan(2)
        dvetan(2)=sn(3)*dNksietan(1)-sn(1)*dNksietan(3)
        dvetan(3)=sn(1)*dNksietan(2)-sn(2)*dNksietan(1)

        vn(1)=mn(2)*sn(3)-mn(3)*sn(2)
        vn(2)=mn(3)*sn(1)-mn(1)*sn(3)
        vn(3)=mn(1)*sn(2)-mn(2)*sn(1)



        dbarvn=sqrt(vn(1)**2+vn(2)**2+vn(3)**2)

        vn(1)=vn(1)/dbarv
        vn(2)=vn(2)/dbarv
        vn(3)=vn(3)/dbarv

        call BRMUL(dNk,nodenew,3,12,1,dNkvn)
        call BRMUL(dNe,nodenew,3,12,1,dNevn)

        dF(1,1)=dNkvn(1,1)+(dNksi(1)+dNksi(3)+dNksi(2)+dNksi(4))*vn(1)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dvksin(1)*g

        dF(1,2)=dNkvn(2,1)+(dNksi(1)+dNksi(3)+dNksi(2)+dNksi(4))*vn(2)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dvksin(2)*g

        dF(1,3)=dNkvn(3,1)+(dNksi(1)+dNksi(3)+dNksi(2)+dNksi(4))*vn(3)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dvksin(3)*g

        dF(2,1)=dNevn(1,1)+(dNeta(1)+dNeta(3)+dNeta(2)+dNeta(4))*vn(1)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dvetan(1)*g

        dF(2,2)=dNevn(2,1)+(dNeta(1)+dNeta(3)+dNeta(2)+dNeta(4))*vn(2)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dvetan(2)*g

        dF(2,3)=dNevn(3,1)+(dNeta(1)+dNeta(3)+dNeta(2)+dNeta(4))*vn(3)*g+
        1           (N(1)+N(3)+N(2)+N(4))*dvetan(3)*g

        dF(3,1)=(N(1)+N(3)+N(2)+N(4))*vn(1)

        dF(3,2)=(N(1)+N(3)+N(2)+N(4))*vn(2)

        dF(3,3)=(N(1)+N(3)+N(2)+N(4))*vn(3)

        call KMINV(J,INVJ)  求逆
      CALL BRMUL(INVJ,dF,3,3,3,dF1) 矩阵乘法

        do k1=1,3
         do k2=1,3
          F(k1,k2)=dF1(k2,k1)
         enddo
        enddo
               
        RETURN
        END

chenqian1982

这个求出来的变形梯度是关于整体坐标下的。应为结构单元一般用的是随体坐标，所以要求随体坐标的变形梯度还要做一个转动。那个应该比较容易了吧。

felixyu623

chenqian1982 发表于 2008-2-7 17:19
这个求出来的变形梯度是关于整体坐标下的。应为结构单元一般用的是随体坐标，所以要求随体坐标的变形梯度还 ...

你好，请问你所发的程序是用在VUEL里面的吗？