不知道怎样输入这些参数

wj-6608 · 发表于 2020-10-25 15:48:35

本帖最后由 wj-6608 于 2020-10-25 17:18 编辑

[Fortran] 纯文本查看 复制代码

1层数 3
2层序数    厚度m   E1(Pa)   E2(Pa)   v12    v21  G12(Pa)    角度（°）
31        1.00E-03     9.60E+10        0.10    0.40   1.00E+10   0.00
42        1.00E-03     9.60E+10        0.10    0.40   1.00E+10   90.00
53        1.00E-03     9.60E+10        0.10    0.40   1.00E+10   0.00

[Fortran] 纯文本查看 复制代码

001PROGRAM COMPOSITE
002IMPLICIT NONE
003REAL(8)::A(3,3),B(3,3),D(3,3),MC(5),TEMP,ROT(3,3)
004!A拉伸刚度；B耦合刚度；D弯曲刚度；
005!MC读入材料常数；ROT旋转矩阵
006REAL(8)::TOTAL_TH,HALF_TH  !总厚度；半厚度
007REAL(8),ALLOCATABLE::Q(:,:,:),AL(:),T(:),Z(:),Z1(:),Z2(:),Z3(:)
008!Q每层板相应刚度；AL转角；T每层板的厚度；Z坐标量
009INTEGER(4)::N,I,J,K,SEQ,L
010!____IJK循环变量；N板的层数；SEQ序数
011CHARACTER(50)::CHR(8),TEMPC,filename1,filename2
012!CHR、TEMPC:character variables
013WRITE(*,*)"Please insert the INP file name(a.txt for example):"
014READ(*,*)filename1
015OPEN(8,file=filename1)!Open data file
016!Read in data
017READ(8,*)TEMPC,N
018ALLOCATE(Q(3,3,N),AL(N),T(N),Z(N+1),Z1(N),Z2(N),Z3(N))
019READ(8,*)CHR(1:8)
020DO I=1,N
021READ(8,*)SEQ,T(I),MC(1:5),AL(I)
022Q(:,:,I)=0!Calculate stiffness of each layer for the principal axis
023TEMP=1./(1-MC(3)*MC(4))
024Q(1,1,I)=MC(1)*TEMP
025Q(2,2,I)=MC(2)*TEMP
026Q(3,3,I)=MC(5)
027Q(1,2,I)=MC(4)*MC(2)*TEMP
028Q(2,1,I)=Q(1,2,I)
029AL(I)=AL(I)*3.1415926535898/180
030ROT(1,1)=(cos(AL(I)))**2!Work out Rot Matrix
031ROT(2,2)=ROT(1,1)
032ROT(3,3)=cos(2*AL(I))
033ROT(2,1)=1-ROT(1,1)
034ROT(1,2)=ROT(2,1)
035ROT(3,1)=0.5*sin(2*AL(I))
036ROT(3,2)=-ROT(3,1)
037ROT(1,3)=-2*ROT(3,1)
038ROT(2,3)=-2*ROT(3,2)
039Q(:,:,I)=MATMUL(MATMUL(ROT,Q(:,:,I)),TRANSPOSE(ROT))
040ENDDO
041TOTAL_TH=sum(T)
042HALF_TH=TOTAL_TH/2
043Z(1)=-HALF_TH
044!Work out Z
045DO I=1,N
046Z(I+1)=Z(I)+T(I)
047ENDDO
048!calculate tensor A、B and D
049DO K=1,N
050Z1(K)=(Z(K+1)-Z(K))
051Z2(K)=(Z(K+1)-Z(K))*(Z(K+1)+Z(K))/2
052Z3(K)=(Z(K+1)**3-Z(K)**3)/3
053ENDDO
054A=0;B=0;D=0
055WRITE(*,*)"Please insert the OUP file name(b.txt for example):"
056READ(*,*)filename2
057OPEN(9,file=filename2)
058!Write in stiffness tensor for each single ply
059DO K=1,N
060WRITE(9,100)K
061100 FORMAT("The stiffness of number",1X,I2,2X,"ply is:")
062DO I=1,3
063WRITE(9,200)Q(I,:,K)
064200 FORMAT(ES12.4,6X,ES12.4,6X,ES12.4)
065ENDDO
066WRITE(9,"(/)")
067A=A+Q(:,:,K)*Z1(K)
068B=B+Q(:,:,K)*Z2(K)
069D=D+Q(:,:,K)*Z3(K)
070ENDDO
071!Output the actual stiffness of the laminate
072WRITE(9,"(/)");WRITE(9,"(/)")
073WRITE(9,*)"The ACTUAL stiffness tensor of the laminate:"
074WRITE(9,"(/)")
075WRITE(9,*)"The extension stiffness A equals:"
076DO I=1,3
077WRITE(9,200)A(I,1:3)
078ENDDO
079WRITE(9,"(/)")
080WRITE(9,*)"The coupling stiffness B equals:"
081DO I=1,3
082WRITE(9,200)B(I,1:3)
083ENDDO
084WRITE(9,"(/)")
085WRITE(9,*)"The bending stiffness D equals:"
086DO I=1,3
087WRITE(9,200)D(I,1:3)
088ENDDO
089!Normalised tensor output
090WRITE(9,"(/)");WRITE(9,"(/)")
091WRITE(9,*)"The NORMALISED stiffness tensor of the laminate:"
092WRITE(9,"(/)")
093WRITE(9,*)"The NORMALISED extension stiffness A* equals:"
094DO I=1,3
095WRITE(9,200)A(I,1:3)/TOTAL_TH
096ENDDO
097WRITE(9,"(/)")
098WRITE(9,*)"The NORMALISED coupling tensor C equals:"
099DO I=1,3
100WRITE(9,200)A(I,1:3)/TOTAL_TH-12*D(I,1:3)/TOTAL_TH**3
101ENDDO
102WRITE(9,"(/)")
103WRITE(9,*)"The NORMALISED bending stiffness D* equals:"
104DO I=1,3
105WRITE(9,200)12*D(I,1:3)/TOTAL_TH**3
106ENDDO
107WRITE(*,*)"OUTPUT successfully,please press any key to end program!"
108READ(*,*)
109END PROGRAM COMPOSITE