/************************************************************************************** * 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 = 1.0000e+06 NY = 5.0000e+05 NXY = 2.5000e+05 * * MX = 0.00000000 MY = 0.00000000 MXY = 0.00000000 * * * * 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: * * 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: * * 35.2120 4.6931 6.7316 * * 4.6931 9.4731 0.6986 * * 6.7316 0.6986 5.3416 * * * * Layer #:2 * * * * REDUCED STIFFNESS MATRIX IN GPa: * * 26.2896 8.1763 -9.4512 * * 8.1763 11.4292 -3.4183 * * -9.4512 -3.4183 8.8247 * * * * 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 * * --/-- * * STIFFNESS R MATRIX: * * 1.5822e+008 3.0432e+007 0.0000e+000 -1.3384e+004 5.2247e+003 -1.6154e+003 * * 3.0432e+007 5.1278e+007 0.0000e+000 5.2247e+003 2.9342e+003 5.9258e+003 * * 0.0000e+000 0.0000e+000 3.3674e+007 -1.6154e+003 5.9258e+003 5.2247e+003 * * -1.3384e+004 5.2247e+003 -1.6154e+003 3.2738e+002 6.4271e+001 6.0686e+001 * * 5.2247e+003 2.9342e+003 5.9258e+003 6.4271e+001 1.0732e+002 1.5438e+001 * * -1.6154e+003 5.9258e+003 5.2247e+003 6.0686e+001 1.5438e+001 7.1025e+001 * * * * DETERMINANT OF R: 4.2688e+029 * * * * INVERSE MATRIX OF R: * * 7.2071e+000 -4.3217e+000 6.9155e-002 4.1544e+002 -5.2541e+002 2.7865e+002 * * -4.3217e+000 2.2297e+001 2.7923e-001 -1.8738e+002 -4.2073e+001 -1.8099e+003 * * 6.9155e-002 2.7923e-001 3.0508e+001 1.0321e+003 -1.9209e+003 -2.7302e+003 * * 4.1544e+002 -1.8738e+002 1.0321e+003 4.0518e+006 -2.0576e+006 -3.0655e+006 * * -5.2541e+002 -4.2073e+001 -1.9209e+003 -2.0576e+006 1.0747e+007 -4.4512e+005 * * 2.7865e+002 -1.8099e+003 -2.7302e+003 -3.0655e+006 -4.4512e+005 1.7154e+007 * * --/-- * * Stresses in each layer in main axes: * * * * Layer #1 * * For z = -2.50000 mm: * * sigma L : 282.73562 MPa * * sigma T : 78.07995 MPa * * sigma LT: 45.15936 MPa * * For z = -1.00000 mm: * * sigma L : 288.58230 MPa * * sigma T : 70.53304 MPa * * sigma LT: 34.70107 MPa * * * * Layer #2 * * For z = -1.00000 mm: * * sigma L : 86.43422 MPa * * sigma T : 105.75930 MPa * * sigma LT: 5.73649 MPa * * For z = 0.00000 mm: * * sigma L : 111.92641 MPa * * sigma T : 96.96499 MPa * * sigma LT: 8.38896 MPa * * * * Layer #3 * * For z = 0.00000 mm: * * sigma L : 151.46586 MPa * * sigma T : 90.07486 MPa * * sigma LT: 21.12949 MPa * * For z = 1.50000 mm: * * sigma L : 192.64459 MPa * * sigma T : 76.37100 MPa * * sigma LT: 19.34600 MPa * * * * Layer #4 * * For z = 1.50000 mm: * * sigma L : 333.21160 MPa * * sigma T : 51.87584 MPa * * sigma LT: 8.77294 MPa * * For z = 2.50000 mm: * * sigma L : 317.90586 MPa * * sigma T : 50.19097 MPa * * sigma LT: 1.40860 MPa * * * * ----------------------------------------------------------------------------------- * * 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 #define MACH_EPS 2e-16 //index 0 not used except for inversion of R typedef double MAT[NMAX+1][NMAX+1]; typedef double MAT1[MAXCOU+1][NMAX+1][NMAX+1]; typedef double MAT2[MAXCOU+1][NMAX+1]; typedef double MAT6[7][7]; double EL,ET,HM,NULT,GLT,PI; // TH, EP = angle and thickness of layers double TH[MAXCOU+1], EP[MAXCOU+1], H[MAXCOU+1]; MAT A,B,D, Q0, T; MAT1 Q; MAT6 R, X; double N[7], E0[7], EDEB[7], EFIN[7]; double DET,NX,NY,NXY,MX,MY,MXY; MAT2 AK, AK1, BK, BK1; double TEMP[NMAX], TEMP1[NMAX]; int I,J,K, NCOUCHES; FILE *fp; // set a matrix (N,N) to zero void MATNUL(MAT A,int N) { int I,J; for (I=1; I<=N; I++) for (J=1; J<=N; J++) A[I][J]=0.0; } // A(K,3,3) = B(3,3) void MATEGAL(int K, int N, MAT1 A, MAT B) { int I,J; for (I=1; I<=N; I++) for (J=1; J<=N; J++) A[K][I][J]=B[I][J]; } // multiply a 3x3 matrix by a 3x1 vector void MATMUL(MAT A, double *B, double *C,int N) { 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; } } // multiply a 6x6 matrix by a 6x1 vector void MATMUL6(MAT6 A, double *B, double *C,int N) { 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; fprintf(fp,"\n%s\n", title); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) fprintf(fp," %10.4f", A[I][J]); fprintf(fp,"\n"); } } void MATPR1(char *title, int K, int N, MAT1 A) { int I,J; fprintf(fp,"\n%s\n", title); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) fprintf(fp," %10.4f", A[K][I][J]); fprintf(fp,"\n"); } } void MATPR6(char *title, int N, MAT6 A) { int I,J; fprintf(fp,"\n%s\n", title); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) fprintf(fp," %12.4e", A[I][J]); fprintf(fp,"\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)); } /****************************************** * SOLVING A LINEAR MATRIX SYSTEM AX = B * * with Gauss-Jordan method using full * * pivoting at each step. During the pro- * * cess, original AA and BB matrices are * * destroyed to spare storage location. * * --------------------------------------- * * INPUTS: AA MATRIX N*N * * BB MATRIX N*M * * --------------------------------------- * * OUTPUS: AA INVERSE OF AA N*N * * DET DETERMINANT OF AA * * BB SOLUTION MATRIX N*M * * --------------------------------------- * * NOTA - If M=0 inversion of AA matrix * * only. * ******************************************/ void MATINV(int N,int M,MAT6 AA,MAT6 BB,double *DET) { double PC[7],PL[7],CS[7]; double PV,PAV,temp,TT; int I,IK,J,JK,K; // Initializations : *DET= 1.0; for (I=0; I PAV) { PV=AA[I][J]; PAV= fabs(PV); IK=I; JK=J; } } // Search terminated, the pivot is in location I=IK, J=JK. // Memorizing pivot location: PC[K]=JK; PL[K]=IK; // DETERMINANT DET is actualised // If DET=0, ERROR MESSAGE and STOP // Machine dependant EPSMACH equals here 1e-20 if (IK!=K) *DET=-*DET; if (JK!=K) *DET=-*DET; *DET=*DET*PV; temp= fabs(*DET); if (temp < MACH_EPS) { fprintf(fp,"\n The determinant equals ZERO !!!\n"); return; } // POSITIONNING PIVOT IN K,K: if(IK!=K) for (I=0; I=0; I--) { IK=(int) PC[I]; if (IK==I) goto fin1; // EXCHANGE LINES I and PC(I) of matrix AA: for (J=0; J=0; J--) { JK=(int) PL[J]; if (JK==J) goto fin2; // EXCHANGE COLUMNS J ET PL(J) of matrix AA: for (I=0; I 0 for (I=1; I<=NCOUCHES; I++) { if (TH[I] == 0.0) MATEGAL(I,NMAX,Q,Q0); else CAL_Q(I,Q0,Q,TH[I]); fprintf(fp,"\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(A, NMAX); 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 (fabs(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(B,NMAX); 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 (fabs(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(D,NMAX); 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 (fabs(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); // start index in R from zero before using MATINV for (I=0; I<6; I++) for (J=0; J<6; J++) R[I][J] = R[I+1][J+1]; // Direct inversion of R matrix, stored in R MATINV(6,0,R,X,&DET); fprintf(fp,"\n DETERMINANT OF R: %12.4e\n\n", DET); // restart index in R from one to continue for (I=6; I>0; I--) for (J=6; J>0; J--) R[I][J] = R[I-1][J-1]; for (I=1; I<7; I++) for (J=1; J<7; J++) R[I][J]=1e9*R[I][J]; MATPR6(" INVERSE MATRIX OF R:", 6, R); // Calculate déformations eps0xx, eps0yy, gam0xy,kx,ky,kxy // in main axes, stored in matrix E0(6) MATMUL6(R,N,E0,6); fprintf(fp,"\n Membrane deformations (in mm) and curbatures:\n\n"); fprintf(fp," eps0xx: %12.5f\n", 1e3*E0[1]); fprintf(fp," eps0yy: %12.5f\n", 1e3*E0[2]); fprintf(fp," eps0xy: %12.5f\n", 1e3*E0[3]); fprintf(fp," kx : %12.5f\n", E0[4]); fprintf(fp," ky : %12.5f\n", E0[5]); fprintf(fp," kxy : %12.5f\n", E0[6]); /* Deformations in reference axes (x,y) ! epsxx ! ! eps0xx ! ! kx ! ! epsyy ! = ! eps0yy ! +z ! ky ! ! gamxy ! ! gam0xy ! ! kxy ! */ fprintf(fp,"\n Deformations in reference axes (x,y):\n"); fprintf(fp," For z = %12.5f mm:\n", 1e3*H[1]); for (I=1; I<=NMAX; I++) EDEB[I]=E0[I]+H[1]*E0[I+3]; fprintf(fp," epsxx: %12.5f mm\n", 1e3*EDEB[1]); fprintf(fp," epsyy: %12.5f mm\n", 1e3*EDEB[2]); fprintf(fp," epsxy: %12.5f mm\n", 1e3*EDEB[3]); fprintf(fp," For z = %12.5f mm:\n", 1e3*H[NCOUCHES+1]); for (I=1; I<=NMAX; I++) EFIN[I]=E0[I]+H[NCOUCHES+1]*E0[I+3]; fprintf(fp," epsxx: %12.5f mm\n", 1e3*EFIN[1]); fprintf(fp," epsyy: %12.5f mm\n", 1e3*EFIN[2]); fprintf(fp," epsxy: %12.5f mm\n", 1e3*EFIN[3]); // Deformations in each layer in main axes of layer fprintf(fp,"\n Deformations in each layer in main axes:\n"); for (K=1; K<=NCOUCHES; K++) { INITT(T,TH[K]); for (I=1; I<=NMAX; I++) TEMP[I]=E0[I]; MATMUL(T,TEMP,TEMP1,NMAX); for (I=1; I<=NMAX; I++) AK[K][I]=TEMP1[I]; for (I=1; I<=NMAX; I++) TEMP[I]=E0[I+3]; MATMUL(T,TEMP,TEMP1,NMAX); for (I=1; I<=NMAX; I++) BK[K][I]=TEMP1[I]; fprintf(fp,"\n Layer #%d\n", K); fprintf(fp," For z = %12.5f mm:\n", 1e3*H[K]); for (I=1; I<=NMAX; I++) EDEB[I]=AK[K][I]+H[K]*BK[K][I]; fprintf(fp," epsxx: %12.5f mm\n", 1e3*EDEB[1]); fprintf(fp," epsyy: %12.5f mm\n", 1e3*EDEB[2]); fprintf(fp," epsxy: %12.5f mm\n", 1e3*EDEB[3]); fprintf(fp," For z = %12.5f mm:\n", 1e3*H[K+1]); for (I=1; I<=NMAX; I++) EFIN[I]=AK[K][I]+H[K+1]*BK[K][I]; fprintf(fp," epsxx: %12.5f mm\n", 1e3*EFIN[1]); fprintf(fp," epsyy: %12.5f mm\n", 1e3*EFIN[2]); fprintf(fp," epsxy: %12.5f mm\n", 1e3*EFIN[3]); } // Stresses in each layer in reference axes fprintf(fp,"\n Stresses in each layer in reference axes:\n"); for (K=1; K<=NCOUCHES; K++) { for (I=1; I<=NMAX; I++) TEMP[I]=E0[I]; for (I=1; I<=NMAX; I++) for (J=1; J<=NMAX; J++) A[I][J]=Q[K][I][J]; MATMUL(A,TEMP,TEMP1,NMAX); for (I=1; I<=NMAX; I++) AK[K][I]=TEMP1[I]; for (I=1; I<=NMAX; I++) TEMP[I]=E0[I+3]; MATMUL(A,TEMP,TEMP1,NMAX); for (I=1; I<=NMAX; I++) BK[K][I]=TEMP1[I]; fprintf(fp,"\n Layer #%d\n", K); fprintf(fp," For z = %12.5f mm:\n", 1e3*H[K]); for (I=1; I<=NMAX; I++) EDEB[I]=1e3*(AK[K][I]+H[K]*BK[K][I]); fprintf(fp," sigma L : %12.5f MPa\n", EDEB[1]); fprintf(fp," sigma T : %12.5f MPa\n", EDEB[2]); fprintf(fp," sigma LT: %12.5f MPa\n", EDEB[3]); fprintf(fp," For z = %12.5f mm:\n", 1e3*H[K+1]); for (I=1; I<=NMAX; I++) EFIN[I]=1e3*(AK[K][I]+H[K+1]*BK[K][I]); fprintf(fp," sigma L : %12.5f MPa\n", EFIN[1]); fprintf(fp," sigma T : %12.5f MPa\n", EFIN[2]); fprintf(fp," sigma LT: %12.5f MPa\n", EFIN[3]); } // Stresses in each layer in main axes fprintf(fp,"\n Stresses in each layer in main axes:\n"); for (K=1; K<=NCOUCHES; K++) { INITT1(T,TH[K]); for (I=1; I<=NMAX; I++) TEMP[I]=AK[K][I]; MATMUL(T,TEMP,TEMP1,NMAX); for (I=1; I<=NMAX; I++) AK1[K][I]=TEMP1[I]; for (I=1; I<=NMAX; I++) TEMP[I]=BK[K][I]; MATMUL(T,TEMP,TEMP1,NMAX); for (I=1; I<=NMAX; I++) BK1[K][I]=TEMP1[I]; fprintf(fp,"\n Layer #%d\n", K); fprintf(fp," For z = %12.5f mm:\n", 1e3*H[K]); for (I=1; I<=NMAX; I++) EDEB[I]=1e3*(AK1[K][I]+H[K]*BK1[K][I]); fprintf(fp," sigma L : %12.5f MPa\n", EDEB[1]); fprintf(fp," sigma T : %12.5f MPa\n", EDEB[2]); fprintf(fp," sigma LT: %12.5f MPa\n", EDEB[3]); fprintf(fp," For z = %12.5f mm:\n", 1e3*H[K+1]); for (I=1; I<=NMAX; I++) EFIN[I]=1E3*(AK1[K][I]+H[K+1]*BK1[K][I]); fprintf(fp," sigma L : %12.5f MPa\n", EFIN[1]); fprintf(fp," sigma T : %12.5f MPa\n", EFIN[2]); fprintf(fp," sigma LT: %12.5f MPa\n", EFIN[3]); } fprintf(fp,"\n"); fclose(fp); printf("\n END OF PROGRAM.\n\n"); } // End of file Compos04.cpp