关于三维梁有限元坐标转换的一点问题

hao1982

最近在编三维梁有限元程序时有点问题,想请知道的指点下.谢谢

在二维梁有限元中,假设X,Y为整体坐标的两个方向,第I个单元的局部坐标X',Y'可以这样确定,设I单元的轴线方向的可以定为X',另一个是Y',这样问题比较简单.
但是在三维中,假设X,Y,Z为整体坐标的三个方向,X',Y',Z'为第I个单元的局部坐标方向,可以将X'仍然定为是第I个单元的轴线方向,但此时Y',Z'怎么确定,如果没确定好,输入抗弯惯性矩时哪个是对Y'轴的抗弯惯性矩,哪个是对Z'轴的抗弯惯性矩,另外坐标转换矩阵中也需要计算方向余弦cos(Y,Y')等,我想应该有一定的规矩方法吧,请大家指教.

欧阳中华

.
    这是编程序必然要遇到的问题。

    因为梁的节点都在一条线上，所以从二维组装到三维是不唯一，简单的办法就是在设一个辅助节点，通常是与梁组成一个唯一的面，那么就可以用这个面与空间三维的关系形成坐标转换矩阵。

    程序编制时，节点编号输入除I和J以外，再加个K，程序中要将K的坐标传进去。

   计算算例时，根据算例的节点情况，选择每个梁单元的合适K点，当然两单元的剖面参数也要与参考点定义的相对应。如果没有合适的参考点可以利用（如悬臂梁），则在结构网格剖分时，专门设个点用来做参考点，但要注意这个点的所有约束都要约束住，不然就会奇异。输出结果时不需要这个点的结果。

   应该说清楚了吧，中午请吃饭喽.. ..

hao1982

:@) :@)
但是感觉还需要好好看看书,一般的有限元书籍鲜有提到此种方法的,我想编程中才会碰到这个问题.
不知道中华兄有没有相关资料推荐或共享下
能把书名告诉我吗,或者有资料的话发到我邮箱的话更是感谢,我的邮箱是 hao19820125@yahoo.com.cn

[ 本帖最后由 hao1982 于 2007-3-10 13:16 编辑 ]

欧阳中华

.
   给段梁的程序你琢磨一下就明白了。

C1998.05.20--------------------------------------------------2001.06.06
C
C       This subroutine is used to carry out beam element matrices.
C
C  Input
C     ID[6,NUMNP] => Boundinary condition
C      X[3,NUMNP] => Coordinates of the nodes
C
C  Workstore
C        MHT[NEQ] => Active column heights
C  SIZE[6,NUMSIZ] => Section properties     Ao  Ay  Az  Jx  Iy  Iz
C     LM[12,NUME] => Storing the assemblage degrees of freedom of ele.
C     XYZ[9,NUME] => XYZ of two nodes of the element
C      MTYPE[NUME] => Size property of the element
C
C----------------------------------------------------------------------
      SUBROUTINE BEAM_INPUT(ID,X,MHT,NPAR,SIZE,LM,XYZ,MTYP,IN,IO,IC,IP)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ID(6,1),X(3,1),MHT(1),NPAR(1)
      DIMENSION SIZE(6,1),LM(12,1),XYZ(9,1),MTYP(1)
      ND = 2*6
C
C        输入当前单元组的几何与材料信息
C
      WRITE(IO,2000) NPAR(1),NPAR(2),NPAR(3)
      WRITE(IO,2030)
      DO 100 I = 1,NPAR(3)
      READ(IN, '(I5,6E10.3)') N,(SIZE(L,N),L=1,6)
      WRITE(IO,'(I6,6E11.3)') N,(SIZE(L,N),L=1,6)
      IF(N.NE.I) STOP
  100 CONTINUE
      READ(IN,*)
