'*************************************************************************************** '* 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 deg. - thickness 1 mm * '* * '* Layer # 2 : angle 15 deg. - thickness 1.5 mm * '* * '* Layer # 3 : angle -15 deg. - thickness 1.5 mm * '* * '* Layer # 4 : angle 30 deg. - thickness 1 mm * '* * '* * '* Material Parameters: * '* EL= 38 ET= 9 * '* NULT= .32 GLT= 3.6 * '* * '* * '* BASIC REDUCED STIFFNESS MATRIX IN GPa: * '* * '* 38.9445 2.9516 0.0000 * '* 2.9516 9.2237 0.0000 * '* 0.0000 0.0000 3.6000 * '* * '* Layer # 1 : * '* * '* REDUCED STIFFNESS MATRIX IN GPa: * '* * '* 26.2896 8.1763 -9.4512 * '* 8.1763 11.4292 -3.4183 * '* -9.4512 -3.4183 8.8247 * '* * '* Layer # 2 : * '* * '* REDUCED STIFFNESS MATRIX IN GPa: * '* * '* 35.2120 4.6931 6.7316 * '* 4.6931 9.4731 0.6986 * '* 6.7316 0.6986 5.3416 * '* * '* Layer # 3 : * '* * '* REDUCED STIFFNESS MATRIX IN GPa: * '* * '* 35.2120 4.6931 -6.7316 * '* 4.6931 9.4731 -0.6986 * '* -6.7316 -0.6986 5.3416 * '* * '* Layer # 4 : * '* * '* REDUCED STIFFNESS MATRIX IN GPa: * '* * '* 26.2896 8.1763 9.4512 * '* 8.1763 11.4292 3.4183 * '* 9.4512 3.4183 8.8247 * '* * '* A MATRIX IN 1E6 N/M: * '* * '* 158.2153 30.4320 0.0000 * '* 30.4320 51.2776 0.0000 * '* 0.0000 0.0000 33.6741 * '* * '* B MATRIX IN 1E3 N: * '* * '* 0.0000 0.0000 22.6588 * '* 0.0000 0.0000 12.1012 * '* 22.6588 12.1012 0.0000 * '* * '* D MATRIX IN N.DM: * '* * '* 293.9255 77.3325 0.0000 * '* 77.3325 114.6529 0.0000 * '* 0.0000 0.0000 84.0869 * '* * '* STIFFNESS R MATRIX: * '* * '* 158215296.03 30432002.46 0.00 0.00 0.00 22658.84 * '* 30432002.46 51277564.50 0.00 0.00 0.00 12101.16 * '* 0.00 0.00 33674084.10 22658.84 12101.16 0.00 * '* 0.00 0.00 22658.84 293.93 77.33 0.00 * '* 0.00 0.00 12101.16 77.33 114.65 0.00 * '* 22658.84 12101.16 0.00 0.00 0.00 84.09 * '* * '* END OF PROGRAM. * '* * '* ----------------------------------------------------------------------------------- * '* REFERENCE: "MATERIAUX COMPOSITES - Comportement mécanique et analyse des structures * '* By J.-M. Berthelot - MASSON 1996". * '* * '* Basic Release By J-P MOREAU - August 1997. * '* (www.jpmoreau.fr) * '*************************************************************************************** ' PROGRAM COMPOS03 DEFDBL A-H, O-Z DEFINT I-N ' MAXCOU = maximum number of layers NMAX = 3 MAXCOU = 10 DIM NULT AS DOUBLE ' layer angle and thickness DIM TH(MAXCOU), EP(MAXCOU) DIM H(MAXCOU) DIM Q0(NMAX, NMAX), Q(MAXCOU, NMAX, NMAX) DIM A(NMAX, NMAX), B(NMAX, NMAX), D(NMAX, NMAX) DIM R(6, 6) DIM Title AS STRING PI = 4# * ATN(1#) F\$ = " ####.####" F1\$ = "##########.##" ' material constants in GPa EL = 38# ET = 9# NULT = .32# GLT = 3.6# ' material data PRINT PRINT " CALCULATE THE STIFFNESS MATRIX OF A LAMINATE" PRINT ' number of layers NCOUCHES = 4 ' define angles and thicknesses of layers 1-4 TH(1) = -30# EP(1) = .001# TH(2) = 15# EP(2) = .0015# TH(3) = -15# EP(3) = .0015# TH(4) = 30# EP(4) = .001# PRINT " Number of layers: "; NCOUCHES PRINT FOR I = 1 TO NCOUCHES PRINT " Layer #"; I; ": angle "; TH(I); " deg. - thickness "; 1000# * EP(I); " mm" PRINT ' put angles in radians TH(I) = TH(I) * PI / 180# NEXT I PRINT PRINT " Material Parameters:" PRINT " EL= "; EL; " ET= "; ET PRINT " NULT= "; NULT; " GLT= "; GLT PRINT ' reduced stiffness matrix (angle = 0) GOSUB 500 'CALL MATNUL(Q0,NMAX) GOSUB 1000 'CALL CAL_Q0(Q0) Title = " BASIC REDUCED STIFFNESS MATRIX IN GPa:" GOSUB 600 'CALL MATPR2 (" BASIC REDUCED STIFFNESS MATRIX IN GPa:",Q0,NMAX) ' calculate reduced stiffness matrices for different angles <> 0 FOR L = 1 TO NCOUCHES IF TH(L) = 0# THEN K = L GOSUB 510 'CALL MATEGAL(K,Q,Q0) ELSE K = L TH = TH(L) GOSUB 2000 'CALL CAL_Q(K,Q0,Q,TH) END IF PRINT PRINT " Layer #"; L; ":" Title = " REDUCED STIFFNESS MATRIX IN GPa:": K = L GOSUB 610 'CALL MATPR1(" REDUCED STIFFNESS MATRIX IN GPa:",K,Q,NMAX) NEXT L ' calculate A matrix ' A = sum of 1e3*EP(k)*Q(K) 'CALL MATNUL(A,NMAX) FOR I = 1 TO NMAX FOR J = 1 TO NMAX A(I, J) = 0# NEXT J NEXT I FOR I = 1 TO NMAX FOR J = 1 TO NMAX FOR K = 1 TO NCOUCHES A(I, J) = A(I, J) + 1000# * (EP(K) * Q(K, I, J)) NEXT K IF ABS(A(I, J)) < .00000001# THEN A(I, J) = 0# NEXT J NEXT I Title = " A MATRIX IN 1E6 N/M:" GOSUB 601 '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) FOR I = 1 TO NMAX FOR J = 1 TO NMAX B(I, J) = 0# NEXT J NEXT I HM = 0# FOR I = 1 TO NCOUCHES HM = HM + EP(I) NEXT I ' average height of stack HM = HM / 2# ' lower height layer #1 H(1) = -HM FOR I = 2 TO NCOUCHES + 1 H(I) = H(I - 1) + EP(I - 1) NEXT I FOR I = 1 TO NMAX FOR J = 1 TO NMAX FOR K = 1 TO NCOUCHES B(I, J) = B(I, J) + 1000000# * .5 * ((H(K + 1) * H(K + 1) - H(K) * H(K)) * Q(K, I, J)) NEXT K IF ABS(B(I, J)) < .00000001# THEN B(I, J) = 0# NEXT J NEXT I Title = " B MATRIX IN 1E3 N:" GOSUB 602 '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) FOR I = 1 TO NMAX FOR J = 1 TO NMAX D(I, J) = 0# NEXT J NEXT I FOR I = 1 TO NMAX FOR J = 1 TO NMAX FOR K = 2 TO NCOUCHES + 1 D(I, J) = D(I, J) + 1000000000# * (1# / 3#) * ((H(K) * H(K) * H(K) - H(K - 1) * H(K - 1) * H(K - 1)) * Q(K - 1, I, J)) NEXT K IF ABS(D(I, J)) < .00000001# THEN D(I, J) = 0# NEXT J NEXT I Title = " D MATRIX IN N.DM:" GOSUB 603 'CALL MATPR2(" D MATRIX IN N.DM:",D,NMAX) ' assemble stiffness matrix of the laminate ' ! A B ! ' R = ! ! ' ! B D ! 'CALL MATNUL(R,6) FOR I = 1 TO 6 FOR J = 1 TO 6 R(I, J) = 0# NEXT J NEXT I FOR I = 1 TO NMAX FOR J = 1 TO NMAX R(I, J) = 1000000# * A(I, J) NEXT J NEXT I FOR I = 1 TO NMAX FOR J = 4 TO NMAX + 3 R(I, J) = 1000# * B(I, J - 3) NEXT J NEXT I FOR I = 4 TO NMAX + 3 FOR J = 1 TO NMAX R(I, J) = 1000# * B(I - 3, J) NEXT J NEXT I FOR I = 4 TO NMAX + 3 FOR J = 4 TO NMAX + 3 R(I, J) = D(I - 3, J - 3) NEXT J NEXT I Title = " STIFFNESS R MATRIX:" GOSUB 620 'CALL MATPR2(" STIFFNESS R MATRIX:",R,6) PRINT PRINT " END OF PROGRAM." END 'of mainprogram ' set matrix Q0(N,N) to zero 500 'SUBROUTINE MATNUL(Q0,NMAX) FOR I = 1 TO NMAX FOR J = 1 TO NMAX Q0(I, J) = 0# NEXT J NEXT I RETURN ' A(K,3,3) = B(3,3) 510 'SUBROUTINE MATEGAL(K,Q,Q0) FOR II = 1 TO NMAX FOR JJ = 1 TO NMAX Q(K, II, JJ) = Q0(II, JJ) NEXT JJ NEXT II RETURN ' print Q0(3,3) matrix 600 'SUBROUTINE MATPR2(TITLE, Q0, NMAX) PRINT PRINT Title PRINT FOR I = 1 TO NMAX FOR J = 1 TO NMAX PRINT USING F\$; Q0(I, J); NEXT J PRINT NEXT I RETURN ' print A(3,3) matrix 601 'SUBROUTINE MATPR2(TITLE, A, NMAX) PRINT PRINT Title PRINT FOR I = 1 TO NMAX FOR J = 1 TO NMAX PRINT USING F\$; A(I, J); NEXT J PRINT NEXT I RETURN ' print B(3,3) matrix 602 'SUBROUTINE MATPR2(TITLE, B, NMAX) PRINT PRINT Title PRINT FOR I = 1 TO NMAX FOR J = 1 TO NMAX PRINT USING F\$; B(I, J); NEXT J PRINT NEXT I RETURN ' print D(3,3) matrix 603 'SUBROUTINE MATPR2(TITLE, B, NMAX) PRINT PRINT Title PRINT FOR I = 1 TO NMAX FOR J = 1 TO NMAX PRINT USING F\$; D(I, J); NEXT J PRINT NEXT I RETURN ' print Q(K,3,3) matrix 610 'SUBROUTINE MATPR1(TITLE, K, Q, NMAX) PRINT PRINT Title PRINT FOR II = 1 TO NMAX FOR JJ = 1 TO NMAX PRINT USING F\$; Q(K, II, JJ); NEXT JJ PRINT NEXT II RETURN ' print R(6,6) matrix 620 'SUBROUTINE MATPR2(TITLE, R, 6) PRINT PRINT Title PRINT FOR I = 1 TO 6 FOR J = 1 TO 6 PRINT USING F1\$; R(I, J); NEXT J PRINT NEXT I RETURN ' calculate cos(t) power 4 700 'REAL*8 FUNCTION DCOS4(T) A = COS(TH) DCOS4 = A * A * A * A RETURN ' calculate sin(t) power 4 701 'REAL*8 FUNCTION DSIN4(T) A = SIN(TH) DSIN4 = A * A * A * A RETURN ' calculate the reduced stiffness matrix for theta=0 1000 'SUBROUTINE CAL_Q0(Q0) Q0(1, 1) = EL / (1# - ((ET / EL) * NULT * NULT)) Q0(2, 2) = (ET / EL) * Q0(1, 1) Q0(1, 2) = NULT * Q0(2, 2) Q0(2, 1) = Q0(1, 2) Q0(3, 3) = GLT RETURN ' calculate the reduced stiffness Q matrix for theta <> 0 ' with the condition that the Q matrix for theta=0 is available ' called here Q2 2000 'SUBROUTINE CAL_Q(K, Q0, Q, TH) GOSUB 700: GOSUB 701 'calculate dcos4, dsin4 TEMP = SIN(TH) * SIN(TH) * COS(TH) * COS(TH) Q(K, 1, 1) = Q0(1, 1) * DCOS4 + Q0(2, 2) * DSIN4 + 2# * (Q0(1, 2) + 2# * Q0(3, 3)) * TEMP Q(K, 1, 2) = (Q0(1, 1) + Q0(2, 2) - 4# * Q0(3, 3)) * TEMP + Q0(1, 2) * (DCOS4 + DSIN4) Q(K, 2, 1) = Q(K, 1, 2) TEMP = SIN(TH) * COS(TH) * COS(TH) * COS(TH) Q(K, 1, 3) = (Q0(1, 1) - Q0(1, 2) - 2# * Q0(3, 3)) * TEMP TEMP = SIN(TH) * SIN(TH) * SIN(TH) * COS(TH) Q(K, 1, 3) = Q(K, 1, 3) + (Q0(1, 2) - Q0(2, 2) + 2# * Q0(3, 3)) * TEMP Q(K, 3, 1) = Q(K, 1, 3) TEMP = SIN(TH) * SIN(TH) * COS(TH) * COS(TH) Q(K, 2, 2) = Q0(1, 1) * DSIN4 + 2# * (Q0(1, 2) + 2# * Q0(3, 3)) * TEMP + Q0(2, 2) * DCOS4 TEMP = SIN(TH) * SIN(TH) * SIN(TH) * COS(TH) Q(K, 2, 3) = (Q0(1, 1) - Q0(1, 2) - 2# * Q0(3, 3)) * TEMP TEMP = COS(TH) * COS(TH) * COS(TH) * SIN(TH) Q(K, 2, 3) = Q(K, 2, 3) + (Q0(1, 2) - Q0(2, 2) + 2# * Q0(3, 3)) * TEMP Q(K, 3, 2) = Q(K, 2, 3) TEMP = SIN(TH) * SIN(TH) * COS(TH) * COS(TH) Q(K, 3, 3) = (Q0(1, 1) + Q0(2, 2) - 2# * (Q0(1, 2) + Q0(3, 3))) * TEMP + Q0(3, 3) * (DSIN4 + DCOS4) RETURN ' End of file compos03.bas