'*************************************************************************************** '* This program calculates the deformations and stresses in a laminated material made * '* of n unidirectional composite layers, knowing the resulting imposed efforts. * '* * '* 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 imposed efforts vector is (Nx,Ny,Nxy,Mx,My,Mxy) * '* * '* The deformation vector is (EPs xx, EPs yy, gam xy, Kx, Ky, Kxy) * '* * '* ----------------------------------------------------------------------------------- * '* SAMPLE RUN: * '* The output file Compos04.lst contains (extract): * '* * '* CALCULATE DEFORMATIONS AND STRESSES * '* IN A LAMINATED MATERIAL * '* * '* Imposed resulting efforts: * '* * '* NX =0.10000E+07 NY =0.50000E+06 NXY =0.25000E+06 * '* MX =0.00000E+00 MY =0.00000E+00 MXY =0.00000E+00 * '* * '* Number of layers: 4 * '* * '* Layer 1 : angle 15.0 deg. - thickness 1.5 mm * '* * '* Layer 2 : angle -30.0 deg. - thickness 1.0 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.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 * '* ---/--- * '* Stresses in each layer in main axes: * '* * '* Layer # 1 * '* * '* For z = -0.25000E+01 mm: * '* * '* sigma L : 0.28274E+03 MPa * '* sigma T : 0.78080E+02 MPa * '* sigma LT: 0.45159E+02 MPa * '* * '* For z = -0.10000E+01 mm: * '* * '* sigma L : 0.28858E+03 MPa * '* sigma T : 0.70533E+02 MPa * '* sigma LT: 0.34701E+02 MPa * '* * '* Layer # 2 * '* * '* For z = -0.10000E+01 mm: * '* * '* sigma L : 0.86434E+02 MPa * '* sigma T : 0.10576E+03 MPa * '* sigma LT: 0.57365E+01 MPa * '* * '* For z = 0.00000E+00 mm: * '* * '* sigma L : 0.11193E+03 MPa * '* sigma T : 0.96965E+02 MPa * '* sigma LT: 0.83890E+01 MPa * '* * '* Layer # 3 * '* * '* For z = 0.00000E+00 mm: * '* * '* sigma L : 0.15147E+03 MPa * '* sigma T : 0.90075E+02 MPa * '* sigma LT: 0.21129E+02 MPa * '* * '* For z = 0.15000E+01 mm: * '* * '* sigma L : 0.19264E+03 MPa * '* sigma T : 0.76371E+02 MPa * '* sigma LT: 0.19346E+02 MPa * '* * '* Layer # 4 * '* * '* For z = 0.15000E+01 mm: * '* * '* sigma L : 0.33321E+03 MPa * '* sigma T : 0.51876E+02 MPa * '* sigma LT: 0.87729E+01 MPa * '* * '* For z = 0.25000E+01 mm: * '* * '* sigma L : 0.31791E+03 MPa * '* sigma T : 0.50191E+02 MPa * '* sigma LT: 0.14086E+01 MPa * '* ----------------------------------------------------------------------------------- * '* 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 COMPOS04 DEFDBL A-H, O-Z DEFINT I-N ' MAXCOU = maximum number of layers NMAX = 3 MAXCOU = 10 DIM NULT AS DOUBLE ' angle and thickness of layer 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), X(6, 6), E0(6), EDEB(6), EFIN(6) DIM N(6) AS DOUBLE DIM NX AS DOUBLE DIM NY AS DOUBLE DIM NXY AS DOUBLE DIM MX AS DOUBLE DIM MY AS DOUBLE DIM MXY AS DOUBLE DIM T(NMAX, NMAX), AK(MAXCOU, NMAX), BK(MAXCOU, NMAX) DIM AK1(MAXCOU, NMAX), BK1(MAXCOU, NMAX) DIM TEMP(NMAX), TEMP1(NMAX) DIM Title AS STRING PI = 4# * ATN(1#) 'Output formats F\$ = "####.####" F1\$ = "##########.##" F2\$ = "####.#### mm" F3\$ = "####.#### MPa" ' Material constants in GPa EL = 38# ET = 9# NULT = .32# GLT = 3.6# ' imposed efforts NX = 1000000# NY = 500000# NXY = 250000# MX = 0# MY = 0# MXY = 0# N(1) = .000000001# * NX N(2) = .000000001# * NY N(3) = .000000001# * NXY N(4) = .000000001# * MX N(5) = .000000001# * MY N(6) = .000000001# * MXY FOR I = 1 TO 6 E0(I) = 0# EDEB(I) = 0# EFIN(I) = 0# NEXT I ' print data to output file OPEN "compos04.lst" FOR OUTPUT AS #2 CLS PRINT #2, PRINT #2, " CALCULATE DEFORMATIONS AND STRESSES" PRINT #2, " IN A LAMINATED MATERIAL" PRINT #2, PRINT #2, " Imposed resulting efforts:" PRINT #2, PRINT #2, " NX="; NX; " NY = "; NY; " NXY = "; NXY PRINT #2, " MX="; MX; " MY = "; MY; " MXY = "; MXY ' number of layers in laminate NCOUCHES = 4 ' define angles and thicknesses of layers 1-4 TH(1) = 15# EP(1) = .0015# TH(2) = -30# EP(2) = .001# TH(3) = -15# EP(3) = .0015# TH(4) = 30# EP(4) = .001# PRINT #2, PRINT #2, " Number of layers: "; NCOUCHES PRINT #2, FOR I = 1 TO NCOUCHES PRINT #2, " Layer #"; I; ": angle "; TH(I); " deg. - thickness "; 1000# * EP(I); " mm" PRINT #2, ' put angles in radians TH(I) = TH(I) * PI / 180# NEXT I PRINT #2, PRINT #2, " Material Parameters:" PRINT #2, PRINT #2, " EL="; EL; " ET="; ET PRINT #2, " NULT="; NULT; " GLT="; GLT PRINT #2, ' calculate the 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 the 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 #2, " 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 #2, ' Direct inversion of R matrix, stored in R N = 6: M = 0 GOSUB 3000 'CALL MATINV(6,0,R,X,DET) PRINT #2, PRINT #2, " DETERMINANT OF R: "; DET FOR I = 1 TO 6 FOR J = 1 TO 6 R(I, J) = 1000000000# * R(I, J) NEXT J NEXT I Title = " INVERSE MATRIX OF R:" GOSUB 620 'CALL MATPR2(" INVERSE MATRIX OF R:", R, 6) ' Calculate deformations eps0xx, eps0yy, gam0xy,kx,ky,kxy ' in main axes, stored in vector E0(6) GOSUB 520 'CALL MATMUL(R,N,E0,6) PRINT #2, PRINT #2, " Membrane deformations (in mm) and curbatures:" PRINT #2, PRINT #2, " eps0xx: "; : PRINT #2, USING F\$; 1000# * E0(1) PRINT #2, " eps0yy: "; : PRINT #2, USING F\$; 1000# * E0(2) PRINT #2, " gam0xy: "; : PRINT #2, USING F\$; 1000# * E0(3) PRINT #2, " kx : "; : PRINT #2, USING F\$; E0(4) PRINT #2, " ky : "; : PRINT #2, USING F\$; E0(5) PRINT #2, " kxy : "; : PRINT #2, USING F\$; E0(6) ' Deformations in reference axes (x,y) ' ! epsxx ! ! eps0xx ! ! kx ! ' ! epsyy ! = ! eps0yy ! + z ! ky ! ' ! gamxy ! ! gam0xy ! ! kxy ! PRINT #2, PRINT #2, " Deformations in reference axes (x,y):" PRINT #2, " For z = "; 1000# * H(1); " mm:" FOR I = 1 TO NMAX EDEB(I) = E0(I) + H(1) * E0(I + 3) NEXT I PRINT #2, " epsxx: "; : PRINT #2, USING F2\$; 1000# * EDEB(1) PRINT #2, " epsyy: "; : PRINT #2, USING F2\$; 1000# * EDEB(2) PRINT #2, " gamxy: "; : PRINT #2, USING F2\$; 1000# * EDEB(3) PRINT #2, " For z = "; 1000# * H(NCOUCHES + 1); " mm:" FOR I = 1 TO NMAX EFIN(I) = E0(I) + H(NCOUCHES + 1) * E0(I + 3) NEXT I PRINT #2, " epsxx: "; : PRINT #2, USING F2\$; 1000# * EFIN(1) PRINT #2, " epsyy: "; : PRINT #2, USING F2\$; 1000# * EFIN(2) PRINT #2, " gamxy: "; : PRINT #2, USING F2\$; 1000# * EFIN(3) ' Deformations in each layer in main axes of layer PRINT #2, PRINT #2, " Deformations in each layer in main axes:" FOR K = 1 TO NCOUCHES TH = TH(K) GOSUB 4000 'CALL INITT(T,TH(K)) FOR I = 1 TO NMAX TEMP(I) = E0(I) NEXT I GOSUB 530 'CALL MATMUL(T,TEMP,TEMP1,3) FOR I = 1 TO NMAX AK(K, I) = TEMP1(I) NEXT I FOR I = 1 TO NMAX TEMP(I) = E0(I + 3) NEXT I GOSUB 530 'CALL MATMUL(T,TEMP,TEMP1,3) FOR I = 1 TO NMAX BK(K, I) = TEMP1(I) NEXT I PRINT #2, PRINT #2, " Layer #"; K; ":" PRINT #2, " For z = "; 1000# * H(K); " mm:" FOR I = 1 TO NMAX EDEB(I) = AK(K, I) + H(K) * BK(K, I) NEXT I PRINT #2, " epsxx: "; : PRINT #2, USING F2\$; 1000# * EDEB(1) PRINT #2, " epsyy: "; : PRINT #2, USING F2\$; 1000# * EDEB(2) PRINT #2, " gamxy: "; : PRINT #2, USING F2\$; 1000# * EDEB(3) PRINT #2, " For z = "; 1000# * H(K + 1); " mm:" FOR I = 1 TO NMAX EFIN(I) = AK(K, I) + H(K + 1) * BK(K, I) NEXT I PRINT #2, " epsxx: "; : PRINT #2, USING F2\$; 1000# * EFIN(1) PRINT #2, " epsyy: "; : PRINT #2, USING F2\$; 1000# * EFIN(2) PRINT #2, " gamxy: "; : PRINT #2, USING F2\$; 1000# * EFIN(3) NEXT K ' Stresses in each layer in reference axes PRINT #2, PRINT #2, " Stresses in each layer in reference axes:" FOR K = 1 TO NCOUCHES FOR I = 1 TO NMAX TEMP(I) = E0(I) NEXT I FOR I = 1 TO NMAX FOR J = 1 TO NMAX A(I, J) = Q(K, I, J) NEXT J NEXT I GOSUB 531 'CALL MATMUL(A,TEMP,TEMP1,3) FOR I = 1 TO NMAX AK(K, I) = TEMP1(I) NEXT I FOR I = 1 TO NMAX TEMP(I) = E0(I + 3) NEXT I GOSUB 531 'CALL MATMUL(A,TEMP,TEMP1,3) FOR I = 1 TO NMAX BK(K, I) = TEMP1(I) NEXT I PRINT #2, PRINT #2, " Layer #"; K; ":" PRINT #2, " For z = "; 1000# * H(K); " mm:" FOR I = 1 TO NMAX EDEB(I) = 1000# * (AK(K, I) + H(K) * BK(K, I)) NEXT I PRINT #2, " sigma L : "; : PRINT #2, USING F3\$; EDEB(1) PRINT #2, " sigma T : "; : PRINT #2, USING F3\$; EDEB(2) PRINT #2, " sigma LT: "; : PRINT #2, USING F3\$; EDEB(3) PRINT #2, " For z = "; 1000# * H(K + 1); " mm:" FOR I = 1 TO NMAX EFIN(I) = 1000# * (AK(K, I) + H(K + 1) * BK(K, I)) NEXT I PRINT #2, " sigma L : "; : PRINT #2, USING F3\$; EFIN(1) PRINT #2, " sigma T : "; : PRINT #2, USING F3\$; EFIN(2) PRINT #2, " sigma LT: "; : PRINT #2, USING F3\$; EFIN(3) NEXT K ' Stresses in each layer in main axes PRINT #2, PRINT #2, " Stresses in each layer in main axes:" FOR K = 1 TO NCOUCHES TH = TH(K) GOSUB 4100 'CALL INITT1(T,TH(K)) FOR I = 1 TO NMAX TEMP(I) = AK(K, I) NEXT I GOSUB 530 'CALL MATMUL(T,TEMP,TEMP1,3) FOR I = 1 TO NMAX AK1(K, I) = TEMP1(I) NEXT I FOR I = 1 TO NMAX TEMP(I) = BK(K, I) NEXT I GOSUB 530 'CALL MATMUL(T,TEMP,TEMP1,3) FOR I = 1 TO NMAX BK1(K, I) = TEMP1(I) NEXT I PRINT #2, PRINT #2, " Layer #"; K; ":" PRINT #2, " For z = "; 1000# * H(K); " mm:" FOR I = 1 TO NMAX EDEB(I) = 1000# * (AK1(K, I) + H(K) * BK1(K, I)) NEXT I PRINT #2, " sigma L : "; : PRINT #2, USING F3\$; EDEB(1) PRINT #2, " sigma T : "; : PRINT #2, USING F3\$; EDEB(2) PRINT #2, " sigma LT: "; : PRINT #2, USING F3\$; EDEB(3) PRINT #2, " For z = "; 1000# * H(K + 1); " mm:" FOR I = 1 TO NMAX EFIN(I) = 1000# * (AK1(K, I) + H(K + 1) * BK1(K, I)) NEXT I PRINT #2, " sigma L : "; : PRINT #2, USING F3\$; EFIN(1) PRINT #2, " sigma T : "; : PRINT #2, USING F3\$; EFIN(2) PRINT #2, " sigma LT: "; : PRINT #2, USING F3\$; EFIN(3) NEXT K PRINT #2, CLOSE #2 PRINT " END OF PROGRAM." END 'of main program 'UTILITY ROUTINES ' 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 ' multiply a 6x6 matrix by a 6x1 vector 520 'SUBROUTINE MATMUL(R, N, E0, 6) FOR I = 1 TO 6 SUM = 0# FOR K = 1 TO 6 SUM = SUM + R(I, K) * N(K) NEXT K E0(I) = SUM NEXT I RETURN ' multiply a 3x3 matrix by a 3x1 vector 530 'SUBROUTINE MATMUL(T, TEMP, TEMP1, 3) FOR I = 1 TO 3 SUM = 0# FOR K1 = 1 TO 3 SUM = SUM + T(I, K1) * TEMP(K1) NEXT K1 TEMP1(I) = SUM NEXT I RETURN 531 'SUBROUTINE MATMUL(A, TEMP, TEMP1, 3) FOR I = 1 TO 3 SUM = 0# FOR K1 = 1 TO 3 SUM = SUM + A(I, K1) * TEMP(K1) NEXT K1 TEMP1(I) = SUM NEXT I RETURN ' print Q0(3,3) matrix 600 'SUBROUTINE MATPR2(TITLE, Q0, NMAX) PRINT #2, PRINT #2, Title PRINT #2, FOR I = 1 TO NMAX FOR J = 1 TO NMAX PRINT #2, USING F\$; Q0(I, J); NEXT J PRINT #2, NEXT I RETURN ' print A(3,3) matrix 601 'SUBROUTINE MATPR2(TITLE, A, NMAX) PRINT #2, PRINT #2, Title PRINT #2, FOR I = 1 TO NMAX FOR J = 1 TO NMAX PRINT #2, USING F\$; A(I, J); NEXT J PRINT #2, NEXT I RETURN ' print B(3,3) matrix 602 'SUBROUTINE MATPR2(TITLE, B, NMAX) PRINT #2, PRINT #2, Title PRINT #2, FOR I = 1 TO NMAX FOR J = 1 TO NMAX PRINT #2, USING F\$; B(I, J); NEXT J PRINT #2, NEXT I RETURN ' print D(3,3) matrix 603 'SUBROUTINE MATPR2(TITLE, B, NMAX) PRINT #2, PRINT #2, Title PRINT #2, FOR I = 1 TO NMAX FOR J = 1 TO NMAX PRINT #2, USING F\$; D(I, J); NEXT J PRINT #2, NEXT I RETURN ' print Q(K,3,3) matrix 610 'SUBROUTINE MATPR1(TITLE, K, Q, NMAX) PRINT #2, PRINT #2, Title PRINT #2, FOR II = 1 TO NMAX FOR JJ = 1 TO NMAX PRINT #2, USING F\$; Q(K, II, JJ); NEXT JJ PRINT #2, NEXT II RETURN ' print R(6,6) matrix 620 'SUBROUTINE MATPR2(TITLE, R, 6) PRINT #2, PRINT #2, Title PRINT #2, FOR I = 1 TO 6 FOR J = 1 TO 6 PRINT #2, USING F1\$; R(I, J); NEXT J PRINT #2, 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 3000 'SUB MATINV(N,M,R,X,DET) '******************************************* '* SOLVING A LINEAR MATRIX SYSTEM AX = B * '* with Gauss-Jordan method using full * '* pivoting at each step. During the pro- * '* cess, original A and B matrices are * '* destroyed to spare storage location. * '* --------------------------------------- * '* INPUTS: A MATRIX N*N * '* B MATRIX N*M * '* --------------------------------------- * '* OUTPUS: A INVERSE OF A N*N * '* DET DETERMINANT OF A * '* B SOLUTION MATRIX N*M * '* --------------------------------------- * '* NOTA - If M=0 inversion of A matrix * '* only. * '******************************************* DIM PC(N) AS INTEGER DIM PL(N) AS INTEGER DIM CS(N) EPSMACH = .0000000001# ' INITIALIZATIONS: DET = 1# FOR I = 1 TO N PC(I) = 0 PL(I) = 0 CS(I) = 0# NEXT I ' Main loop FOR K = 1 TO N ' Searching greatest pivot: PV = R(K, K) IK = K JK = K PAV = ABS(PV) FOR I = K TO N FOR J = K TO N IF ABS(R(I, J)) > PAV THEN PV = R(I, J) PAV = ABS(PV) IK = I JK = J END IF NEXT J NEXT I ' Search terminated, Pivot is PV at location I=IK, J=JK. ' Save Pivot location: PC(K) = JK PL(K) = IK 'Determinant DET is actualised 'If DET=0, ERROR MESSAGE and STOP 'Machine dependant EPSMACH equals here 1D-10 IF IK <> K THEN DET = -DET IF JK <> K THEN DET = -DET DET = DET * PV IF ABS(DET) < EPSMACH THEN ' ERROR MESSAGE AND STOP PRINT " STOP IN MATINV: DETERMINANT TOO SMALL!" END END IF ' Locate Pivot in K,K: IF IK <> K THEN FOR I = 1 TO N ' EXCHANGE LINES IK and K of R: TT = R(IK, I) R(IK, I) = R(K, I) R(K, I) = TT NEXT I END IF IF M <> 0 THEN FOR I = 1 TO M TT = X(IK, I) X(IK, I) = X(K, I) X(K, I) = TT NEXT I END IF ' Pivot is located at right line. IF JK <> K THEN FOR I = 1 TO N ' EXCHANGE COLUMNS JK and K of R : TT = R(I, JK) R(I, JK) = R(I, K) R(I, K) = TT NEXT I END IF ' Pivot is located at right column. ' Pivot is now located at K,K. ' STORE COLUMN K OF R IN VECTOR CS ' ET PUT COLUMN TO ZERO FOR I = 1 TO N CS(I) = R(I, K) R(I, K) = 0# NEXT I CS(K) = 0# R(K, K) = 1# ' MODIFY LINE K OF R: IF ABS(PV) < EPSMACH THEN PRINT #2, " STOP MATINV: PIVOT TOO SMALL!" END END IF FOR I = 1 TO N R(K, I) = R(K, I) / PV NEXT I IF M <> 0 THEN FOR I = 1 TO M X(K, I) = X(K, I) / PV NEXT I END IF ' MODIFY OTHER LINES OF R: FOR J = 1 TO N IF J = K THEN GOTO 20 FOR I = 1 TO N ' MODIFY LINE J OF R : R(J, I) = R(J, I) - CS(J) * R(K, I) NEXT I IF M <> 0 THEN FOR I = 1 TO M X(J, I) = X(J, I) - CS(J) * X(K, I) NEXT I END IF 20 NEXT J ' LINE K OVER. NEXT K 'End K loop ' MATRIX IS NOW INVERSED - MATRIX MUST BE REARRANGED ' EXCHANGE LINES FOR I = N TO 1 STEP -1 IK = PC(I) IF IK = I THEN GOTO 30 ' EXCHANGE LINES I and PC(I) of R: FOR J = 1 TO N TT = R(I, J) R(I, J) = R(IK, J) R(IK, J) = TT NEXT J IF M <> 0 THEN FOR J = 1 TO M TT = X(I, J) X(I, J) = X(IK, J) X(IK, J) = TT NEXT J END IF ' NO MORE EXCHANGE NECESSARY, GO TO NEXT LINE 30 NEXT I ' EXCHANGE COLUMNS FOR J = N TO 1 STEP -1 JK = PL(J) IF JK = J THEN GOTO 40 ' EXCHANGE COLUMNS J and PL(J) of R: FOR I = 1 TO N TT = R(I, J) R(I, J) = R(I, JK) R(I, JK) = TT NEXT I ' NO MORE EXCHANGE NECESSARY, GO TO NEXT COLUMN 40 NEXT J ' REARRANGEMENT TERMINATED. RETURN 'MATINV ' Matrix to change axes (x,y) for deformations 4000 'SUB INITT(T,TH) T(1, 1) = COS(TH) * COS(TH) T(1, 2) = SIN(TH) * SIN(TH) T(1, 3) = SIN(TH) * COS(TH) T(2, 1) = SIN(TH) * SIN(TH) T(2, 2) = COS(TH) * COS(TH) T(2, 3) = -SIN(TH) * COS(TH) T(3, 1) = -2# * SIN(TH) * COS(TH) T(3, 2) = 2# * SIN(TH) * COS(TH) T(3, 3) = COS(TH) * COS(TH) - SIN(TH) * SIN(TH) RETURN ' Matrix to change axes (x,y) for stresses 4100 'SUB INITT1(T,TH) T(1, 1) = COS(TH) * COS(TH) T(1, 2) = SIN(TH) * SIN(TH) T(1, 3) = 2# * SIN(TH) * COS(TH) T(2, 1) = SIN(TH) * SIN(TH) T(2, 2) = COS(TH) * COS(TH) T(2, 3) = -2# * SIN(TH) * COS(TH) T(3, 1) = -SIN(TH) * COS(TH) T(3, 2) = SIN(TH) * COS(TH) T(3, 3) = COS(TH) * COS(TH) - SIN(TH) * SIN(TH) RETURN ' End of file Compos04.bas