关于三维梁有限元坐标转换的一点问题
最近在编三维梁有限元程序时有点问题,想请知道的指点下.谢谢在二维梁有限元中,假设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点，当然两单元的剖面参数也要与参考点定义的相对应。如果没有合适的参考点可以利用（如悬臂梁），则在结构网格剖分时，专门设个点用来做参考点，但要注意这个点的所有约束都要约束住，不然就会奇异。输出结果时不需要这个点的结果。
应该说清楚了吧，中午请吃饭喽.. ..
说的很详细了,谢谢中华兄
:@) :@)但是感觉还需要好好看看书,一般的有限元书籍鲜有提到此种方法的,我想编程中才会碰到这个问题.
不知道中华兄有没有相关资料推荐或共享下
能把书名告诉我吗,或者有资料的话发到我邮箱的话更是感谢,我的邮箱是 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
CInput
C ID => Boundinary condition
C X => Coordinates of the nodes
C
CWorkstore
C MHT => Active column heights
CSIZE => Section properties AoAyAzJxIyIz
C LM => Storing the assemblage degrees of freedom of ele.
C XYZ => XYZ of two nodes of the element
C MTYPE => 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
CInput
C MAXA =>Addresses of diagonal elements
C GK =>Global stiffness matrix
C GM =>Global mass matrix
C
CWorkstore
C MHT => Active column heights
CSIZE => Section properties AoAyAzJxIyIz
C LM => Storing the assemblage degrees of freedom of ele.
C XYZ => XYZ of two nodes of the element
C MTYPE => 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---------------------------------------------------------------------
THANK YOU
谢谢
页:
[1]