!*************************************************************************************** !* This program calculates the matrix linking the stresses to deformamations in a * !* laminated material made of n unidirectional composite layers. * !* * !* The 6x6 stiffness matrix to define has the following coefficients: * !* * !* A11 A12 A16 B11 B12 B16 * !* A12 A22 A26 B12 B22 B26 ! A B ! * !* A16 A26 A66 B16 B26 B66 or ! ! * !* B11 B12 B16 D11 D12 D16 ! B D ! * !* B12 B22 B26 D12 D22 D26 * !* B16 B26 B66 D16 D26 D66 * !* * !* The stress vector is (Nx,Ny,Nxy,Mx,My,Mxy) * !* * !* The deformation vector is (EPs xx, EPs yy, gam xy, Kx, Ky, Kxy) * !* * !* ----------------------------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* CALCULATE THE STIFFNESS MATRIX OF A LAMINATE * !* * !* Number of layers: 4 * !* * !* Layer 1: angle -30.0 deg. - thickness 1.0 mm * !* * !* Layer 2: angle 15.0 deg. - thickness 1.5 mm * !* * !* Layer 3: angle -15.0 deg. - thickness 1.5 mm * !* * !* Layer 4: angle 30.0 deg. - thickness 1.0 mm * !* * !* Material Parameters: * !* * !* El = 38.0 Et = 9.0 * !* NUlt = 0.32 Glt= 3.6 * !* * !* * !* BASIC REDUCED STIFFNESS MATRIX IN GPa: * !* * !* 0.38945E+02 0.29516E+01 0.00000E+00 * !* 0.29516E+01 0.92237E+01 0.00000E+00 * !* 0.00000E+00 0.00000E+00 0.36000E+01 * !* * !* Layer #: 1 * !* * !* REDUCED STIFFNESS MATRIX IN GPa: * !* * !* 0.26290E+02 0.81763E+01 -0.94512E+01 * !* 0.81763E+01 0.11429E+02 -0.34183E+01 * !* -0.94512E+01 -0.34183E+01 0.88247E+01 * !* * !* Layer #: 2 * !* * !* REDUCED STIFFNESS MATRIX IN GPa: * !* * !* 0.35212E+02 0.46931E+01 0.67316E+01 * !* 0.46931E+01 0.94731E+01 0.69862E+00 * !* 0.67316E+01 0.69862E+00 0.53416E+01 * !* * !* Layer #: 3 * !* * !* REDUCED STIFFNESS MATRIX IN GPa: * !* * !* 0.35212E+02 0.46931E+01 -0.67316E+01 * !* 0.46931E+01 0.94731E+01 -0.69862E+00 * !* -0.67316E+01 -0.69862E+00 0.53416E+01 * !* * !* Layer #: 4 * !* * !* REDUCED STIFFNESS MATRIX IN GPa: * !* * !* 0.26290E+02 0.81763E+01 0.94512E+01 * !* 0.81763E+01 0.11429E+02 0.34183E+01 * !* 0.94512E+01 0.34183E+01 0.88247E+01 * !* * !* A MATRIX IN 1E6 N/M: * !* * !* 0.15822E+03 0.30432E+02 0.00000E+00 * !* 0.30432E+02 0.51278E+02 0.00000E+00 * !* 0.00000E+00 0.00000E+00 0.33674E+02 * !* * !* B MATRIX IN 1E3 N: * !* * !* 0.00000E+00 0.00000E+00 0.22659E+02 * !* 0.00000E+00 0.00000E+00 0.12101E+02 * !* 0.22659E+02 0.12101E+02 0.00000E+00 * !* * !* D MATRIX IN N.DM: * !* * !* 0.29393E+03 0.77333E+02 0.00000E+00 * !* 0.77333E+02 0.11465E+03 0.00000E+00 * !* 0.00000E+00 0.00000E+00 0.84087E+02 * !* * !* STIFFNESS R MATRIX: * !* * !* 0.15822E+09 0.30432E+08 0.00000E+00 0.00000E+00 0.00000E+00 0.22659E+05 * !* 0.30432E+08 0.51278E+08 0.00000E+00 0.00000E+00 0.00000E+00 0.12101E+05 * !* 0.00000E+00 0.00000E+00 0.33674E+08 0.22659E+05 0.12101E+05 0.00000E+00 * !* 0.00000E+00 0.00000E+00 0.22659E+05 0.29393E+03 0.77333E+02 0.00000E+00 * !* 0.00000E+00 0.00000E+00 0.12101E+05 0.77333E+02 0.11465E+03 0.00000E+00 * !* 0.22659E+05 0.12101E+05 0.00000E+00 0.00000E+00 0.00000E+00 0.84087E+02 * !* * !* END OF PROGRAM. * !* ----------------------------------------------------------------------------------- * !* REFERENCE: "MATERIAUX COMPOSITES - Comportement mécanique et analyse des structures * !* By J.-M. Berthelot - MASSON 1996". * !* * !* F90 Release By J-P MOREAU - August 1997. * !* (www.jpmoreau.fr) * !*************************************************************************************** PROGRAM COMPOS03 ! MAXCOU = maximum number of layers PARAMETER (NMAX=3,MAXCOU=10) COMMON /A/ EL,ET,NULT,GLT REAL*8 EL,ET,NULT,GLT,PI ! layer angle and thickness REAL*8 TH(MAXCOU),EP(MAXCOU) REAL*8 H(MAXCOU),HM REAL*8 Q0(NMAX,NMAX),Q(MAXCOU,NMAX,NMAX) REAL*8 A(NMAX,NMAX),B(NMAX,NMAX),D(NMAX,NMAX) REAL*8 R(6,6) PI = 4.D0*DATAN(1.D0) ! material constants in GPa EL = 38.D0 ET = 9.D0 NULT = 0.32D0 GLT = 3.6D0 ! material data WRITE(*,*) ' ' WRITE(*,*) ' CALCULATE THE STIFFNESS MATRIX OF A LAMINATE' WRITE(*,*) ' ' ! number of layers NCOUCHES=4 ! define angles and thicknesses of layers 1-4 TH(1)=-30.D0 EP(1)=1.0D-3 TH(2)= 15.D0 EP(2)=1.5D-3 TH(3)=-15.D0 EP(3)=1.5D-3 TH(4)= 30.D0 EP(4)=1.0D-3 WRITE(*,*) ' Number of layers: ', NCOUCHES DO I=1,NCOUCHES WRITE(*,10) I,TH(I),1.D3*EP(I) ! put angles in radians TH(I) = TH(I)*PI/180.D0 ENDDO WRITE(*,*) ' ' WRITE(*,*) ' Material Parameters:' WRITE(*,20) EL,ET WRITE(*,21) NULT,GLT WRITE(*,*) ' ' ! reduced stiffness matrix (angle = 0) CALL MATNUL(Q0,NMAX) CALL CAL_Q0(Q0) CALL MATPR2 (' BASIC REDUCED STIFFNESS MATRIX IN GPa:',Q0,NMAX) ! calculate reduced stiffness matrices for different angles <> 0 DO I=1,NCOUCHES IF (TH(I).EQ.0.D0) THEN CALL MATEGAL(I,Q,Q0) ELSE CALL CAL_Q(I,Q0,Q,TH(I)) ENDIF WRITE(*,*) ' ' WRITE(*,*) ' Layer #:',I CALL MATPR1(' REDUCED STIFFNESS MATRIX IN GPa:',I,Q,NMAX) ENDDO ! calculate A matrix ! A = sum of 1e3*EP(k)*Q(K) CALL MATNUL(A,NMAX) DO I=1,NMAX DO J=1,NMAX DO K=1,NCOUCHES A(I,J)=A(I,J)+1000.D0*(EP(K)*Q(K,I,J)) ENDDO IF (A(I,J).LT.1.D-8) A(I,J)=0.D0 ENDDO ENDDO CALL MATPR2(' A MATRIX IN 1E6 N/M:',A,3) ! calculate B matrix ! B = sum of 1e6*(1/2)*(h(k)*h(k)-h(k-1)*h(k-1))*Q(k) CALL MATNUL(B,NMAX) HM=0.D0 DO I=1,NCOUCHES HM=HM+EP(I) ENDDO ! average height of stack HM=HM/2.D0 ! lower height layer #1 H(1) = -HM DO I=2,NCOUCHES+1 H(I)=H(I-1)+EP(I-1) ENDDO DO I=1,NMAX DO J=1,NMAX DO K=1,NCOUCHES B(I,J)=B(I,J)+1.D6*0.5D0*((H(K+1)*H(K+1)-H(K)*H(K))*Q(K,I,J)) ENDDO IF (B(I,J).LT.1.D-8) B(I,J)=0.D0 ENDDO ENDDO CALL MATPR2(' B MATRIX IN 1E3 N:',B,NMAX) ! calculate D matrix ! D = sum of 1e9*(1/3)*(h(k)*h(k)*h(k)-h(k-1)*h(k-1)*h(k-1))*Q(k) CALL MATNUL(D,NMAX) DO I=1,NMAX DO J=1,NMAX DO K=2,NCOUCHES+1 D(I,J)=D(I,J)+1.0D9*(1.D0/3.D0)*((H(K)*H(K)*H(K)-H(K-1)*H(K-1)*H(K-1))*Q(K-1,I,J)) ENDDO IF (D(I,J).LT.1.D-8) D(I,J)=0.D0 ENDDO ENDDO CALL MATPR2(' D MATRIX IN N.DM:',D,NMAX) ! assemble stiffness matrix of the laminate ! ! A B ! ! R = ! ! ! ! B D ! CALL MATNUL(R,6) DO I=1,NMAX DO J=1,NMAX R(I,J) = 1.D6*A(I,J) ENDDO ENDDO DO I=1,NMAX DO J=4,NMAX+3 R(I,J) = 1.D3*B(I,J-3) ENDDO ENDDO DO I=4,NMAX+3 DO J=1,NMAX R(I,J) = 1.D3*B(I-3,J) ENDDO ENDDO DO I=4,NMAX+3 DO J=4,NMAX+3 R(I,J) = D(I-3,J-3) ENDDO ENDDO CALL MATPR2(' STIFFNESS R MATRIX:',R,6) WRITE(*,*) ' ' STOP ' END OF PROGRAM.' 10 FORMAT(/' Layer',I3,': angle',F6.1,' deg. - thickness',F6.1,' mm') 20 FORMAT(/' El =',F6.1,' Et =',F6.1) 21 FORMAT(' NUlt =',F6.2,' Glt=',F6.1) END ! set a matrix (N,N) to zero SUBROUTINE MATNUL(A,N) REAL*8 A(N,N) DO I=1,N DO J=1,N A(I,J)=0.D0 ENDDO ENDDO END ! A(3,3) = B(3,3) SUBROUTINE MATEGAL(K,A,B) PARAMETER (N=3,NCOUCHES=10) REAL*8 A(NCOUCHES,N,N), B(N,N) DO I=1,N DO J=1,N A(K,I,J)=B(I,J) ENDDO ENDDO END ! multiply a 3x3 matrixpar by a 3x1 vector SUBROUTINE MATMULT(A, B, C) PARAMETER (N=3) REAL*8 A(N,N),B(N),C(N),SUM DO I=1,N SUM= 0.D0 DO K=1,N SUM=SUM+A(I,K)*B(K) ENDDO C(I)=SUM ENDDO END ! print a NxN matrix SUBROUTINE MATPR2(TITRE, A, N) CHARACTER*(*) TITRE REAL*8 A(N,N) WRITE(*,*) ' ' WRITE(*,*) titre WRITE(*,*) ' ' IF (N.EQ.3) THEN WRITE(*,10) ((A(I,J),J=1,N),I=1,N) ELSE WRITE(*,20) ((A(I,J),J=1,N),I=1,N) ENDIF 10 FORMAT(3(' ',E12.5)) 20 FORMAT(6(' ',E12.5)) END SUBROUTINE MATPR1(TITRE, K, A, N) PARAMETER(NCOU=10) CHARACTER*(*) TITRE REAL*8 A(NCOU,N,N) WRITE(*,*) ' ' WRITE(*,*) TITRE WRITE(*,*) ' ' WRITE(*,10) ((A(K,I,J),J=1,N),I=1,N) 10 FORMAT(3(' ',E12.5)) END ! calculate cos(t) power 4 REAL*8 FUNCTION DCOS4(T) REAL*8 A,T A=DCOS(T) DCOS4=A*A*A*A END ! calculate sin(t) power 4 REAL*8 FUNCTION DSIN4(T) REAL*8 A,T A=DSIN(T) DSIN4=A*A*A*A END ! calculate the reduced stiffness matrix for theta=0 SUBROUTINE CAL_Q0(Q) COMMON /A/ EL,ET,NULT,GLT REAL*8 EL,ET,NULT,GLT REAL*8 Q(3,3) Q(1,1) = EL/(1.D0-((ET/EL)*NULT*NULT)) Q(2,2) = (ET/EL)*Q(1,1) Q(1,2) = NULT*Q(2,2) Q(2,1) = Q(1,2) Q(3,3) = GLT END ! calculate the reduced stiffness Q matrix for theta <> 0 ! with the condition that the Q matrix for theta=0 is available ! called here Q2 SUBROUTINE CAL_Q(K, Q2, Q, TH) PARAMETER (N=3,NCOUCHES=10) REAL*8 Q(NCOUCHES,N,N),Q2(N,N),TH,TEMP REAL*8 DCOS4, DSIN4 TEMP = DSIN(TH)*DSIN(TH)*DCOS(TH)*DCOS(TH) Q(K,1,1) = Q2(1,1)*DCOS4(TH)+Q2(2,2)*DSIN4(TH)+2.D0*(Q2(1,2)+2.D0*Q2(3,3))*TEMP Q(K,1,2) = (Q2(1,1)+Q2(2,2)-4.D0*Q2(3,3))*TEMP+Q2(1,2)*(DCOS4(TH)+DSIN4(TH)) Q(K,2,1)=Q(K,1,2) TEMP = DSIN(TH)*DCOS(TH)*DCOS(TH)*DCOS(TH) Q(K,1,3)=(Q2(1,1)-Q2(1,2)-2.D0*Q2(3,3))*TEMP TEMP = DSIN(TH)*DSIN(TH)*DSIN(TH)*DCOS(TH) Q(K,1,3)=Q(K,1,3)+(Q2(1,2)-Q2(2,2)+2.D0*Q2(3,3))*TEMP Q(K,3,1)=Q(K,1,3) TEMP = DSIN(TH)*DSIN(TH)*DCOS(TH)*DCOS(TH) Q(K,2,2) = Q2(1,1)*DSIN4(TH)+2.D0*(Q2(1,2)+2.D0*Q2(3,3))*TEMP+Q2(2,2)*DCOS4(TH) TEMP = DSIN(TH)*DSIN(TH)*DSIN(TH)*DCOS(TH) Q(K,2,3)=(Q2(1,1)-Q2(1,2)-2.D0*Q2(3,3))*TEMP TEMP = DCOS(TH)*DCOS(TH)*DCOS(TH)*DSIN(TH) Q(K,2,3)=Q(K,2,3)+(Q2(1,2)-Q2(2,2)+2.D0*Q2(3,3))*TEMP Q(K,3,2)=Q(K,2,3) TEMP = DSIN(TH)*DSIN(TH)*DCOS(TH)*DCOS(TH) Q(K,3,3) = (Q2(1,1)+Q2(2,2)-2.D0*(Q2(1,2)+Q2(3,3)))*TEMP+Q2(3,3)*(DSIN4(TH)+DCOS4(TH)) END ! end of file compos3.f90