C
C    输入当前单元组各单元信息，并确定单元各结点自由度对应的方程编号
C
      WRITE(IO,2040)
        DO 110 N_ELEMENT=1,NPAR(2)
      READ( IN,'(6I5)') N,I,J,K,MTYP(N),MTYP(N)
      WRITE(IP,'(10X,2I5)') I,J
      WRITE(IO,'(5X,5(5X,I5))') N,I,J,K,MTYP(N)
        IF(N.NE.N_ELEMENT) STOP
      XYZ(1,N) = X(1,I)
      XYZ(2,N) = X(2,I)
      XYZ(3,N) = X(3,I)
      XYZ(4,N) = X(1,J)
      XYZ(5,N) = X(2,J)
      XYZ(6,N) = X(3,J)
      XYZ(7,N) = X(1,K)
      XYZ(8,N) = X(2,K)
      XYZ(9,N) = X(3,K)
      LM( 1,N) = ID(1,I)
      LM( 2,N) = ID(2,I)
      LM( 3,N) = ID(3,I)
      LM( 4,N) = ID(4,I)
      LM( 5,N) = ID(5,I)
      LM( 6,N) = ID(6,I)
      LM( 7,N) = ID(1,J)
      LM( 8,N) = ID(2,J)
      LM( 9,N) = ID(3,J)
      LM(10,N) = ID(4,J)
      LM(11,N) = ID(5,J)
      LM(12,N) = ID(6,J)
      CALL COLHT(MHT,LM(1,N),ND)  !计算各单元对系统总体矩阵的列高的影响
C      IF(N.EQ.1) THEN
C         WRITE(IC,'(/5X,20HNUMBER OF ELEMENT : ,I4/)') N
C       WRITE(IC,'(5X,6HNi :  ,2I8/)')  I,J
C       WRITE(IC,'(5X,6HUi :  ,2I8)') LM(1,N),LM( 7,N)
C       WRITE(IC,'(5X,6HVi :  ,2I8)') LM(2,N),LM( 8,N)
C       WRITE(IC,'(5X,6HWi :  ,2I8)') LM(3,N),LM( 9,N)
C       WRITE(IC,'(5X,6H0x :  ,2I8)') LM(4,N),LM(10,N)
C       WRITE(IC,'(5X,6H0y :  ,2I8)') LM(5,N),LM(11,N)
C       WRITE(IC,'(5X,6H0z :  ,2I8)') LM(6,N),LM(12,N)
C        ENDIF
  110        CONTINUE
      RETURN
C
2000 FORMAT(
     &  5X,'Elements type                 = ',I5/
     &  5X,'Number of elements            = ',I5/
     &  5X,'Number of different sizes     = ',I5/)
2030 FORMAT(//3X,'Sets',5X,'Ao',9X,'Ay',9X,'Az',9X,'Jx',9X,'Iy',
     &  9X,'Iz')
2040 FORMAT(/12X,
     &  'Number    Node-i    Node-j    Node-k    Size')
        END
C1998.05.20-------------------------------------------------2001.06.06
C
C       This subroutine is used to carry out beam element matrices .
C
C  Input
C       MAXA[NEQ] =>  Addresses of diagonal elements
C         GK[NWK] =>  Global stiffness matrix
C         GM[NWM] =>  Global mass matrix
C
C  Workstore
C        MHT[NEQ] => Active column heights
C  SIZE[6,NUMSIZ] => Section properties     Ao  Ay  Az  Jx  Iy  Iz
C     LM[12,NUME] => Storing the assemblage degrees of freedom of ele.
C     XYZ[9,NUME] => XYZ of two nodes of the element
C     MTYPE[NUME] => Size property of the element
C
C---------------------------------------------------------------------
      SUBROUTINE BEAM_FORM(NDYN,MAXA,GK,GM,NPAR,
     &          SIZE,LM,XYZ,MTYPE,ES,EK,ESK,EM,EG,EMG)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION MAXA(1),GK(1),GM(1),NPAR(1)
        DIMENSION SIZE(6,1),LM(12,1),XYZ(9,1),MTYPE(1),ET(3,3)
      DIMENSION ES(12,1),EK(12,1),ESK(1)
      DIMENSION EM(12,1),EG(12,1),EMG(1)
      ND = 2*6
