clt.c (31131B)
1 /* 2 * CLASSICAL LAMINATE THEORY 3 * clt.h 4 */ 5 6 #include <stdio.h> 7 #include <stdlib.h> 8 #include <string.h> 9 #include "utils.c" 10 #include "rw.c" 11 #include "clt.h" 12 13 14 void 15 assemble_Z(double *Z, Laminate *lam, bool debugZ) 16 { 17 double total_thk = 0.0; 18 double h; 19 for(i = 0; i < lam->num_layers; i++) h += lam->thk[i]; 20 21 for (i = 0; i < lam->num_layers; ++i) { 22 total_thk += lam->thk[i]; 23 } 24 25 Z[0] = -total_thk / 2.0; 26 for (i = 1; i <= lam->num_layers; ++i) { 27 Z[i] = Z[i-1] + lam->thk[i-1]; 28 } 29 30 if (debugZ) { 31 printf("h = %0.4e mm\n", h); 32 printf("Layer Interface Positions (Z):\n"); 33 for (i = 0; i <= lam->num_layers; ++i) { 34 printf("Z[%d] = %.4e mm\n", i, Z[i]); 35 } 36 } 37 } 38 39 void 40 calc_thermal_forces(double *Nt, MaterialProperties *mat_list, 41 Laminate *lam, double *Z, int *fail_list, double dT) 42 { 43 double a[3], a_LCS[3]; 44 double Q[3][3], Qbar[3][3], T[3][3], Ti[3][3], temp[3][3]; 45 46 for (i = 0; i < lam->num_layers; ++i) { 47 MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]]; 48 a[0] = mat_prop->a1; a[1] = mat_prop->a2; a[2] = 0.0; 49 50 assemble_T(T, lam->ang[i]); 51 52 /* Transform to Laminate Coordinate System */ 53 for (j = 0; j < 3; ++j) { 54 for (k = 0; k < 3; ++k) { 55 a_LCS[j] += T[k][j] * a[k]; 56 } 57 } 58 59 /* Assemble Q matrix */ 60 assemble_Q(Q, mat_prop); 61 62 /* Calculate Qbar */ 63 transpose(T,Ti); 64 matmul(Ti,Q,temp); 65 matmul(Q,T,Qbar); 66 67 /* Calculate force resultants */ 68 for (j = 0; j < 3; ++j) { 69 for (k = 0; k < 3; ++k) { 70 double force = Qbar[j][k] * a_LCS[k] * dT; 71 Nt[j] += force * (Z[i+1] - Z[i]); 72 Nt[j+3] += force * 0.5 * 73 (Z[i+1]*Z[i+1] - Z[i]*Z[i]); 74 } 75 } 76 } 77 } 78 79 void 80 calc_moisture_forces(double *Nm, MaterialProperties *mat_list, 81 Laminate *lam, double *Z, int *fail_list, double dM) 82 { 83 double b[3], b_LCS[3]; 84 double Q[3][3], Qbar[3][3], T[3][3], Ti[3][3], temp[3][3]; 85 86 for (i = 0; i < lam->num_layers; ++i) { 87 MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]]; 88 b[0] = mat_prop->b1; b[1] = mat_prop->b2; b[2] = 0.0; 89 90 /* Apply transformed matrix */ 91 assemble_T(T, lam->ang[i]); 92 93 /* Transform to Laminate Coordinate System */ 94 for (j = 0; j < 3; ++j) { 95 for (k = 0; k < 3; ++k) { 96 b_LCS[j] += T[k][j] * b[k]; 97 } 98 } 99 100 /* Assemble Q matrix */ 101 assemble_Q(Q, mat_prop); 102 103 /* Calculate Qbar */ 104 transpose(T,Ti); 105 matmul(Ti,Q,temp); 106 matmul(Q,T,Qbar); 107 108 /* Calculate force resultants */ 109 for (j = 0; j < 3; ++j) { 110 for (k = 0; k < 3; ++k) { 111 double force = Qbar[j][k] * b_LCS[k] * dM; 112 Nm[j] += force * (Z[i+1] - Z[i]); 113 Nm[j+3] += force * 0.5 * 114 (Z[i+1]*Z[i+1] - Z[i]*Z[i]); 115 } 116 } 117 } 118 } 119 120 void 121 assemble_Q(double Q[3][3], MaterialProperties *mat_prop) 122 { 123 124 if (mat_prop->isotropic) { 125 double E, nu, E_div; 126 E = mat_prop->E1; 127 nu = mat_prop->nu12; 128 129 /*if (fail_type == 1 || fail_type == 2 || fail_type == 3) {*/ 130 /* double df = 0.001; */ 131 /* E *= df;*/ 132 /* nu *= df;*/ 133 /*}*/ 134 /**/ 135 E_div = E / (1 - nu * nu); 136 137 Q[0][0] = Q[1][1] = E_div; 138 Q[0][1] = Q[1][0] = nu * E_div; 139 Q[2][2] = E_div * (1 - nu) / 2.0; // Since G = E / (2(1+nu)) for isotropic 140 Q[0][2] = Q[1][2] = Q[2][0] = Q[2][1] = 0.0; 141 } else { 142 double E1 = mat_prop->E1, 143 E2 = mat_prop->E2, 144 nu12 = mat_prop->nu12, 145 G12 = mat_prop->G12; 146 double df = 0.001; /* Degradation factor */ 147 148 /* fiber or shear failure */ 149 /*if (fail_type == 1 || fail_type == 2) { */ 150 /* E1 *= df; */ 151 /* E2 *= df; */ 152 /* nu12 *= df; */ 153 /* G12 *= df;*/ 154 /*} else if (fail_type == 3) { */ 155 /* matrix failure */ 156 /* E2 *= df; */ 157 /* nu12 *= df; */ 158 /* G12 *= df;*/ 159 /*}*/ 160 161 double n21 = nu12 * E2 / E1; 162 163 Q[0][0] = E1 / (1 - nu12 * n21); 164 Q[0][1] = Q[1][0] = nu12 * E1 * E2 / (E1 - nu12 * nu12 * E2); 165 Q[1][1] = E1 * E2 / (E1 - nu12 * nu12 * E2); 166 Q[2][2] = G12; 167 Q[0][2] = Q[1][2] = Q[2][0] = Q[2][1] = 0.0; 168 } 169 } 170 171 void 172 assemble_T(double T[3][3], double angle) 173 { 174 double rad = angle * DEGTORAD; 175 176 double cos = Cos(rad), sin = Sin(rad); 177 double cs = cos * sin, cc = cos * cos, ss = sin * sin; 178 179 T[0][0] = cc; T[0][1] = ss; T[0][2] = cs; 180 T[1][0] = ss; T[1][1] = cc; T[1][2] = -cs; 181 T[2][0] = -2*cs; T[2][1] = 2*cs; T[2][2] = cc - ss; 182 } 183 184 void 185 assemble_ABD(double ABD[6][6], MaterialProperties *mat_list, 186 Laminate *lam, double *Z, bool debugABD) 187 { 188 double A[3][3] = {0}, B[3][3] = {0}, D[3][3] = {0}; 189 double Q[3][3], T[3][3], Ti[3][3]; 190 double temp[3][3], Qbar[3][3]; 191 192 for (i = 0; i < lam->num_layers; ++i) { 193 MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]]; 194 195 assemble_Q(Q, mat_prop); 196 assemble_T(T, lam->ang[i]); 197 198 /* Transpose T to get Ti */ 199 transpose(T,Ti); 200 201 /* Qbar = T^T * Q * T */ 202 matmul(Ti, Q, temp); 203 matmul(temp, T, Qbar); 204 205 /* B and D matrices */ 206 //double h = Z[i+1] - Z[i]; /* thickness of the layer */ 207 for (j = 0; j < 3; ++j) { 208 for (k = 0; k < 3; ++k) { 209 A[j][k] += Qbar[j][k] * (Z[i+1] - Z[i]); 210 B[j][k] += 0.5 * Qbar[j][k] * (Z[i+1] * Z[i+1] - Z[i] * Z[i]); 211 D[j][k] += CONST13 * Qbar[j][k] * (Z[i+1] * Z[i+1] * Z[i+1] - Z[i] * Z[i] * Z[i]); 212 } 213 } 214 } 215 216 /* Combine into ABD matrix */ 217 for (i = 0; i < 3; i++) { 218 for (j = 0; j < 3; j++) { 219 ABD[i][j] = A[i][j]; /* A matrix */ 220 ABD[i][j+3] = ABD[i+3][j] = B[i][j]; /* B matrix */ 221 ABD[i+3][j+3] = D[i][j]; /* D matrix */ 222 } 223 } 224 if (debugABD) { 225 dispmat(&ABD[0][0], 6, 6, 0); 226 } 227 } 228 229 /* Tsai-Wu criterion */ 230 void 231 fs_tsaiwu_2D(MaterialProperties *mat_prop, 232 double sig1, double sig2, double tau, 233 FSResult *result) 234 { 235 double f11 = 1 / (mat_prop->Xt * mat_prop->Xc); 236 double f22 = 1 / (mat_prop->Yt * mat_prop->Yc); 237 double f12 = -1 / (2 * sqrt(mat_prop->Xt * mat_prop->Xc * mat_prop->Yt * mat_prop->Yc)); 238 double f66 = 1 / (mat_prop->S12 * mat_prop->S12); 239 double f1 = 1/mat_prop->Xt - 1/mat_prop->Xc; 240 double f2 = 1/mat_prop->Yt - 1/mat_prop->Yc; 241 242 double a = f11*sig1*sig1 + f22*sig2*sig2 + f66*tau*tau + 2*f12*sig1*sig2; 243 double b = f1*sig1 + f2*sig2; 244 245 result->fs = (-b + sqrt(b*b + 4*a)) / (2*a); 246 247 double H1 = fabs(f1*sig1 + f11*sig1*sig1); 248 double H2 = fabs(f2*sig2 + f22*sig2*sig2); 249 double H6 = fabs(f66*tau*tau); 250 251 if (H1 >= H2 && H1 >= H6) { 252 strcpy(result->mode, "fiber"); 253 } else if (H2 >= H1 && H2 >= H6) { 254 strcpy(result->mode, "matrix"); 255 } else { 256 strcpy(result->mode, "shear"); 257 } 258 } 259 260 void 261 tsaiwu_2D(MaterialProperties *mat_list, Laminate *lam, double **stress_inf, 262 double **stress_sup, FSResult *fs_inf, FSResult *fs_sup) 263 { 264 for (i = 0; i < lam->num_layers; i++) { 265 MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]]; 266 267 /* Superior Stresses */ 268 fs_tsaiwu_2D(mat_prop, stress_sup[0][i], stress_sup[1][i], stress_sup[2][i], &fs_sup[i]); 269 270 /* Inferior Stresses */ 271 fs_tsaiwu_2D(mat_prop, stress_inf[0][i], stress_inf[1][i], stress_inf[2][i], &fs_inf[i]); 272 } 273 } 274 275 /* Maximum Stress criterion */ 276 void 277 fs_maxstress_2D(MaterialProperties *mat_prop, double sig1, double sig2, 278 double tau, FSResult *result) 279 { 280 double f_1 = sig1 > 0 ? sig1 / mat_prop->Xt : -sig1 / mat_prop->Xc; 281 double f_2 = sig2 > 0 ? sig2 / mat_prop->Yt : -sig2 / mat_prop->Yc; 282 double f_s = fabs(tau) / mat_prop->S12; 283 284 double f_max = fmax(fmax(f_1, f_2), f_s); 285 286 result->fs = 1 / f_max; 287 288 if (f_max == f_1) { 289 strcpy(result->mode, "fiber"); 290 } else if (f_max == f_2) { 291 strcpy(result->mode, "matrix"); 292 } else { 293 strcpy(result->mode, "shear"); 294 } 295 } 296 297 void 298 maxstress_2D(MaterialProperties *mat_list, Laminate *lam, 299 double **stress_inf, double **stress_sup, 300 FSResult *fs_inf, FSResult *fs_sup) 301 { 302 for (i = 0; i < lam->num_layers; i++) { 303 MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]]; 304 305 /* Superior Stresses */ 306 fs_maxstress_2D(mat_prop, stress_sup[0][i], stress_sup[1][i], 307 stress_sup[2][i], &fs_sup[i]); 308 309 /* Inferior Stresses */ 310 fs_maxstress_2D(mat_prop, stress_inf[0][i], stress_inf[1][i], 311 stress_inf[2][i], &fs_inf[i]); 312 } 313 } 314 315 /* Maximum Strain criterion is not directly translated due to lack of strain input in the given code. */ 316 317 /* Hashin criterion */ 318 void 319 fs_hashin_2D(MaterialProperties *mat_prop, 320 double sig1, double sig2, double tau, 321 FSResult *result) 322 { 323 double f_1 = sig1 >= 0 ? sig1 / mat_prop->Xt : -sig1 / mat_prop->Xc; 324 double f_2 = sig2 >= 0 ? sqrt((sig2/mat_prop->Yt)*(sig2/mat_prop->Yt) 325 + (tau/mat_prop->S12)*(tau/mat_prop->S12)) : 326 sqrt((sig2/mat_prop->Yc)*(sig2/mat_prop->Yc) + 327 (tau/mat_prop->S12)*(tau/mat_prop->S12)); 328 329 double f_max = fmax(f_1, f_2); 330 331 result->fs = 1 / f_max; 332 333 if (f_max == f_1) { 334 strcpy(result->mode, "fiber"); 335 } else { 336 strcpy(result->mode, "matrix"); 337 } 338 } 339 340 void 341 hashin_2D(MaterialProperties *mat_list, Laminate *lam, 342 double **stress_inf, double **stress_sup, 343 FSResult *fs_inf, FSResult *fs_sup) 344 { 345 for (i = 0; i < lam->num_layers; i++) { 346 MaterialProperties *mat_prop = &mat_list[lam->mat_id[i]]; 347 348 /* Superior Stresses */ 349 fs_hashin_2D(mat_prop, stress_sup[0][i], stress_sup[1][i], stress_sup[2][i], &fs_sup[i]); 350 351 /* Inferior Stresses */ 352 fs_hashin_2D(mat_prop, stress_inf[0][i], stress_inf[1][i], stress_inf[2][i], &fs_inf[i]); 353 } 354 } 355 356 double E_x(double ABD[6][6], double h) { 357 double matrix1[6][6], matrix2[5][5]; 358 double det1, det2, Ex; 359 360 // Copy relevant parts of ABD to matrix1 for the first determinant 361 for (i = 0; i < 6; i++) { 362 for (j = 0; j < 6; j++) { 363 matrix1[i][j] = ABD[i][j]; 364 } 365 } 366 367 // Perform LU decomposition for matrix1 368 if (LUDecompose(&matrix1[0][0], 6, 6) < 0) { 369 fprintf(stderr, "Error during LU decomposition for matrix1.\n"); 370 return -1; 371 } 372 373 // Calculate determinant for matrix1 374 det1 = LUDeterminant(&matrix1[0][0], 6, 6); 375 376 matrix2[0][0] = ABD[1][1]; 377 matrix2[0][1] = ABD[1][2]; 378 matrix2[0][2] = ABD[1][3]; 379 matrix2[0][3] = ABD[1][4]; 380 matrix2[0][4] = ABD[1][5]; 381 matrix2[1][0] = ABD[2][1]; 382 matrix2[1][1] = ABD[2][2]; 383 matrix2[1][2] = ABD[2][3]; 384 matrix2[1][3] = ABD[2][4]; 385 matrix2[1][4] = ABD[2][5]; 386 matrix2[2][0] = ABD[3][1]; 387 matrix2[2][1] = ABD[3][2]; 388 matrix2[2][2] = ABD[3][3]; 389 matrix2[2][3] = ABD[3][4]; 390 matrix2[2][4] = ABD[3][5]; 391 matrix2[3][0] = ABD[4][1]; 392 matrix2[3][1] = ABD[4][2]; 393 matrix2[3][2] = ABD[4][3]; 394 matrix2[3][3] = ABD[4][4]; 395 matrix2[3][4] = ABD[4][5]; 396 matrix2[4][0] = ABD[5][1]; 397 matrix2[4][1] = ABD[5][2]; 398 matrix2[4][2] = ABD[5][3]; 399 matrix2[4][3] = ABD[5][4]; 400 matrix2[4][4] = ABD[5][5]; 401 402 // Perform LU decomposition for matrix2 403 if (LUDecompose(&matrix2[0][0], 5, 5) < 0) { 404 fprintf(stderr, "Error during LU decomposition for matrix2.\n"); 405 return -1; 406 } 407 408 // Calculate determinant for matrix2 409 det2 = LUDeterminant(&matrix2[0][0], 5, 5); 410 411 // Check for division by zero 412 if (det2 == 0) { 413 fprintf(stderr, "Error: Division by zero in E_x calculation, determinant of the submatrix is zero.\n"); 414 return -1; 415 } 416 417 // Compute E_x 418 Ex = (det1 / det2) / h; 419 return Ex; 420 } 421 422 void 423 E_y(double ABD[6][6], double h, double Ey) 424 { 425 double matrix1[6][6], matrix2[5][5]; 426 double det1, det2; 427 428 // Copy relevant parts of ABD to matrix1 for the first determinant 429 for (i = 0; i < 6; i++) { 430 for (j = 0; j < 6; j++) { 431 matrix1[i][j] = ABD[i][j]; 432 } 433 } 434 435 // Perform LU decomposition for matrix1 436 if (LUDecompose(&matrix1[0][0], 6, 6) < 0) { 437 fprintf(stderr, "Error during LU decomposition for matrix1.\n"); 438 return; 439 } 440 441 // Calculate determinant for matrix1 442 det1 = LUDeterminant(&matrix1[0][0], 6, 6); 443 444 matrix2[0][0] = ABD[0][0]; 445 matrix2[0][1] = ABD[0][2]; 446 matrix2[0][2] = ABD[0][3]; 447 matrix2[0][3] = ABD[0][4]; 448 matrix2[0][4] = ABD[0][5]; 449 matrix2[1][0] = ABD[2][0]; 450 matrix2[1][1] = ABD[2][2]; 451 matrix2[1][2] = ABD[2][3]; 452 matrix2[1][3] = ABD[2][4]; 453 matrix2[1][4] = ABD[2][5]; 454 matrix2[2][0] = ABD[3][0]; 455 matrix2[2][1] = ABD[3][2]; 456 matrix2[2][2] = ABD[3][3]; 457 matrix2[2][3] = ABD[3][4]; 458 matrix2[2][4] = ABD[3][5]; 459 matrix2[3][0] = ABD[4][0]; 460 matrix2[3][1] = ABD[4][2]; 461 matrix2[3][2] = ABD[4][3]; 462 matrix2[3][3] = ABD[4][4]; 463 matrix2[3][4] = ABD[4][5]; 464 matrix2[4][0] = ABD[5][0]; 465 matrix2[4][1] = ABD[5][2]; 466 matrix2[4][2] = ABD[5][3]; 467 matrix2[4][3] = ABD[5][4]; 468 matrix2[4][4] = ABD[5][5]; 469 470 // Perform LU decomposition for matrix2 471 if (LUDecompose(&matrix2[0][0], 5, 5) < 0) { 472 fprintf(stderr, "Error during LU decomposition for matrix2.\n"); 473 return; 474 } 475 476 // Calculate determinant for matrix2 477 det2 = LUDeterminant(&matrix2[0][0], 5, 5); 478 479 // Check for division by zero 480 if (det2 == 0) { 481 fprintf(stderr, "Error: Division by zero in E_x calculation, determinant of the submatrix is zero.\n"); 482 return; 483 } 484 485 // Compute E_x 486 Ey = (det1 / det2) / h; 487 } 488 489 double 490 G_xy(double ABD[6][6], double h) 491 { 492 double matrix1[6][6], matrix2[5][5]; 493 double det1, det2, Gxy; 494 495 // Copy relevant parts of ABD to matrix1 for the first determinant 496 for (i = 0; i < 6; i++) { 497 for (j = 0; j < 6; j++) { 498 matrix1[i][j] = ABD[i][j]; 499 } 500 } 501 502 // Perform LU decomposition for matrix1 503 if (LUDecompose(&matrix1[0][0], 6, 6) < 0) { 504 fprintf(stderr, "Error during LU decomposition for matrix1.\n"); 505 return -1; 506 } 507 508 // Calculate determinant for matrix1 509 det1 = LUDeterminant(&matrix1[0][0], 6, 6); 510 511 matrix2[0][0] = ABD[0][0]; 512 matrix2[0][1] = ABD[0][1]; 513 matrix2[0][2] = ABD[0][3]; 514 matrix2[0][3] = ABD[0][4]; 515 matrix2[0][4] = ABD[0][5]; 516 517 matrix2[1][0] = ABD[1][0]; 518 matrix2[1][1] = ABD[1][1]; 519 matrix2[1][2] = ABD[1][3]; 520 matrix2[1][3] = ABD[1][4]; 521 matrix2[1][4] = ABD[1][5]; 522 523 matrix2[2][0] = ABD[3][0]; 524 matrix2[2][1] = ABD[3][1]; 525 matrix2[2][2] = ABD[3][3]; 526 matrix2[2][3] = ABD[3][4]; 527 matrix2[2][4] = ABD[3][5]; 528 529 matrix2[3][0] = ABD[4][0]; 530 matrix2[3][1] = ABD[4][1]; 531 matrix2[3][2] = ABD[4][3]; 532 matrix2[3][3] = ABD[4][4]; 533 matrix2[3][4] = ABD[4][5]; 534 535 matrix2[4][0] = ABD[5][0]; 536 matrix2[4][1] = ABD[5][1]; 537 matrix2[4][2] = ABD[5][3]; 538 matrix2[4][3] = ABD[5][4]; 539 matrix2[4][4] = ABD[5][5]; 540 541 // Perform LU decomposition for matrix2 542 if (LUDecompose(&matrix2[0][0], 5, 5) < 0) { 543 fprintf(stderr, "Error during LU decomposition for matrix2.\n"); 544 return -1; 545 } 546 547 // Calculate determinant for matrix2 548 det2 = LUDeterminant(&matrix2[0][0], 5, 5); 549 550 // Check for division by zero 551 if (det2 == 0) { 552 fprintf(stderr, "Error: Division by zero in E_x calculation, determinant of the submatrix is zero.\n"); 553 return -1; 554 } 555 556 // Compute E_x 557 Gxy = (det1 / det2) / h; 558 return Gxy; 559 } 560 561 double 562 nu_xy(double ABD[6][6], double h) 563 { 564 double matrix1[5][5], matrix2[5][5]; 565 double det1, det2, nxy; 566 567 matrix1[0][0] = ABD[0][1]; 568 matrix1[0][1] = ABD[1][2]; 569 matrix1[0][2] = ABD[1][3]; 570 matrix1[0][3] = ABD[1][4]; 571 matrix1[0][4] = ABD[1][5]; 572 matrix1[1][0] = ABD[0][2]; 573 matrix1[1][1] = ABD[2][2]; 574 matrix1[1][2] = ABD[2][3]; 575 matrix1[1][3] = ABD[2][4]; 576 matrix1[1][4] = ABD[2][5]; 577 matrix1[2][0] = ABD[0][3]; 578 matrix1[2][1] = ABD[3][2]; 579 matrix1[2][2] = ABD[3][3]; 580 matrix1[2][3] = ABD[3][4]; 581 matrix1[2][4] = ABD[3][5]; 582 matrix1[3][0] = ABD[0][4]; 583 matrix1[3][1] = ABD[4][2]; 584 matrix1[3][2] = ABD[4][3]; 585 matrix1[3][3] = ABD[4][4]; 586 matrix1[3][4] = ABD[4][5]; 587 matrix1[4][0] = ABD[0][5]; 588 matrix1[4][1] = ABD[5][2]; 589 matrix1[4][2] = ABD[5][3]; 590 matrix1[4][3] = ABD[5][4]; 591 matrix1[4][4] = ABD[5][5]; 592 593 /* Perform LU decomposition for matrix1 */ 594 if (LUDecompose(&matrix1[0][0], 5, 5) < 0) { 595 fprintf(stderr, "Error during LU decomposition for matrix1.\n"); 596 return -1; 597 } 598 599 /* Calculate determinant for matrix1 */ 600 det1 = LUDeterminant(&matrix1[0][0], 5, 5); 601 602 matrix2[0][0] = ABD[1][1]; 603 matrix2[0][1] = ABD[1][2]; 604 matrix2[0][2] = ABD[1][3]; 605 matrix2[0][3] = ABD[1][4]; 606 matrix2[0][4] = ABD[1][5]; 607 matrix2[1][0] = ABD[2][1]; 608 matrix2[1][1] = ABD[2][2]; 609 matrix2[1][2] = ABD[2][3]; 610 matrix2[1][3] = ABD[2][4]; 611 matrix2[1][4] = ABD[2][5]; 612 matrix2[2][0] = ABD[3][1]; 613 matrix2[2][1] = ABD[3][2]; 614 matrix2[2][2] = ABD[3][3]; 615 matrix2[2][3] = ABD[3][4]; 616 matrix2[2][4] = ABD[3][5]; 617 matrix2[3][0] = ABD[4][1]; 618 matrix2[3][1] = ABD[4][2]; 619 matrix2[3][2] = ABD[4][3]; 620 matrix2[3][3] = ABD[4][4]; 621 matrix2[3][4] = ABD[4][5]; 622 matrix2[4][0] = ABD[5][1]; 623 matrix2[4][1] = ABD[5][2]; 624 matrix2[4][2] = ABD[5][3]; 625 matrix2[4][3] = ABD[5][4]; 626 matrix2[4][4] = ABD[5][5]; 627 628 // Perform LU decomposition for matrix2 629 if (LUDecompose(&matrix2[0][0], 5, 5) < 0) { 630 fprintf(stderr, "Error during LU decomposition for submatrix.\n"); 631 return -1; 632 } 633 634 // Calculate determinant for matrix2 635 det2 = LUDeterminant(&matrix2[0][0], 5, 5); 636 637 // Check for division by zero 638 if (det2 == 0) { 639 fprintf(stderr, "Error: Division by zero in nu_xy calculation, determinant of the submatrix is zero.\n"); 640 return -1; 641 } 642 643 // Compute E_x 644 nxy = det1 / det2; 645 return nxy; 646 } 647 648 double 649 nu_yx(double ABD[6][6], double h) 650 { 651 double matrix1[5][5], matrix2[5][5]; 652 double det1, det2, nyx; 653 654 matrix1[0][0] = ABD[0][1]; 655 matrix1[0][1] = ABD[0][2]; 656 matrix1[0][2] = ABD[0][3]; 657 matrix1[0][3] = ABD[0][4]; 658 matrix1[0][4] = ABD[0][5]; 659 matrix1[1][0] = ABD[2][0]; 660 matrix1[1][1] = ABD[2][2]; 661 matrix1[1][2] = ABD[2][3]; 662 matrix1[1][3] = ABD[2][4]; 663 matrix1[1][4] = ABD[2][5]; 664 matrix1[2][0] = ABD[3][0]; 665 matrix1[2][1] = ABD[3][2]; 666 matrix1[2][2] = ABD[3][3]; 667 matrix1[2][3] = ABD[3][4]; 668 matrix1[2][4] = ABD[3][5]; 669 matrix1[3][0] = ABD[4][0]; 670 matrix1[3][1] = ABD[4][2]; 671 matrix1[3][2] = ABD[4][3]; 672 matrix1[3][3] = ABD[4][4]; 673 matrix1[3][4] = ABD[4][5]; 674 matrix1[4][0] = ABD[5][0]; 675 matrix1[4][1] = ABD[5][2]; 676 matrix1[4][2] = ABD[5][3]; 677 matrix1[4][3] = ABD[5][4]; 678 matrix1[4][4] = ABD[5][5]; 679 680 // Perform LU decomposition for matrix1 681 if (LUDecompose(&matrix1[0][0], 5, 5) < 0) { 682 fprintf(stderr, "Error during LU decomposition for matrix1.\n"); 683 return -1; 684 } 685 686 // Calculate determinant for matrix1 687 det1 = LUDeterminant(&matrix1[0][0], 5, 5); 688 689 matrix2[0][0] = ABD[0][0]; 690 matrix2[0][1] = ABD[0][2]; 691 matrix2[0][2] = ABD[0][3]; 692 matrix2[0][3] = ABD[0][4]; 693 matrix2[0][4] = ABD[0][5]; 694 matrix2[1][0] = ABD[2][0]; 695 matrix2[1][1] = ABD[2][2]; 696 matrix2[1][2] = ABD[2][3]; 697 matrix2[1][3] = ABD[2][4]; 698 matrix2[1][4] = ABD[2][5]; 699 matrix2[2][0] = ABD[3][0]; 700 matrix2[2][1] = ABD[3][2]; 701 matrix2[2][2] = ABD[3][3]; 702 matrix2[2][3] = ABD[3][4]; 703 matrix2[2][4] = ABD[3][5]; 704 matrix2[3][0] = ABD[4][0]; 705 matrix2[3][1] = ABD[4][2]; 706 matrix2[3][2] = ABD[4][3]; 707 matrix2[3][3] = ABD[4][4]; 708 matrix2[3][4] = ABD[4][5]; 709 matrix2[4][0] = ABD[5][0]; 710 matrix2[4][1] = ABD[5][2]; 711 matrix2[4][2] = ABD[5][3]; 712 matrix2[4][3] = ABD[5][4]; 713 matrix2[4][4] = ABD[5][5]; 714 715 // Perform LU decomposition for matrix2 716 if (LUDecompose(&matrix2[0][0], 5, 5) < 0) { 717 fprintf(stderr, "Error during LU decomposition for matrix2.\n"); 718 return -1; 719 } 720 721 // Calculate determinant for matrix2 722 det2 = LUDeterminant(&matrix2[0][0], 5, 5); 723 724 // Check for division by zero 725 if (det2 == 0) { 726 fprintf(stderr, "Error: Division by zero in n_yx calculation, determinant of the submatrix is zero.\n"); 727 return -1; 728 } 729 730 // Compute E_x 731 nyx = det1 / det2; 732 return nyx; 733 } 734 735 void 736 clt(Laminate *lam, MaterialProperties *materials, 737 int num_materials, double *F, 738 bool print, bool debugABD, bool debugZ) 739 { 740 /* Allocate and initialize Laminate System strain vectors */ 741 double LS_strain_sup[3][lam->num_layers]; 742 double LS_strain_inf[3][lam->num_layers]; 743 double MS_strain_sup[3][lam->num_layers]; 744 double MS_strain_inf[3][lam->num_layers]; 745 double MS_stress_sup[3][lam->num_layers]; 746 double MS_stress_inf[3][lam->num_layers]; 747 double T[3][3], Q[3][3]; 748 double *Z = malloc((lam->num_layers + 1) * sizeof(double)); 749 assemble_Z(Z, lam, debugZ); 750 751 double ABD[N][N]; 752 assemble_ABD(ABD, materials, lam, Z, debugABD); 753 754 double h = 0; 755 for(i = 0; i < lam->num_layers; i++) h += lam->thk[i]; 756 757 if (LUDecompose(&ABD[0][0], 6, 6) < 0) { 758 fprintf(stderr, "Error during LU decomposition for matrix1.\n"); 759 return; 760 } 761 LUSolve(&ABD[0][0], F, 6, 6); 762 763 double *strain = malloc(3 * sizeof(double)); 764 double *curvature = malloc(3 * sizeof(double)); 765 for (i = 0; i < 3; i++) strain[i] = F[i]; 766 for (i = 0; i < 3; i++) curvature[i] = F[i+3]; 767 768 for (i = 0; i < lam->num_layers; i++) { 769 for(j = 0; j < 3; j++) { 770 LS_strain_inf[j][i] = strain[j] + curvature[j] * Z[i]; 771 LS_strain_sup[j][i] = strain[j] + curvature[j] * Z[i + 1]; 772 } 773 } 774 775 for (i = 0; i < lam->num_layers; i++) { 776 MaterialProperties *mat_prop = &materials[lam->mat_id[i]]; 777 assemble_Q(Q, mat_prop); 778 assemble_T(T, lam->ang[i]); 779 double temp[3]; 780 for (j = 0; j < 3; j++) { 781 MS_strain_sup[j][i] = 0; 782 for (k = 0; k < 3; k++) 783 MS_strain_sup[j][i] += T[j][k] * LS_strain_sup[k][i]; 784 } 785 for (j = 0; j < 3; j++) { 786 MS_strain_inf[j][i] = 0; 787 for (k = 0; k < 3; k++) 788 MS_strain_inf[j][i] += T[j][k] * LS_strain_inf[k][i]; 789 } 790 for (j = 0; j < 3; j++) { 791 MS_stress_sup[j][i] = 0; 792 for (k = 0; k < 3; k++) 793 MS_stress_sup[j][i] += Q[j][k] * MS_strain_sup[k][i]; 794 } 795 for (j = 0; j < 3; j++) { 796 MS_stress_inf[j][i] = 0; 797 for (k = 0; k < 3; k++) 798 MS_stress_inf[j][i] += Q[j][k] * MS_strain_inf[k][i]; 799 } 800 } 801 802 if (print) { 803 print_equivalent_properties(ABD, h); 804 printf("\n"); 805 printf("Strain vector:\n"); 806 for (i = 0; i < 6; i++) { 807 printf("%e\n", F[i]); 808 } 809 print_lcs_mcs(lam->num_layers, 810 LS_strain_inf, LS_strain_sup, 811 MS_strain_sup, MS_strain_inf, 812 MS_stress_sup, MS_stress_inf); 813 } 814 815 free(Z); 816 free(strain); 817 free(curvature); 818 } 819 820 void 821 print_lcs_mcs(int num_layers, 822 double LS_strain_inf[3][MAX_LAYERS], 823 double LS_strain_sup[3][MAX_LAYERS], 824 double MS_strain_sup[3][MAX_LAYERS], 825 double MS_strain_inf[3][MAX_LAYERS], 826 double MS_stress_sup[3][MAX_LAYERS], 827 double MS_stress_inf[3][MAX_LAYERS]) 828 { 829 printf("\n"); 830 printf("+------------------------------------+\n"); 831 printf("| Laminate Coordinate System Strains |\n"); 832 printf("+------------------------------------+\n"); 833 printf("\n"); 834 printf("Layer | Inferior (Z-) | Superior (Z+)\n"); 835 printf("------|--------------------------------|----------------------------\n"); 836 837 for (j = 0; j < num_layers; j++) { 838 printf("%4d |", j + 1); 839 for (i = 0; i < 3; i++) { 840 printf("%10.2e", *(&LS_strain_inf[0][0] + i * num_layers + j)); 841 } 842 printf(" |"); 843 for (i = 0; i < 3; i++) { 844 printf("%10.2e", *(&LS_strain_sup[0][0] + i * num_layers + j)); 845 } 846 printf("\n"); 847 } 848 849 printf("\n"); 850 printf("+------------------------------------+\n"); 851 printf("| Material Coordinate System Strains |\n"); 852 printf("+------------------------------------+\n"); 853 printf("\n"); 854 printf("Layer | Inferior (Z-) | Superior (Z+)\n"); 855 printf("------|--------------------------------|----------------------------\n"); 856 857 for (j = 0; j < num_layers; j++) { 858 printf("%4d |", j + 1); 859 for (i = 0; i < 3; i++) { 860 printf("%10.2e", *(&MS_strain_inf[0][0] + i * num_layers + j)); 861 } 862 printf(" |"); 863 for (i = 0; i < 3; i++) { 864 printf("%10.2e", *(&MS_strain_sup[0][0] + i * num_layers + j)); 865 } 866 printf("\n"); 867 } 868 869 printf("\n"); 870 printf("+-----------------------------------+\n"); 871 printf("| Material Coordinate System Stress |\n"); 872 printf("+-----------------------------------+\n"); 873 printf("\n"); 874 printf("Layer | Inferior (Z-) | Superior (Z+)\n"); 875 printf("------|--------------------------------|----------------------------\n"); 876 for (j = 0; j < num_layers; j++) { 877 printf("%4d |", j + 1); 878 for (i = 0; i < 3; i++) { 879 printf("%10.2e", *(&MS_stress_inf[0][0] + i * num_layers + j)); 880 } 881 printf(" |"); 882 for (i = 0; i < 3; i++) { 883 printf("%10.2e", *(&MS_stress_sup[0][0] + i * num_layers + j)); 884 } 885 printf("\n"); 886 } 887 printf("\n"); 888 } 889 890 void 891 print_equivalent_properties(double ABD[6][6], double h) 892 { 893 double Ex = E_x(ABD, h); 894 double Ey = 0; E_y(ABD, h, Ey); 895 double Gxy = G_xy(ABD, h); 896 double nxy = nu_xy(ABD, h); 897 double nyx = nu_yx(ABD, h); 898 899 printf("\n"); 900 printf("+-----------------------+\n"); 901 printf("| Engineering Constants |\n"); 902 printf("+-----------------------+\n"); 903 printf("\n"); 904 905 printf("Equivalent elastic modulus Ex:\t%.4e\n", Ex); 906 printf("Equivalent elastic modulus Ey:\t%.4e\n", Ey); 907 printf("Equivalent shear modulus Gxy:\t%.4e\n", Gxy); 908 printf("Equivalent poisson ratio nu_xy:\t%.3f\n", nxy); 909 printf("Equivalent poisson ratio nu_yx:\t%.3f\n", nyx); 910 }