请教怎么求 四点薄膜或薄壳单元的变形梯度
怎么求四点薄膜或薄壳单元的变形梯度？以及在厚度方向的积分点得值是多少？（在整体坐标的） 终于搞定了 那就给大伙介绍一下吧，谢谢 好啊，发一个子程序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 发表于 2008-2-7 17:19
这个求出来的变形梯度是关于整体坐标下的。应为结构单元一般用的是随体坐标，所以要求随体坐标的变形梯度还 ...
你好，请问你所发的程序是用在VUEL里面的吗？
页:
[1]