C
C        计 算 当 前 单 元 组 各 单 元 的 单 元 矩 阵      
C
        EN = 2.1E11          ! 结构材料的弹性模量
        PO = 0.3             ! 结构材料的泊松比
        GN = EN/(2.*(1.+PO)) ! 结构材料的剪切模量
        PN = 7800.0          ! 结构材料的质量密度
      DO 200 N=1,NPAR(2)
      AX = SIZE(1,MTYPE(N))
      AY = SIZE(2,MTYPE(N))
      AZ = SIZE(3,MTYPE(N))
      EIX= SIZE(4,MTYPE(N))
      EIY= SIZE(5,MTYPE(N))
      EIZ= SIZE(6,MTYPE(N))
      A1 = XYZ(4,N) - XYZ(1,N)
      A2 = XYZ(5,N) - XYZ(2,N)
      A3 = XYZ(6,N) - XYZ(3,N)
      XL = SQRT(A1*A1 + A2*A2 + A3*A3)
      IF(XL.EQ.0.) THEN
        WRITE(*,'(5X,15HERROR IN *BEAM*,I5,18H HAS ZERO LENGTH !//)') N
        STOP
      ENDIF
      FY = 0.
      FZ = 0.
      IF(AY.GT.0.) FY = 12.*EN*EIZ/(GN*AY*XL*XL)
      IF(AZ.GT.0.) FZ = 12.*EN*EIY/(GN*AZ*XL*XL)
      CALL CONSD(ES,144)
      CALL CONSD(EM,144)
C
      ES( 1, 1) = EN*AX/XL
      ES( 7, 1) =-ES( 1, 1)
      ES( 2, 2) = 12.*EN*EIZ/(XL*XL*XL*(1.+FY))
      ES( 6, 2) =  6.*EN*EIZ/(XL*XL*(1.+FY))
      ES( 8, 2) =-ES( 2, 2)
      ES(12, 2) = ES( 6, 2)
      ES( 3, 3) = 12.*EN*EIY/(XL*XL*XL*(1.+FZ))
      ES( 5, 3) = -6.*EN*EIY/(XL*XL*(1.+FZ))
      ES( 9, 3) =-ES( 3, 3)
      ES(11, 3) = ES( 5, 3)
      ES( 4, 4) = GN*EIX/XL
      ES(10, 4) =-ES( 4, 4)
      ES( 5, 5) = (4.+FZ)*EN*EIY/(XL*(1.+FZ))
      ES( 9, 5) = 6.*EN*EIY/(XL*XL*(1.+FZ))
      ES(11, 5) = (2.-FZ)*EN*EIY/(XL*(1.+FZ))
      ES( 6, 6) = (4.+FY)*EN*EIZ/(XL*(1.+FY))
      ES( 8, 6) =-ES( 6, 2)
      ES(12, 6) = (2.-FY)*EN*EIZ/(XL*(1.+FY))
      ES( 7, 7) = ES( 1, 1)
      ES( 8, 8) = ES( 2, 2)
      ES(12, 8) = ES( 8, 6)
      ES( 9, 9) = ES( 3, 3)
      ES(11, 9) = ES( 9, 5)
      ES(10,10) = ES( 4, 4)
      ES(11,11) = ES( 5, 5)
      ES(12,12) = ES( 6, 6)
C
      EM( 1, 1) = 1./3.
      EM( 7, 1) = 1./6.
      EM( 2, 2) = 13./35.
      EM( 6, 2) = 11.*XL/210.
      EM( 8, 2) = 9./70.
      EM(12, 2) =-13.*XL/420.
      EM( 3, 3) = 13./35.
      EM( 5, 3) =-11.*XL/210.
      EM( 9, 3) = 9./70.
      EM(11, 3) = 13.*XL/420.
      EM( 4, 4) = EIX/(3.*AX)
      EM(10, 4) = EIX/(6.*AX)
      EM( 5, 5) = XL*XL/105.
      EM( 9, 5) =-13.*XL/420.
      EM(11, 5) =-XL*XL/140.
      EM( 6, 6) = XL*XL/105.
      EM( 8, 6) = 13.*XL/420.
      EM(12, 6) =-XL*XL/140.
      EM( 7, 7) = 1./3.
      EM( 8, 8) = 13./35.
      EM(12, 8) =-11.*XL/210.
      EM( 9, 9) = 13./35.
      EM(11, 9) = 11.*XL/210.
      EM(10,10) = EM( 4, 4)
      EM(11,11) = XL*XL/105.
      EM(12,12) = XL*XL/105.
C
      DO 100 I=1,12
      DO 100 J=1,I
      EM(I,J) = PN*AX*XL*EM(I,J)
      EM(J,I) = EM(I,J)
  100 ES(J,I) = ES(I,J)
