1 /* 2 * This program is free software; you can redistribute it and/or 3 * modify it under the terms of the GNU General Public License as 4 * published by the Free Software Foundation; either version 2 of the 5 * License, or (at your option) any later version. 6 * 7 * This program is distributed in the hope that it will be useful, but 8 * WITHOUT ANY WARRANTY; without even the implied warranty of 9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOUSE. See the GNU 10 * General Public License for more details. 11 * 12 * You should have recieved a copy of the GNU General Public License 13 * along with this program; if not write to the Free Software 14 * Foundation, inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 15 * USA 16 */ 17 18 package run.materials; 19 20 import run.*; 21 22 23 /** 24 * Insert the type's description here. Creation date: (09/09/01 %T) 25 * 26 * @author: Jonas Forssell 27 */ 28 public class Elastic extends Material implements Cloneable { 29 private Jama.Matrix new_stress; 30 private boolean nu_initialized = false; 31 private boolean youngs_modulus_initialized = false; 32 private Jama.Matrix stiffness_matrix_3d; 33 private Jama.Matrix stiffness_matrix_plane_stress; 34 private boolean E_is_set; 35 private boolean NU_is_set; 36 37 /** 38 * This material law calculates stresses based on strain, strain increase 39 * and old stress as input. Since this material law is isotropic and purly 40 * elastic, the history is not important. Therefore, the old stresses are 41 * never used, but are kept here for consistency. The constitutive 42 * matrices are defined in the setInitialConditions method. 43 */ Elastic()44 public Elastic() { 45 type = new String("ELASTIC"); 46 new_stress = new Jama.Matrix(6, 1); 47 stiffness_matrix_3d = new Jama.Matrix(6, 6); 48 stiffness_matrix_plane_stress = new Jama.Matrix(6, 6); 49 } 50 51 /** 52 * This method calculates a stress based on a given strain, strain increase 53 * and stress for an element This is the default input for a material law. 54 * In this case, the job is simple. Just to multiply with youngs modulus 55 * and then we are OK. Note that the incoming objects are modified. No 56 * return on the method is needed. All Matrices are assumed to be of the 57 * standard 6x1 type: strain = [epsilonx, epsilony, epsilonz, gammaxy, 58 * gammayz, gammazx]T dstrain = same format as strain but latest increment 59 * stress = [sigmax, sigmay, sigmaz, tauxy, tauyz, tauzx]T Creation date: 60 * (17/12/01 %T) 61 * 62 * @param strain double 63 */ calculateStressOneDimensional( Jama.Matrix strain, Jama.Matrix dstrain, Jama.Matrix stress, double timestep )64 public void calculateStressOneDimensional( 65 Jama.Matrix strain, Jama.Matrix dstrain, Jama.Matrix stress, 66 double timestep 67 ) 68 { 69 strain.plusEquals(dstrain); 70 71 // Update only one value in stress, since this is one dimensional 72 stress.set(0, 0, strain.get(0, 0) * youngs_modulus); 73 } 74 75 /** 76 * This method calculates a stress based on a given strain, strain increase 77 * and stress for an element This is the default input for a material law. 78 * In this case, the job is simple. Just to multiply with the constitutive 79 * matrix and then we are OK. Note that the incoming objects are modified. 80 * No return on the method is needed. All Matrices are assumed to be of 81 * the standard 6x1 type: strain = [epsilonx, epsilony, epsilonz, gammaxy, 82 * gammayz, gammazx]T dstrain = same format as strain but latest increment 83 * stress = [sigmax, sigmay, sigmaz, tauxy, tauyz, tauzx]T Creation date: 84 * (13/12/01 %T) 85 * 86 * @param stress Jama.Matrix 87 * 88 * @return Jama.Matrix 89 * 90 */ calculateStressThreeDimensional( Jama.Matrix strain, Jama.Matrix dstrain, Jama.Matrix stress, double timestep )91 public void calculateStressThreeDimensional( 92 Jama.Matrix strain, Jama.Matrix dstrain, Jama.Matrix stress, 93 double timestep 94 ) 95 { 96 strain.plusEquals(dstrain); 97 stress.setMatrix(0, 5, 0, 0, stiffness_matrix_3d.times(strain)); 98 } 99 100 /** 101 * This method calculates a stress based on a given strain, strain increase 102 * and stress for an element This is the default input for a material law. 103 * In this case, the job is simple. Just to multiply with the constitutive 104 * matrix and then we are OK. Note that the incoming objects are modified. 105 * No return on the method is needed. All Matrices are assumed to be of 106 * the standard 6x1 type: strain = [epsilonx, epsilony, epsilonz, gammaxy, 107 * gammayz, gammazx]T dstrain = same format as strain but latest increment 108 * stress = [sigmax, sigmay, sigmaz, tauxy, tauyz, tauzx]T Creation date: 109 * (13/12/01 %T) 110 * 111 * @param stress Jama.Matrix 112 * 113 * @return Jama.Matrix 114 * 115 */ calculateStressTwoDimensionalPlaneStress( Jama.Matrix strain, Jama.Matrix dstrain, Jama.Matrix stress, double timestep )116 public void calculateStressTwoDimensionalPlaneStress( 117 Jama.Matrix strain, Jama.Matrix dstrain, Jama.Matrix stress, 118 double timestep 119 ) 120 { 121 /* This is the simple and clean way of calculating the stress 122 stress.setMatrix(0,5,0,0,stiffness_matrix_plane_stress.times(strain)); 123 * 124 * but there are many zeroes there, and speed by can be improved 37% by just 125 * multiplying the terms that are nonzero, so here it is: 126 * stress = stress + dstress 127 */ 128 stress.set( 129 0, 0, 130 stress.get(0, 0) + 131 (stiffness_matrix_plane_stress.get(0, 0) * dstrain.get(0, 0)) + 132 (stiffness_matrix_plane_stress.get(0, 1) * dstrain.get(1, 0)) 133 ); 134 stress.set( 135 1, 0, 136 stress.get(1, 0) + 137 (stiffness_matrix_plane_stress.get(1, 0) * dstrain.get(0, 0)) + 138 (stiffness_matrix_plane_stress.get(1, 1) * dstrain.get(1, 0)) 139 ); 140 stress.set( 141 3, 0, 142 stress.get(3, 0) + 143 (stiffness_matrix_plane_stress.get(3, 3) * dstrain.get(3, 0)) 144 ); 145 stress.set( 146 4, 0, 147 stress.get(4, 0) + 148 (stiffness_matrix_plane_stress.get(4, 4) * dstrain.get(4, 0)) 149 ); 150 stress.set( 151 5, 0, 152 stress.get(5, 0) + 153 (stiffness_matrix_plane_stress.get(5, 5) * dstrain.get(5, 0)) 154 ); 155 156 // To predict thinning of element, update the strain increment in thickness direction as an isotropic 157 dstrain.set( 158 2, 0, -(nu / (1 - nu)) * (dstrain.get(0, 0) + dstrain.get(1, 0)) 159 ); 160 161 // Update the strain 162 strain.plusEquals(dstrain); 163 } 164 165 /** 166 * Insert the method's description here. Creation date: (13/09/01 %T) 167 * 168 * @param arg1 java.lang.String 169 * @param arg2 java.lang.String 170 * @param arg3 java.lang.String 171 * 172 * @exception java.lang.IllegalArgumentException The exception description. 173 */ parse_Fembic(Token[] arg, int lineno)174 public void parse_Fembic(Token[] arg, int lineno) 175 throws java.text.ParseException 176 { 177 // The Youngs modulus of the material is defined 178 int i = 0; 179 180 while (i < arg.length) { 181 if ( 182 arg[i].getw().toUpperCase().equals("E") && 183 arg[i + 1].getw().toUpperCase().equals("=") 184 ) { 185 // The value of the cross section area is in arg3. Set this in the element 186 youngs_modulus = arg[i + 2].getn(); 187 youngs_modulus_initialized = true; 188 i += 3; 189 E_is_set = true; 190 } else 191 // The density of the material is defined 192 if ( 193 arg[i].getw().toUpperCase().equals("RHO") && 194 arg[i + 1].getw().toUpperCase().equals("=") 195 ) { 196 // The value of the cross section area is in arg3. Set this in the element 197 density = arg[i + 2].getn(); 198 density_is_set = true; 199 i += 3; 200 } else 201 // Poissons constant of the material is defined 202 if ( 203 arg[i].getw().toUpperCase().equals("NU") && 204 arg[i + 1].getw().toUpperCase().equals("=") 205 ) { 206 // The value of the cross section area is in arg3. Set this in the element 207 nu = arg[i + 2].getn(); 208 nu_initialized = true; 209 i += 3; 210 NU_is_set = true; 211 } else 212 // Failure strain for the material is set 213 if ( 214 arg[i].getw().toUpperCase().equals("FAILURE_STRAIN") && 215 arg[i + 1].getw().toUpperCase().equals("=") 216 ) { 217 failure_strain = arg[i + 2].getn(); 218 i += 3; 219 failure_strain_is_set = true; 220 } else 221 // Failure strain for the material is set 222 if ( 223 arg[i].getw().toUpperCase().equals("FAILURE_STRESS") && 224 arg[i + 1].getw().toUpperCase().equals("=") 225 ) { 226 failure_stress = arg[i + 2].getn(); 227 i += 3; 228 failure_stress_is_set = true; 229 } else { 230 throw new java.text.ParseException( 231 "Unknown material parameter in line ", lineno 232 ); 233 } 234 } 235 } 236 237 /** 238 * Insert the method's description here. Creation date: (13/09/01 %T) 239 * 240 * @param arg1 java.lang.String 241 * @param arg2 java.lang.String 242 * @param arg3 java.lang.String 243 * 244 * @exception java.lang.IllegalArgumentException The exception description. 245 */ parse_Nastran(Token[] param, int lineno)246 public void parse_Nastran(Token[] param, int lineno) 247 throws java.text.ParseException 248 { 249 int i = 0; 250 int index; 251 Constraint temp_constraint; 252 Load temp_load; 253 254 if (param[1].is_a_number()) { 255 name = new String("Mat_" + (int)param[1].getn()); 256 name_is_set = true; 257 } 258 else 259 throw new java.text.ParseException( 260 "Illegal identification number of MAT1: " + param[1].getn(),lineno 261 ); 262 263 if (param[2].is_a_number()) { 264 youngs_modulus = param[1].getn(); 265 E_is_set = true; 266 youngs_modulus_initialized = true; 267 268 } else { 269 youngs_modulus = 2*(1+param[4].getn())*param[3].getn(); 270 E_is_set = true; 271 youngs_modulus_initialized = true; 272 273 } 274 275 if (param[4].is_a_number()) { 276 nu = param[4].getn(); 277 NU_is_set = true; 278 nu_initialized = true; 279 280 281 } else { 282 nu = youngs_modulus / (2 * param[3].getn()) -1; 283 NU_is_set = true; 284 nu_initialized = true; 285 } 286 287 if (param[5].is_a_number()) { 288 density = param[5].getn(); 289 density_is_set = true; 290 291 } 292 else 293 throw new java.text.ParseException( 294 "Illegal density for MAT1: " + param[1].getn(),lineno 295 ); 296 297 } 298 299 /** 300 * Insert the method's description here. Creation date: (13/09/01 %T) 301 * 302 * @param arg1 java.lang.String 303 * @param arg2 java.lang.String 304 * @param arg3 java.lang.String 305 * 306 * @exception java.lang.IllegalArgumentException The exception description. 307 */ parse_Gmsh(Token[] param, int lineno)308 public void parse_Gmsh(Token[] param, int lineno) 309 throws java.text.ParseException 310 { 311 int index; 312 313 if (param[2].is_a_number()) { 314 // Extract number of material 315 index = (int)(param[2].getn()/1E3); 316 index %= 100; 317 318 name = new String("MAT_" + index); 319 name_is_set = true; 320 } 321 else 322 throw new java.text.ParseException( 323 "Illegal identification number of material " + param[1].getn(),lineno 324 ); 325 326 327 // Hard code material propoerties 328 youngs_modulus = 210; 329 E_is_set = true; 330 youngs_modulus_initialized = true; 331 332 nu = 0.3; 333 NU_is_set = true; 334 nu_initialized = true; 335 336 density = 7.8e-6; 337 density_is_set = true; 338 339 } 340 341 342 343 344 345 346 347 /** 348 * Checks that all mandatory variables have been set 349 */ checkIndata()350 public void checkIndata() 351 throws IllegalArgumentException 352 { 353 if (! E_is_set) { 354 throw new IllegalArgumentException( 355 "No Youngs modulus defined for this material: " + name 356 ); 357 } 358 359 if (! density_is_set) { 360 throw new IllegalArgumentException("No Density defined for this material: " + name); 361 } 362 363 if (! NU_is_set) { 364 throw new IllegalArgumentException( 365 "No Poissons_constant defined for this material: " + name 366 ); 367 } 368 } 369 370 /** 371 * This method is used to create the lines needed in the result file. The 372 * method returns a string which is printed directly. However, due to the 373 * fact that the line may be different depending on what is requested to 374 * be printed and that the number of methods should be kept down, the 375 * first parameter here is a control parameter. This parameter describes 376 * what should be printed. The second parameter is a required input when 377 * gauss point results are to be printed. Creation date: (09/12/01 %T) 378 * 379 * @param ctrl The control number to say if a header of result file is to 380 * be printed. 381 * 382 * @return java.lang.String 383 */ print_Fembic(int ctrl)384 public String print_Fembic(int ctrl) { 385 String out; 386 387 switch (ctrl) { 388 389 case Element.MESH: 390 391 /* Print the element number and connected nodes */ 392 out = new String( 393 name + "\tE = " + youngs_modulus 394 + "\tRHO = " + density 395 + "\tNU = " + nu 396 ); 397 398 399 if (failure_strain_is_set) 400 out += " failure_strain = " + failure_strain; 401 402 if (failure_stress_is_set) 403 out += " failure_stress = " + failure_stress; 404 405 out += "\n"; 406 407 return out; 408 409 410 default: 411 return new String(""); 412 } 413 } 414 415 416 417 418 /** 419 * Insert the method's description here. Creation date: (16/12/01 %T) 420 */ setInitialConditions()421 public void setInitialConditions() { 422 double c; 423 double c2; 424 double G; 425 426 // Check all indata is there 427 if (! density_is_set) { 428 throw new IllegalArgumentException("Material density not intialized"); 429 } else if (! name_is_set) { 430 throw new IllegalArgumentException("Material name not intialized"); 431 } else if (! nu_initialized) { 432 throw new IllegalArgumentException("Material Poisson constant (nu) not intialized"); 433 } else if (! youngs_modulus_initialized) { 434 throw new IllegalArgumentException("Material youngs_modulus not intialized"); 435 } 436 437 // Initialize all the stiffness matrices 438 c = youngs_modulus / ((1.0 + nu) * (1.0 - (2.0 * nu))); 439 c2 = youngs_modulus / (1.0 - (nu * nu)); 440 G = youngs_modulus / (2.0 * (1.0 + nu)); 441 442 // 443 stiffness_matrix_3d.set(0, 0, (1.0 - nu) * c); 444 stiffness_matrix_3d.set(0, 1, nu * c); 445 stiffness_matrix_3d.set(0, 2, nu * c); 446 stiffness_matrix_3d.set(0, 3, 0); 447 stiffness_matrix_3d.set(0, 4, 0); 448 stiffness_matrix_3d.set(0, 5, 0); 449 450 // 451 stiffness_matrix_3d.set(1, 0, nu * c); 452 stiffness_matrix_3d.set(1, 1, (1.0 - nu) * c); 453 stiffness_matrix_3d.set(1, 2, nu * c); 454 stiffness_matrix_3d.set(1, 3, 0); 455 stiffness_matrix_3d.set(1, 4, 0); 456 stiffness_matrix_3d.set(1, 5, 0); 457 458 // 459 stiffness_matrix_3d.set(2, 0, nu * c); 460 stiffness_matrix_3d.set(2, 1, nu * c); 461 stiffness_matrix_3d.set(2, 2, (1.0 - nu) * c); 462 stiffness_matrix_3d.set(2, 3, 0); 463 stiffness_matrix_3d.set(2, 4, 0); 464 stiffness_matrix_3d.set(2, 5, 0); 465 466 // 467 stiffness_matrix_3d.set(3, 0, 0); 468 stiffness_matrix_3d.set(3, 1, 0); 469 stiffness_matrix_3d.set(3, 2, 0); 470 stiffness_matrix_3d.set(3, 3, G / 2.0); 471 stiffness_matrix_3d.set(3, 4, 0); 472 stiffness_matrix_3d.set(3, 5, 0); 473 474 // 475 stiffness_matrix_3d.set(4, 0, 0); 476 stiffness_matrix_3d.set(4, 1, 0); 477 stiffness_matrix_3d.set(4, 2, 0); 478 stiffness_matrix_3d.set(4, 3, 0); 479 stiffness_matrix_3d.set(4, 4, G / 2.0); 480 stiffness_matrix_3d.set(4, 5, 0); 481 482 // 483 stiffness_matrix_3d.set(5, 0, 0); 484 stiffness_matrix_3d.set(5, 1, 0); 485 stiffness_matrix_3d.set(5, 2, 0); 486 stiffness_matrix_3d.set(5, 3, 0); 487 stiffness_matrix_3d.set(5, 4, 0); 488 stiffness_matrix_3d.set(5, 5, G / 2.0); 489 490 // Plane stress matrix 491 stiffness_matrix_plane_stress.set(0, 0, c2); 492 stiffness_matrix_plane_stress.set(0, 1, nu * c2); 493 stiffness_matrix_plane_stress.set(0, 2, 0); 494 stiffness_matrix_plane_stress.set(0, 3, 0); 495 stiffness_matrix_plane_stress.set(0, 4, 0); 496 stiffness_matrix_plane_stress.set(0, 5, 0); 497 498 // 499 stiffness_matrix_plane_stress.set(1, 0, nu * c2); 500 stiffness_matrix_plane_stress.set(1, 1, c2); 501 stiffness_matrix_plane_stress.set(1, 2, 0); 502 stiffness_matrix_plane_stress.set(1, 3, 0); 503 stiffness_matrix_plane_stress.set(1, 4, 0); 504 stiffness_matrix_plane_stress.set(1, 5, 0); 505 506 // 507 stiffness_matrix_plane_stress.set(2, 0, 0); 508 stiffness_matrix_plane_stress.set(2, 1, 0); 509 stiffness_matrix_plane_stress.set(2, 2, 0); 510 stiffness_matrix_plane_stress.set(2, 3, 0); 511 stiffness_matrix_plane_stress.set(2, 4, 0); 512 stiffness_matrix_plane_stress.set(2, 5, 0); 513 514 // 515 stiffness_matrix_plane_stress.set(3, 0, 0); 516 stiffness_matrix_plane_stress.set(3, 1, 0); 517 stiffness_matrix_plane_stress.set(3, 2, 0); 518 stiffness_matrix_plane_stress.set(3, 3, ((1 - nu) * c2) / 2); 519 stiffness_matrix_plane_stress.set(3, 4, 0); 520 stiffness_matrix_plane_stress.set(3, 5, 0); 521 522 // 523 stiffness_matrix_plane_stress.set(4, 0, 0); 524 stiffness_matrix_plane_stress.set(4, 1, 0); 525 stiffness_matrix_plane_stress.set(4, 2, 0); 526 stiffness_matrix_plane_stress.set(4, 3, 0); 527 stiffness_matrix_plane_stress.set( 528 4, 4, ((5.0 / 6.0) * (1 - nu) * c2) / 2 529 ); 530 stiffness_matrix_plane_stress.set(4, 5, 0); 531 532 // 533 stiffness_matrix_plane_stress.set(5, 0, 0); 534 stiffness_matrix_plane_stress.set(5, 1, 0); 535 stiffness_matrix_plane_stress.set(5, 2, 0); 536 stiffness_matrix_plane_stress.set(5, 3, 0); 537 stiffness_matrix_plane_stress.set(5, 4, 0); 538 stiffness_matrix_plane_stress.set( 539 5, 5, ((5.0 / 6.0) * (1 - nu) * c2) / 2 540 ); 541 } 542 543 /** 544 * Insert the method's description here. Creation date: (06/01/02 %T) 545 * 546 * @return double 547 */ wavespeedOneDimensional(double param, double param2)548 public double wavespeedOneDimensional(double param, double param2) { 549 return Math.sqrt(youngs_modulus / density); 550 } 551 552 /** 553 * Insert the method's description here. Creation date: (06/01/02 %T) 554 * 555 * @return double 556 */ wavespeedThreeDimensional(double param, double param2)557 public double wavespeedThreeDimensional(double param, double param2) { 558 return Math.sqrt( 559 (youngs_modulus * (1 - nu)) / (density * (1 + nu) * (1 - (2 * nu))) 560 ); 561 } 562 563 /** 564 * Insert the method's description here. Creation date: (06/01/02 %T) 565 * 566 * @return double 567 */ wavespeedTwoDimensional(double param, double param2)568 public double wavespeedTwoDimensional(double param, double param2) { 569 return Math.sqrt(youngs_modulus / (density * (1 - (nu * nu)))); 570 } 571 } 572 573