/************************************************************************************** * 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: * * 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: * * 1.5822e+008 3.0432e+007 0.0000e+000 0.0000e+000 0.0000e+000 2.2659e+004 * * 3.0432e+007 5.1278e+007 0.0000e+000 0.0000e+000 0.0000e+000 1.2101e+004 * * 0.0000e+000 0.0000e+000 3.3674e+007 2.2659e+004 1.2101e+004 0.0000e+000 * * 0.0000e+000 0.0000e+000 2.2659e+004 2.9393e+002 7.7333e+001 0.0000e+000 * * 0.0000e+000 0.0000e+000 1.2101e+004 7.7333e+001 1.1465e+002 0.0000e+000 * * 2.2659e+004 1.2101e+004 0.0000e+000 0.0000e+000 0.0000e+000 8.4087e+001 * * * * END OF PROGRAM. * * ----------------------------------------------------------------------------------- * * REFERENCE: "MATERIAUX COMPOSITES - Comportement mécanique et analyse des structures * * By J.-M. Berthelot - MASSON 1996". * * * * C++ Release By J-P MOREAU - August 1997. * * (www.jpmoreau.fr) * **************************************************************************************/ #include #include // MAXCOU = maximum number of layers #define NMAX 3 #define MAXCOU 10 typedef double MAT[NMAX+1][NMAX+1]; typedef double MAT1[MAXCOU+1][NMAX+1][NMAX+1]; typedef double MAT6[7][7]; double EL,ET,HM,NULT,GLT,PI; // TH, EP = layer angle and thickness double TH[MAXCOU+1], EP[MAXCOU+1]; double H[MAXCOU+1]; MAT Q0, A, B, D; MAT1 Q; MAT6 R; int I,J,K, NCOUCHES; // set a matrix (N,N) to zero void MATNUL(int N, MAT A) { int I,J; for (I=1; I<=N; I++) for (J=1; J<=N; J++) A[I][J]=0.0; } // A(3,3) = B(3,3) void MATEGAL(int K, MAT1 A, MAT B) { int I,J; for (I=1; I<=NMAX; I++) for (J=1; J<=NMAX; J++) A[K][I][J]=B[I][J]; } // multiply a 3x3 matrixpar by a 3x1 vector void MATMULT(int N, MAT A, double *B, double *C) { double SUM; int I,K; for (I=1; I<=N; I++) { SUM= 0.0; for (K=1; K<=N; K++) SUM += A[I][K]*B[K]; C[I]=SUM; } } // print a NxN matrix void MATPR2(char *title, int N, MAT A) { int I,J; printf("\n%s\n", title); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) printf(" %10.4f", A[I][J]); printf("\n"); } } void MATPR1(char *title, int K, int N, MAT1 A) { int I,J; printf("\n%s\n", title); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) printf(" %10.4f", A[K][I][J]); printf("\n"); } } void MATPR6(char *title, int N, MAT6 A) { int I,J; printf("\n%s\n", title); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) printf(" %12.4e", A[I][J]); printf("\n"); } } double dcos4(double t) { double a; a=cos(t); return (a*a*a*a); } double dsin4(double t) { double a; a=sin(t); return (a*a*a*a); } // calculate the reduced stiffness matrix for theta=0 void CAL_Q0(MAT Q) { Q[1][1] = EL/(1.0-((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; } // calculate the reduced stiffness Q matrix for theta <> 0 // with the condition that the Q matrix for theta=0 is available // called here Q2 void CAL_Q(int K, MAT Q2, MAT1 Q, double TH) { double TEMP; TEMP = sin(TH)*sin(TH)*cos(TH)*cos(TH); Q[K][1][1]= Q2[1][1]*dcos4(TH)+Q2[2][2]*dsin4(TH)+2.0*(Q2[1][2]+2.0*Q2[3][3])*TEMP; Q[K][1][2]= (Q2[1][1]+Q2[2][2]-4.0*Q2[3][3])*TEMP+Q2[1][2]*(dcos4(TH)+dsin4(TH)); Q[K][2][1]=Q[K][1][2]; TEMP = sin(TH)*cos(TH)*cos(TH)*cos(TH); Q[K][1][3]=(Q2[1][1]-Q2[1][2]-2.0*Q2[3][3])*TEMP; TEMP = sin(TH)*sin(TH)*sin(TH)*cos(TH); Q[K][1][3]=Q[K][1][3]+(Q2[1][2]-Q2[2][2]+2.0*Q2[3][3])*TEMP; Q[K][3][1]=Q[K][1][3]; TEMP = sin(TH)*sin(TH)*cos(TH)*cos(TH); Q[K][2][2]= Q2[1][1]*dsin4(TH)+2.0*(Q2[1][2]+2.0*Q2[3][3])*TEMP+Q2[2][2]*dcos4(TH); TEMP = sin(TH)*sin(TH)*sin(TH)*cos(TH); Q[K][2][3]=(Q2[1][1]-Q2[1][2]-2.0*Q2[3][3])*TEMP; TEMP = cos(TH)*cos(TH)*cos(TH)*sin(TH); Q[K][2][3]=Q[K][2][3]+(Q2[1][2]-Q2[2][2]+2.0*Q2[3][3])*TEMP; Q[K][3][2]=Q[K][2][3]; TEMP = sin(TH)*sin(TH)*cos(TH)*cos(TH); Q[K][3][3]= (Q2[1][1]+Q2[2][2]-2.0*(Q2[1][2]+Q2[3][3]))*TEMP+Q2[3][3]*(dsin4(TH)+dcos4(TH)); } void main() { PI = 4.0*atan(1.0); // material constants in GPa EL = 38.0; ET = 9.0; NULT = 0.32; GLT = 3.6; // material data printf("\n CALCULATE THE STIFFNESS MATRIX OF A LAMINATE\n\n"); // number of layers NCOUCHES=4; // define angles and thicknesses of layers 1-4 TH[1]=-30.0; EP[1]=1e-3; TH[2]= 15.0; EP[2]=1.5e-3; TH[3]=-15.0; EP[3]=1.5e-3; TH[4]= 30.0; EP[4]=1e-3; printf(" Number of layers: %d\n\n", NCOUCHES); for (I=1; I<=NCOUCHES; I++) { printf(" Layer #%d: angle %6.1f deg. - thickness %6.1f mm\n", I, TH[I], 1000*EP[I]); //put angles in radians TH[I] = TH[I]*PI/180.0; } printf("\n\n Material Parameters:\n\n"); printf(" El = %6.1f Et = %6.1f\n", EL, ET); printf(" NUlt = %6.2f Glt= %6.1f\n\n", NULT, GLT); // reduced stiffness matrix (angle = 0) MATNUL(NMAX,Q0); CAL_Q0(Q0); MATPR2 (" BASIC REDUCED STIFFNESS MATRIX IN GPa:", NMAX, Q0); // calculate reduced stiffness matrices for different angles <> 0 for (I=1; I<=NCOUCHES; I++) { if (TH[I] == 0.0) MATEGAL(I,Q,Q0); else CAL_Q(I,Q0,Q,TH[I]); printf("\n Layer #%d\n", I); MATPR1(" REDUCED STIFFNESS MATRIX IN GPa:", I, NMAX, Q); } // calculate A matrix // A = sum of 1e3*EP(k)*Q(K) MATNUL(NMAX,A); for (I=1; I<=NMAX; I++) for (J=1; J<=NMAX; J++) { for (K=1; K<=NCOUCHES; K++) A[I][J]=A[I][J]+1000*(EP[K]*Q[K][I][J]); if (A[I][J] < 1e-8) A[I][J]=0.0; } MATPR2(" A MATRIX IN 1E6 N/M:", NMAX, A); // calculate B matrix // B = sum of 1e6*(1/2)*(h(k)*h(k)-h(k-1)*h(k-1))*Q(k) MATNUL(NMAX,B); HM=0.0; for (I=1; I<=NCOUCHES; I++) HM += EP[I]; // average height of stack HM=HM/2.0; // lower height layer #1 H[1] = -HM; for (I=2; I<=NCOUCHES+1; I++) H[I]=H[I-1]+EP[I-1]; for (I=1; I<=NMAX; I++) for (J=1; J<=NMAX; J++) { for (K=1; K<=NCOUCHES; K++) B[I][J]=B[I][J]+1e6*0.5*((H[K+1]*H[K+1]-H[K]*H[K])*Q[K][I][J]); if (B[I][J] < 1e-8) B[I][J]=0.0; } MATPR2(" B MATRIX IN 1E3 N:", NMAX, B); // 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) MATNUL(NMAX,D); for (I=1; I<=NMAX; I++) for (J=1; J<=NMAX; J++) { for (K=2; K<=NCOUCHES+1; K++) D[I][J]=D[I][J]+1e9*(1.0/3.0)*((H[K]*H[K]*H[K]-H[K-1]*H[K-1]*H[K-1])*Q[K-1][I][J]); if (D[I][J] < 1e-8) D[I][J]=0.0; } MATPR2(" D MATRIX IN N.DM:", NMAX, D); /* assemble the stiffness matrix of the laminate ! A B ! R = ! ! ! B D ! */ for (I=1; I<7; I++) for (J=1; J<7; J++) R[I][J]=0.0; for (I=1; I<=NMAX; I++) for (J=1; J<=NMAX; J++) R[I][J] = 1e6*A[I][J]; for (I=1; I<=NMAX; I++) for (J=4; J<=NMAX+3; J++) R[I][J] = 1e3*B[I][J-3]; for (I=4; I<=NMAX+3; I++) for (J=1; J<=NMAX; J++) R[I][J] = 1e3*B[I-3][J]; for (I=4; I<=NMAX+3; I++) for (J=4; J<=NMAX+3; J++) R[I][J] = D[I-3][J-3]; MATPR6(" STIFFNESS R MATRIX:", 6, R); printf("\n END OF PROGRAM.\n\n"); } // end of file compos03.cpp