C
C    利用矩阵分块运算技术，对当前单元组各单元的单元矩阵进行坐标转换
C
      B1 = XYZ(7,N) - XYZ(1,N)
      B2 = XYZ(8,N) - XYZ(2,N)
      B3 = XYZ(9,N) - XYZ(3,N)
      C1 = A2*B3 - A3*B2
      C2 = A3*B1 - A1*B3
      C3 = A1*B2 - A2*B1
      D1 = C2*A3 - C3*A2
      D2 = C3*A1 - C1*A3
      D3 = C1*A2 - C2*A1
      DL = SQRT(D1*D1 + D2*D2 + D3*D3)
      IF(DL.EQ.0.) THEN
        WRITE(*,'(5X,15HERROR IN *BEAM* ,I5,17H U EQUALS ZERO !//)') N
          WRITE(*,'(5X,6F8.3//)') (XYZ(I,N),I=1,3),(XYZ(I,N),I=7,9)
        STOP
      ENDIF
      ET(1,1) = A1/XL
      ET(1,2) = A2/XL
      ET(1,3) = A3/XL
      ET(2,1) = D1/DL
      ET(2,2) = D2/DL
      ET(2,3) = D3/DL
      ET(3,1) = ET(1,2)*ET(2,3) - ET(1,3)*ET(2,2)
      ET(3,2) = ET(1,3)*ET(2,1) - ET(1,1)*ET(2,3)
      ET(3,3) = ET(1,1)*ET(2,2) - ET(1,2)*ET(2,1)
      CALL CONSD(EG,144)
      CALL CONSD(EK,144)
      DO 120 LA=1,10,3
      LB = LA + 2
      DO 120 MA=1,10,3
      MB = MA - 1
      DO 120 I=LA,LB
      DO 120 JM=1,3
      J = JM + MB
      XX = 0.
      YY = 0.
      DO 110 K=1,3
      YY = YY + EM(I,K+MB)*ET(K,JM)
  110 XX = XX + ES(I,K+MB)*ET(K,JM)
      EG(I,J) = YY
  120 EK(I,J) = XX
      CALL CONSD(EM,144)
      CALL CONSD(ES,144)
      DO 140 LA=1,10,3
      LB = LA - 1
      DO 140 MA=1,10,3
      MB = MA + 2
      DO 140 IL=1,3
      I = IL + LB
      DO 140 J=MA,MB
      XX = 0.
      YY = 0.
      DO 130 K=1,3
      YY = YY + ET(K,IL)*EG(K+LB,J)
  130 XX = XX + ET(K,IL)*EK(K+LB,J)
      EM(I,J) = YY
  140 ES(I,J) = XX
C
      K = 1
      DO 150 I=1,12
      DO 150 J=I,12
      ESK(K) = ES(I,J)
  150 K = K + 1
C
        AX = (EM( 1, 1)+EM( 1, 7))/EM( 1, 1)
        AY = (EM( 2, 2)+EM( 8, 8))/EM( 2, 2)
        AZ = (EM( 3, 3)+EM( 9, 9))/EM( 3, 3)
      EMG( 1) = EM( 1, 1)*AX
      EMG( 2) = EM( 2, 2)*AY
        EMG( 3) = EM( 3, 3)*AZ
        EMG( 4) = EM( 4, 4)
        EMG( 5) = EM( 5, 5)*AZ
        EMG( 6) = EM( 6, 6)*AY
        EMG( 7) = EM( 7, 7)*AX
        EMG( 8) = EM( 8, 8)*AY
        EMG( 9) = EM( 3, 3)*AZ
        EMG(10) = EM(10,10)
        EMG(11) = EM(11,11)*AZ
        EMG(12) = EM(12,12)*AY
C      WRITE(5,'(5X,12E12.4)')  (EMG(I),I=1,12)
C
C      将当前单元组各单元的一维存储单元矩阵组装进系统总体矩阵中      
C
        IF(NDYN.EQ.1) THEN
      CALL ADDBAN_K(MAXA,GK,ESK,LM(1,N),ND)
        ENDIF
        IF(NDYN.GE.2) THEN
      CALL ADDBAN_K(MAXA,GK,ESK,LM(1,N),ND)
      CALL ADDBAN_M(MAXA,GM,EMG,LM(1,N),ND)
        ENDIF
  200 CONTINUE
      RETURN
      END
C---------------------------------------------------------------------

hao1982

谢谢