FiniteDeformationEP3D.cppGo to the documentation of this file.00001 //=============================================================================== 00002 00003 //# COPYRIGHT (C): Woody's license (by BJ): 00004 00005 // ``This source code is Copyrighted in 00006 00007 // U.S., for an indefinite period, and anybody 00008 00009 // caught using it without our permission, will be 00010 00011 // mighty good friends of ourn, cause we don't give 00012 00013 // a darn. Hack it. Compile it. Debug it. Run it. 00014 00015 // Yodel it. Enjoy it. We wrote it, that's all we 00016 00017 // wanted to do.'' 00018 00019 // 00020 00021 //# PROJECT: Object Oriented Finite Element Program 00022 00023 //# PURPOSE: Finite Deformation Hyper-Elastic classes 00024 00025 //# CLASS: 00026 00027 //# 00028 00029 //# VERSION: 0.6_(1803398874989) (golden section) 00030 00031 //# LANGUAGE: C++ 00032 00033 //# TARGET OS: all... 00034 00035 //# DESIGN: Zhao Cheng, Boris Jeremic (jeremic@ucdavis.edu) 00036 00037 //# PROGRAMMER(S): Zhao Cheng, Boris Jeremic 00038 00039 //# 00040 00041 //# 00042 00043 //# DATE: July 2004 00044 00045 //# UPDATE HISTORY: July25th2004, Zhao added the Algorithmic Tangent Stiffness 00046 00047 //# for nearly quadratic convergence rate 00048 00049 //=============================================================================== 00050 00051 00052 00053 #ifndef FiniteDeformationEP3D_CPP 00054 00055 #define FiniteDefornationEP3D_CPP 00056 00057 00058 00059 #include "FiniteDeformationEP3D.h" 00060 00061 00062 00063 const int Max_Iter = 30; 00064 00065 const double tolerance = 1.0e-6; 00066 00067 00068 00069 tensor tensorZ2(2, def_dim_2, 0.0); 00070 00071 tensor tensorI2("I", 2, def_dim_2); 00072 00073 tensor tensorZ4(4, def_dim_4, 0.0); 00074 00075 00076 00077 stresstensor FiniteDeformationEP3D::static_stress; 00078 00079 straintensor FiniteDeformationEP3D::static_strain; 00080 00081 00082 00083 //Constructor 00---------------------------------------------------------------------------------- 00084 00085 FiniteDeformationEP3D::FiniteDeformationEP3D() 00086 00087 :NDMaterial(0, 0) 00088 00089 { 00090 00091 fde3d = 0; 00092 00093 fdy = 0; 00094 00095 fdf = 0; 00096 00097 fdEvolutionS = 0; 00098 00099 fdEvolutionT = 0; 00100 00101 fdeps = 0; 00102 00103 00104 00105 int err; 00106 00107 err = this->revertToStart(); 00108 00109 } 00110 00111 00112 00113 // Constructor 01--------------------------------------------------------------------------------- 00114 00115 FiniteDeformationEP3D::FiniteDeformationEP3D(int tag, 00116 00117 NDMaterial *fde3d_in, 00118 00119 fdYield *fdy_in, 00120 00121 fdFlow *fdf_in, 00122 00123 fdEvolution_S *fdEvolutionS_in, 00124 00125 fdEvolution_T *fdEvolutionT_in) 00126 00127 :NDMaterial(tag, ND_TAG_FiniteDeformationEP3D) 00128 00129 { 00130 00131 if (fde3d_in) 00132 00133 fde3d = fde3d_in->getCopy(); 00134 00135 else { 00136 00137 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdElastic3D\n"; 00138 00139 exit (-1); 00140 00141 } 00142 00143 00144 00145 if (fdy_in) 00146 00147 fdy = fdy_in->newObj(); 00148 00149 else { 00150 00151 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdYield\n"; 00152 00153 exit (-1); 00154 00155 } 00156 00157 00158 00159 if (fdf_in) 00160 00161 fdf = fdf_in->newObj(); 00162 00163 else { 00164 00165 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdFlow\n"; 00166 00167 exit (-1); 00168 00169 } 00170 00171 00172 00173 if (fdEvolutionS_in) 00174 00175 fdEvolutionS = fdEvolutionS_in->newObj(); 00176 00177 else 00178 00179 fdEvolutionS = 0; 00180 00181 00182 00183 if (fdEvolutionT_in) 00184 00185 fdEvolutionT = fdEvolutionT_in->newObj(); 00186 00187 else 00188 00189 fdEvolutionT = 0; 00190 00191 00192 00193 //if (fdeps == 0) 00194 00195 fdeps = new FDEPState(); 00196 00197 00198 00199 } 00200 00201 00202 00203 // Constructor 02--------------------------------------------------------------------------------- 00204 00205 FiniteDeformationEP3D::FiniteDeformationEP3D(int tag, 00206 00207 NDMaterial *fde3d_in, 00208 00209 fdYield *fdy_in, 00210 00211 fdFlow *fdf_in, 00212 00213 fdEvolution_S *fdEvolutionS_in) 00214 00215 :NDMaterial(tag, ND_TAG_FiniteDeformationEP3D) 00216 00217 { 00218 00219 if (fde3d_in) 00220 00221 fde3d = fde3d_in->getCopy(); 00222 00223 else { 00224 00225 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdElastic3D\n"; 00226 00227 exit (-1); 00228 00229 } 00230 00231 00232 00233 if (fdy_in) 00234 00235 fdy = fdy_in->newObj(); 00236 00237 else { 00238 00239 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdYield\n"; 00240 00241 exit (-1); 00242 00243 } 00244 00245 00246 00247 if (fdf_in) 00248 00249 fdf = fdf_in->newObj(); 00250 00251 else { 00252 00253 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdFlow\n"; 00254 00255 exit (-1); 00256 00257 } 00258 00259 00260 00261 if (fdEvolutionS_in) 00262 00263 fdEvolutionS = fdEvolutionS_in->newObj(); 00264 00265 else 00266 00267 fdEvolutionS = 0; 00268 00269 00270 00271 fdEvolutionT = 0; 00272 00273 00274 00275 //if (fdeps == 0) 00276 00277 fdeps = new FDEPState(); 00278 00279 00280 00281 } 00282 00283 00284 00285 // Constructor 03--------------------------------------------------------------------------------- 00286 00287 FiniteDeformationEP3D::FiniteDeformationEP3D(int tag, 00288 00289 NDMaterial *fde3d_in, 00290 00291 fdYield *fdy_in, 00292 00293 fdFlow *fdf_in, 00294 00295 fdEvolution_T *fdEvolutionT_in) 00296 00297 :NDMaterial(tag, ND_TAG_FiniteDeformationEP3D) 00298 00299 { 00300 00301 if (fde3d_in) 00302 00303 fde3d = fde3d_in->getCopy(); 00304 00305 else { 00306 00307 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdElastic3D\n"; 00308 00309 exit (-1); 00310 00311 } 00312 00313 00314 00315 if (fdy_in) 00316 00317 fdy = fdy_in->newObj(); 00318 00319 else { 00320 00321 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdYield\n"; 00322 00323 exit (-1); 00324 00325 } 00326 00327 00328 00329 if (fdf_in) 00330 00331 fdf = fdf_in->newObj(); 00332 00333 else { 00334 00335 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdFlow\n"; 00336 00337 exit (-1); 00338 00339 } 00340 00341 00342 00343 fdEvolutionS = 0; 00344 00345 00346 00347 if (fdEvolutionT_in) 00348 00349 fdEvolutionT = fdEvolutionT_in->newObj(); 00350 00351 else 00352 00353 fdEvolutionT = 0; 00354 00355 00356 00357 //if (fdeps == 0) 00358 00359 fdeps = new FDEPState(); 00360 00361 00362 00363 } 00364 00365 00366 00367 // Constructor 04--------------------------------------------------------------------------------- 00368 00369 FiniteDeformationEP3D::FiniteDeformationEP3D(int tag, 00370 00371 NDMaterial *fde3d_in, 00372 00373 fdYield *fdy_in, 00374 00375 fdFlow *fdf_in) 00376 00377 :NDMaterial(tag, ND_TAG_FiniteDeformationEP3D) 00378 00379 { 00380 00381 if (fde3d_in) 00382 00383 fde3d = fde3d_in->getCopy(); 00384 00385 else { 00386 00387 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdElastic3D\n"; 00388 00389 exit (-1); 00390 00391 } 00392 00393 00394 00395 if (fdy_in) 00396 00397 fdy = fdy_in->newObj(); 00398 00399 else { 00400 00401 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdYield\n"; 00402 00403 exit (-1); 00404 00405 } 00406 00407 00408 00409 if (fdf_in) 00410 00411 fdf = fdf_in->newObj(); 00412 00413 else { 00414 00415 opserr << "FiniteDeformationEP3D:: FiniteDeformationEP3D failed to construct the fdFlow\n"; 00416 00417 exit (-1); 00418 00419 } 00420 00421 00422 00423 fdEvolutionS = 0; 00424 00425 00426 00427 fdEvolutionT = 0; 00428 00429 00430 00431 //if (fdeps == 0) 00432 00433 fdeps = new FDEPState(); 00434 00435 00436 00437 } 00438 00439 00440 00441 // Destructor------------------------------------------------------------------------------------- 00442 00443 FiniteDeformationEP3D::~FiniteDeformationEP3D() 00444 00445 { 00446 00447 if (fde3d) 00448 00449 delete fde3d; 00450 00451 00452 00453 if (fdy) 00454 00455 delete fdy; 00456 00457 00458 00459 if (fdf) 00460 00461 delete fdf; 00462 00463 00464 00465 if (fdEvolutionS) 00466 00467 delete fdEvolutionS; 00468 00469 00470 00471 if (fdEvolutionT) 00472 00473 delete fdEvolutionT; 00474 00475 00476 00477 if (fdeps) 00478 00479 delete fdeps; 00480 00481 } 00482 00483 00484 00485 //---------------------------------------------------------------------- 00486 00487 double FiniteDeformationEP3D::getRho(void) 00488 00489 { 00490 00491 return fde3d->getRho(); 00492 00493 } 00494 00495 00496 00497 //---------------------------------------------------------------------- 00498 00499 int FiniteDeformationEP3D::setTrialF(const straintensor &f) 00500 00501 { 00502 00503 F.Initialize(f); 00504 00505 00506 00507 // Choose One: 00508 00509 //return ImplicitAlgorithm(); 00510 00511 return SemiImplicitAlgorithm(); 00512 00513 00514 00515 } 00516 00517 00518 00519 //---------------------------------------------------------------------- 00520 00521 int FiniteDeformationEP3D::setTrialFIncr(const straintensor &df) 00522 00523 { 00524 00525 return this->setTrialF(this->getF() + df); 00526 00527 } 00528 00529 00530 00531 //---------------------------------------------------------------------- 00532 00533 const tensor& FiniteDeformationEP3D::getTangentTensor(void) 00534 00535 { 00536 00537 return iniTangent; 00538 00539 } 00540 00541 00542 00543 //---------------------------------------------------------------------- 00544 00545 const straintensor& FiniteDeformationEP3D::getStrainTensor(void) 00546 00547 { 00548 00549 return iniGreen; 00550 00551 } 00552 00553 00554 00555 //---------------------------------------------------------------------- 00556 00557 const stresstensor& FiniteDeformationEP3D::getStressTensor(void) 00558 00559 { 00560 00561 return iniPK2; 00562 00563 } 00564 00565 00566 00567 //---------------------------------------------------------------------- 00568 00569 const straintensor& FiniteDeformationEP3D::getF(void) 00570 00571 { 00572 00573 return F; 00574 00575 } 00576 00577 00578 00579 //---------------------------------------------------------------------- 00580 00581 const straintensor& FiniteDeformationEP3D::getFp(void) 00582 00583 { 00584 00585 straintensor fpp = fdeps->getFpInVar(); 00586 00587 static_strain.Initialize(fpp); 00588 00589 00590 00591 return static_strain; 00592 00593 } 00594 00595 00596 00597 //---------------------------------------------------------------------- 00598 00599 int FiniteDeformationEP3D::commitState(void) 00600 00601 { 00602 00603 return fdeps->commitState(); 00604 00605 } 00606 00607 00608 00609 //---------------------------------------------------------------------- 00610 00611 int FiniteDeformationEP3D::revertToLastCommit(void) 00612 00613 { 00614 00615 return fdeps->revertToLastCommit(); 00616 00617 } 00618 00619 00620 00621 //---------------------------------------------------------------------- 00622 00623 int FiniteDeformationEP3D::revertToStart(void) 00624 00625 { 00626 00627 return fdeps->revertToStart(); 00628 00629 } 00630 00631 00632 00633 //---------------------------------------------------------------------- 00634 00635 NDMaterial* FiniteDeformationEP3D::getCopy (void) 00636 00637 { 00638 00639 NDMaterial* tmp = new FiniteDeformationEP3D( 00640 00641 this->getTag(), 00642 00643 this->getFDE3D(), 00644 00645 this->getFDY(), 00646 00647 this->getFDF(), 00648 00649 this->getFDEvolutionS(), 00650 00651 this->getFDEvolutionT() ); 00652 00653 00654 00655 return tmp; 00656 00657 } 00658 00659 00660 00661 //---------------------------------------------------------------------- 00662 00663 NDMaterial* FiniteDeformationEP3D::getCopy (const char *code) 00664 00665 { 00666 00667 if ( strcmp(code,"FiniteDeformationEP3D") == 0 00668 00669 || strcmp(code,"FDEP3D") == 0 ) 00670 00671 { 00672 00673 NDMaterial* tmp = new FiniteDeformationEP3D ( 00674 00675 this->getTag(), 00676 00677 this->getFDE3D(), 00678 00679 this->getFDY(), 00680 00681 this->getFDF(), 00682 00683 this->getFDEvolutionS(), 00684 00685 this->getFDEvolutionT() ); 00686 00687 00688 00689 return tmp; 00690 00691 } 00692 00693 else 00694 00695 { 00696 00697 opserr << "FiniteDeformationEP3D::getCopy fainled:" << code << "\n"; 00698 00699 exit(-1); 00700 00701 } 00702 00703 } 00704 00705 00706 00707 //---------------------------------------------------------------------- 00708 00709 const char* FiniteDeformationEP3D::getType (void) const 00710 00711 { 00712 00713 return "ThreeDimentionalFD"; 00714 00715 } 00716 00717 00718 00719 00720 00722 00723 //int FiniteDeformationEP3D::getOrder (void) const 00724 00725 //{ 00726 00727 // return 6; 00728 00729 //} 00730 00731 00732 00733 //---------------------------------------------------------------------- 00734 00735 int FiniteDeformationEP3D::sendSelf(int commitTag, Channel &theChannel) 00736 00737 { 00738 00739 // Not yet implemented 00740 00741 return 0; 00742 00743 } 00744 00745 00746 00747 //---------------------------------------------------------------------- 00748 00749 int FiniteDeformationEP3D::recvSelf(int commitTag, Channel &theChannel, FEM_ObjectBroker &theBroker) 00750 00751 { 00752 00753 // Not yet implemented 00754 00755 return 0; 00756 00757 } 00758 00759 00760 00761 //---------------------------------------------------------------------- 00762 00763 void FiniteDeformationEP3D::Print(OPS_Stream &s, int flag) 00764 00765 { 00766 00767 s << (*this); 00768 00769 } 00770 00771 00772 00773 //---------------------------------------------------------------------- 00774 00775 const stresstensor& FiniteDeformationEP3D::getPK1StressTensor(void) 00776 00777 { 00778 00779 stresstensor thisSPKStress; 00780 00781 00782 00783 thisSPKStress = this->getStressTensor(); 00784 00785 static_stress = F("kJ") * thisSPKStress("IJ"); 00786 00787 00788 00789 return static_stress; 00790 00791 } 00792 00793 00794 00795 //---------------------------------------------------------------------- 00796 00797 const stresstensor& FiniteDeformationEP3D::getCauchyStressTensor(void) 00798 00799 { 00800 00801 stresstensor thisSPKStress; 00802 00803 double JJ = F.determinant(); 00804 00805 00806 00807 thisSPKStress = this->getStressTensor(); 00808 00809 static_stress = F("iK") * thisSPKStress("JK"); 00810 00811 static_stress = static_stress("iJ") * F("mJ"); 00812 00813 static_stress = static_stress *(1.0/JJ); 00814 00815 00816 00817 return static_stress; 00818 00819 } 00820 00821 00822 00823 //---------------------------------------------------------------------- 00824 00825 NDMaterial * FiniteDeformationEP3D::getFDE3D() const 00826 00827 { 00828 00829 return fde3d; 00830 00831 } 00832 00833 00834 00835 //---------------------------------------------------------------------- 00836 00837 fdYield * FiniteDeformationEP3D::getFDY() const 00838 00839 { 00840 00841 return fdy; 00842 00843 } 00844 00845 00846 00847 //---------------------------------------------------------------------- 00848 00849 fdFlow * FiniteDeformationEP3D::getFDF() const 00850 00851 { 00852 00853 return fdf; 00854 00855 } 00856 00857 00858 00859 //---------------------------------------------------------------------- 00860 00861 fdEvolution_S * FiniteDeformationEP3D::getFDEvolutionS() const 00862 00863 { 00864 00865 return fdEvolutionS; 00866 00867 } 00868 00869 00870 00871 //---------------------------------------------------------------------- 00872 00873 fdEvolution_T * FiniteDeformationEP3D::getFDEvolutionT() const 00874 00875 { 00876 00877 return fdEvolutionT; 00878 00879 } 00880 00881 00882 00883 //---------------------------------------------------------------------- 00884 00885 FDEPState * FiniteDeformationEP3D::getFDEPState() const 00886 00887 { 00888 00889 return fdeps; 00890 00891 } 00892 00893 00894 00895 00896 00897 //---------------------------------------------------------------------- 00898 00900 00901 //---------------------------------------------------------------------- 00902 00903 int FiniteDeformationEP3D::ImplicitAlgorithm() 00904 00905 { 00906 00907 // Initializing 00908 00909 double yieldfun = 0.0; // trial value of yield function 00910 00911 double D_gamma = 0.0; // consistency parameter /Delta_{gamma} 00912 00913 double d_gamma = 0.0; // increment of consistency parameter /delta_{gamma} 00914 00915 int iter_counter = 0; 00916 00917 00918 00919 straintensor res_Ee = tensorZ2; // residual of intermediate Ee 00920 00921 straintensor res_eta = tensorZ2; // norm of residual of eta 00922 00923 00924 00925 double res_xi = 0.0; // residual of strain like internal variable 00926 00927 double res_norm_eta = 0.0; // norm of residual of eta 00928 00929 double res_norm_Ee = 0.0; // norm of residual of intermediate Ee 00930 00931 double residual = 0.0; // residual of Ee and xi and eta 00932 00933 00934 00935 straintensor Fp = tensorI2;; // Plastic deformation gradient 00936 00937 double xi = 0.0; // strain like internal isotropic variable 00938 00939 double q = 0.0; // stress like internal isotropic variable 00940 00941 straintensor eta; // strain like internal kinematic variable 00942 00943 stresstensor a; // stress like internal kinematic variable 00944 00945 00946 00947 straintensor Fpinv = tensorI2; // inverse of Fp 00948 00949 straintensor Ce = tensorI2; // intermediate C 00950 00951 straintensor Ee = tensorZ2; // intermediate Ee 00952 00953 00954 00955 straintensor Fp_n = tensorI2; // Plastic deformation gradient at time n 00956 00957 straintensor Fp_ninv = tensorI2; // Plastic deformation gradient at time n 00958 00959 straintensor Ee_n = tensorZ2; // Ee at the incremental step n, calculated from Fp_n 00960 00961 double xi_n = 0.0; // xi at the incremental step n, known 00962 00963 straintensor eta_n = tensorZ2; // eta at the incremental step n, known 00964 00965 00966 00967 stresstensor Mtensor = tensorZ2; // --> dFl/dT 00968 00969 stresstensor MCtensor = tensorZ2; // --> dFl/dS 00970 00971 tensor Ltensor = tensorZ4 ; // Tangent tensor in the intermediate configuration 00972 00973 tensor LATStensor = tensorZ4; // Consistent tangent tensor in the intermediate configuration 00974 00975 00976 00977 double nscalar = 0.0; // --> dFl/d(xi) 00978 00979 double Kscalar = 0.0; // Isotropic hardening modoulus 00980 00981 straintensor ntensor; // --> dFl/d(eta) 00982 00983 tensor Ktensor = tensorZ4; // Kinematic hardening modoulus 00984 00985 00986 00987 stresstensor dyods = tensorZ2 ; // --> dY/d(stress) 00988 00989 double dyodq = 0.0; // --> dY/d(xi) 00990 00991 stresstensor dyoda = tensorZ2 ; // --> dY/d(eta) 00992 00993 00994 00995 tensor dMods = tensorZ4; // --> dM/d(stress), from d2Fl/d(stress)d(stress) 00996 00997 tensor dMCods = tensorZ4; // --> dMs/d(stress) 00998 00999 01000 01001 double dnsodq = 0.0; // --> d2Fl/dqdq 01002 01003 tensor dntoda = tensorZ4; // --> d2Fl/dada 01004 01005 01006 01007 straintensor D_Ee = tensorZ2; 01008 01009 double D_xi = 0.0; 01010 01011 straintensor D_eta = tensorZ2; 01012 01013 01014 01015 tensor tensorI4 = (tensorI2("ij")*tensorI2("kl")).transpose0110(); 01016 01017 01018 01019 tensor A11 = tensorI4; tensor A12 = tensorZ2; tensor A13 = tensorZ4; 01020 01021 tensor A21 = tensorZ2; double a22 = 1.0; tensor A23 = tensorZ2; 01022 01023 tensor A31 = tensorZ4; tensor A32 = tensorZ2; tensor A33 = tensorI4; 01024 01025 01026 01027 BJmatrix C99(9, 9, 0.0); 01028 01029 BJmatrix CCC(19, 19, 0.0); 01030 01031 01032 01033 tensor tensorTemp0; 01034 01035 tensor tensorTemp1; 01036 01037 tensor tensorTemp2; 01038 01039 tensor tensorTemp3; 01040 01041 tensor tensorTemp4; 01042 01043 tensor tensorTemp5; 01044 01045 tensor tensorTemp6; 01046 01047 tensor tensorTemp7; 01048 01049 tensor tensorTemp8; 01050 01051 tensor tensorTemp9; 01052 01053 double temp1 = 0.0; 01054 01055 double temp2 = 0.0; 01056 01057 double temp3 = 0.0; 01058 01059 double temp4 = 0.0; 01060 01061 double temp5 = 0.0; 01062 01063 01064 01065 tensor LM = tensorZ4; // For Mandel Tangent Stiffness 01066 01067 stresstensor B_Mandel = tensorZ2; // For Mandel stress 01068 01069 01070 01071 // Read the previous incremental step history variables 01072 01073 Fp = fdeps->getCommitedFpInVar(); 01074 01075 Fp_n = Fp; 01076 01077 Fpinv = Fp.inverse(); 01078 01079 Fp_ninv = Fp_n.inverse(); 01080 01081 Fe = F("ij")*Fpinv("jk"); Fe.null_indices(); 01082 01083 Ce = Fe("ki")*Fe("kj"); Ce.null_indices(); 01084 01085 Ee = (Ce - tensorI2) * 0.5; 01086 01087 Ee_n = Ee; 01088 01089 01090 01091 if ( fdEvolutionS ) { 01092 01093 xi = fdeps->getCommitedStrainLikeInVar(); 01094 01095 xi_n = xi; 01096 01097 q = fdeps->getCommitedStressLikeInVar(); 01098 01099 } 01100 01101 01102 01103 if ( fdEvolutionT ) { 01104 01105 eta = fdeps->getCommitedStrainLikeKiVar(); 01106 01107 eta_n = eta; 01108 01109 a = fdeps->getCommitedStressLikeKiVar(); 01110 01111 } 01112 01113 01114 01115 // Return stress from finite deformation elastic model 01116 01117 fde3d->setTrialC(Ce); // Note: It is C, not F!!! 01118 01119 B_PK2 = fde3d->getStressTensor(); 01120 01121 01122 01123 B_Mandel = Ce("ik")*B_PK2("kj"); // Mandel Stress 01124 01125 B_Mandel.null_indices(); 01126 01127 01128 01129 // Evaluate the value of yield function 01130 01131 yieldfun = fdy->Yd(B_Mandel, *fdeps); 01132 01133 //printf("\nY0= %e\n", yieldfun); 01134 01135 01136 01137 if ( yieldfun > (fdy->getTolerance()) ) { // start of plastic part 01138 01139 D_gamma = 0.0; 01140 01141 d_gamma = 0.0; 01142 01143 iter_counter = 0; 01144 01145 01146 01147 do { // beginning of do - while 01148 01149 01150 01151 // Return Symmetric tensor M and n_scalar and n_tensor 01152 01153 Mtensor = fdf->dFods(B_Mandel, *fdeps); 01154 01155 Mtensor = Ce("ik")*Mtensor("kj"); 01156 01157 Mtensor.null_indices(); 01158 01159 MCtensor = (Mtensor + Mtensor.transpose11()) * 0.5; 01160 01161 01162 01163 // Return tangent variables 01164 01165 Ltensor = fde3d->getTangentTensor(); 01166 01167 01168 01169 tensorTemp1 = tensorI2("ij") *B_PK2("mn"); 01170 01171 tensorTemp1.null_indices(); 01172 01173 tensorTemp2 = Ce("ik") *Ltensor("kjmn"); 01174 01175 tensorTemp2.null_indices(); 01176 01177 LM = tensorTemp2 + tensorTemp1.transpose0110() + tensorTemp1.transpose0111(); 01178 01179 01180 01181 if ( fdEvolutionS) { 01182 01183 nscalar = fdf->dFodq(B_Mandel, *fdeps); 01184 01185 Kscalar = fdEvolutionS->HModulus(B_Mandel, *fdeps); 01186 01187 dyodq = fdy->dYodq(B_Mandel, *fdeps); 01188 01189 dnsodq = fdf->d2Fodqdq(B_Mandel, *fdeps); 01190 01191 } 01192 01193 01194 01195 if ( fdEvolutionT) { 01196 01197 ntensor = fdf->dFoda(B_Mandel, *fdeps); 01198 01199 Ktensor = fdEvolutionT->HModulus(B_Mandel, *fdeps); 01200 01201 dyoda = fdy->dYoda(B_Mandel, *fdeps); 01202 01203 dntoda = fdf->d2Fodada(B_Mandel, *fdeps); 01204 01205 } 01206 01207 01208 01209 // Return tensor Axx 01210 01211 dyods = fdy->dYods(B_Mandel, *fdeps); 01212 01213 dMods = fdf->d2Fodsds(B_Mandel, *fdeps); 01214 01215 dMCods = Ce("ik")*dMods("kjmn"); 01216 01217 dMCods.null_indices(); 01218 01219 dMCods = (dMCods + dMCods.transpose1100()) * 0.5; 01220 01221 01222 01223 A11 = dMCods("ijmn")*LM("mnkl"); 01224 01225 A11.null_indices(); 01226 01227 A11 = (A11 *D_gamma) + tensorI4; //A11 01228 01229 01230 01231 if ( fdEvolutionS ) { 01232 01233 a22 += (dnsodq *Kscalar*D_gamma); //A22 01234 01235 01236 01237 tensorTemp1 = fdf->d2Fodsdq(B_Mandel, *fdeps); 01238 01239 A12 = Ce("ik") * tensorTemp1("kj"); 01240 01241 A12.null_indices(); 01242 01243 A12 = (A12 + A12.transpose11()) * (0.5*Kscalar*D_gamma); //A12 01244 01245 01246 01247 A21 = tensorTemp1("ij") * LM("ijkl"); 01248 01249 A21.null_indices(); 01250 01251 A21 = A21 *D_gamma; //A21 01252 01253 } 01254 01255 01256 01257 if ( fdEvolutionT ) { 01258 01259 A33 = dntoda("ijmn")*Ktensor("mnkl"); 01260 01261 A33.null_indices(); 01262 01263 A33 += (A33 *D_gamma); //A33 01264 01265 01266 01267 tensorTemp2 = fdf->d2Fodsda(B_Mandel, *fdeps); 01268 01269 A13 = Ce("ik") * tensorTemp2("kjmn"); 01270 01271 A13.null_indices(); 01272 01273 A13 = (A13 + A13.transpose1100()) * (0.5*D_gamma); 01274 01275 A13 = A13("ijmn") * Ktensor("mnkl"); 01276 01277 A13.null_indices(); //A13 01278 01279 A31 = tensorTemp2("ijmn") * LM("mnkl"); 01280 01281 A31.null_indices(); 01282 01283 A31 = A31 *D_gamma; //A31 01284 01285 } 01286 01287 01288 01289 if ( fdEvolutionS && fdEvolutionT ) { 01290 01291 tensorTemp3 = fdf->d2Fodqda(B_Mandel, *fdeps); 01292 01293 A23 = tensorTemp3("ij") * Ktensor("ijkl"); 01294 01295 A23.null_indices(); 01296 01297 A23 = A23 *D_gamma; //A23 01298 01299 tensorTemp2 = fdf->d2Fodsda(B_Mandel, *fdeps); 01300 01301 A32 = tensorTemp2 * (Kscalar*D_gamma); //A32 01302 01303 } 01304 01305 01306 01307 int i, j; 01308 01309 // CCC: Tensor -> Matrix 01310 01311 C99 = A11.BJtensor2BJmatrix_2(); 01312 01313 for (i=10; i<=19; i++) 01314 01315 CCC.val(i,i) = 1.0; 01316 01317 01318 01319 if ( fdEvolutionS ) { 01320 01321 CCC.val(1,10) = A12.cval(1,1); 01322 01323 CCC.val(2,10) = A12.cval(1,2); 01324 01325 CCC.val(3,10) = A12.cval(1,3); 01326 01327 CCC.val(4,10) = A12.cval(2,1); 01328 01329 CCC.val(5,10) = A12.cval(2,2); 01330 01331 CCC.val(6,10) = A12.cval(2,3); 01332 01333 CCC.val(7,10) = A12.cval(3,1); 01334 01335 CCC.val(8,10) = A12.cval(3,2); 01336 01337 CCC.val(9,10) = A12.cval(3,3); 01338 01339 01340 01341 CCC.val(10,1) = A21.cval(1,1); 01342 01343 CCC.val(10,2) = A21.cval(1,2); 01344 01345 CCC.val(10,3) = A21.cval(1,3); 01346 01347 CCC.val(10,4) = A21.cval(2,1); 01348 01349 CCC.val(10,5) = A21.cval(2,2); 01350 01351 CCC.val(10,6) = A21.cval(2,3); 01352 01353 CCC.val(10,7) = A21.cval(3,1); 01354 01355 CCC.val(10,8) = A21.cval(3,2); 01356 01357 CCC.val(10,9) = A21.cval(3,3); 01358 01359 01360 01361 CCC.val(10,10) = a22; 01362 01363 } 01364 01365 01366 01367 if ( fdEvolutionT ) { 01368 01369 C99 = A13.BJtensor2BJmatrix_2(); 01370 01371 for (i =1; i <=9; i++) { 01372 01373 for (j =1; j <=9; j++) { 01374 01375 CCC.val(10+i,j) = C99.cval(i,j); 01376 01377 CCC.val(j,10+i) = C99.cval(i,j); 01378 01379 } 01380 01381 } 01382 01383 C99 = A33.BJtensor2BJmatrix_2(); 01384 01385 for (i =1; i <=9; i++) { 01386 01387 for (j =1; j <=9; j++) { 01388 01389 CCC.val(10+i,10+j) = C99.cval(i,j); 01390 01391 } 01392 01393 } 01394 01395 } 01396 01397 01398 01399 if ( fdEvolutionS && fdEvolutionT ) { 01400 01401 CCC.val(10,11) = A23.cval(1,1); 01402 01403 CCC.val(10,12) = A23.cval(1,2); 01404 01405 CCC.val(10,13) = A23.cval(1,3); 01406 01407 CCC.val(10,14) = A23.cval(2,1); 01408 01409 CCC.val(10,15) = A23.cval(2,2); 01410 01411 CCC.val(10,16) = A23.cval(2,3); 01412 01413 CCC.val(10,17) = A23.cval(3,1); 01414 01415 CCC.val(10,18) = A23.cval(3,2); 01416 01417 CCC.val(10,19) = A23.cval(3,3); 01418 01419 01420 01421 CCC.val(11,10) = A32.cval(1,1); 01422 01423 CCC.val(12,10) = A32.cval(1,2); 01424 01425 CCC.val(13,10) = A32.cval(1,3); 01426 01427 CCC.val(14,10) = A32.cval(2,1); 01428 01429 CCC.val(15,10) = A32.cval(2,2); 01430 01431 CCC.val(16,10) = A32.cval(2,3); 01432 01433 CCC.val(17,10) = A32.cval(3,1); 01434 01435 CCC.val(18,10) = A32.cval(3,2); 01436 01437 CCC.val(19,10) = A32.cval(3,3); 01438 01439 } 01440 01441 01442 01443 // Inverse of CCC 01444 01445 //CCC.print("C","\n"); 01446 01447 CCC = CCC.inverse(); 01448 01449 //CCC.print(); 01450 01451 01452 01453 // CCC: Matrix -> Tensor 01454 01455 for (i =1; i <=9; i++) { 01456 01457 for (j =1; j <=9; j++) { 01458 01459 C99.val(i,j) = CCC.cval(i,j); 01460 01461 } 01462 01463 } 01464 01465 01466 01467 A11 = C99.BJmatrix2BJtensor_2(); 01468 01469 01470 01471 A12.val(1,1)=CCC.cval(1,10); 01472 01473 A12.val(1,2)=CCC.cval(2,10); 01474 01475 A12.val(1,3)=CCC.cval(3,10); 01476 01477 A12.val(2,1)=CCC.cval(4,10); 01478 01479 A12.val(2,2)=CCC.cval(5,10); 01480 01481 A12.val(2,3)=CCC.cval(6,10); 01482 01483 A12.val(3,1)=CCC.cval(7,10); 01484 01485 A12.val(3,2)=CCC.cval(8,10); 01486 01487 A12.val(3,3)=CCC.cval(9,10); 01488 01489 01490 01491 A21.val(1,1)=CCC.cval(10,1); 01492 01493 A21.val(1,2)=CCC.cval(10,2); 01494 01495 A21.val(1,3)=CCC.cval(10,3); 01496 01497 A21.val(2,1)=CCC.cval(10,4); 01498 01499 A21.val(2,2)=CCC.cval(10,5); 01500 01501 A21.val(2,3)=CCC.cval(10,6); 01502 01503 A21.val(3,1)=CCC.cval(10,7); 01504 01505 A21.val(3,2)=CCC.cval(10,8); 01506 01507 A21.val(3,3)=CCC.cval(10,9); 01508 01509 01510 01511 a22 = CCC.cval(10,10); 01512 01513 01514 01515 for (i =1; i <=9; i++) { 01516 01517 for (j =1; j <=9; j++) { 01518 01519 C99.val(i,j) = CCC.cval(i,10+j); 01520 01521 } 01522 01523 } 01524 01525 A13 = C99.BJmatrix2BJtensor_2(); 01526 01527 01528 01529 for (i =1; i <=9; i++) { 01530 01531 for (j =1; j <=9; j++) { 01532 01533 C99.val(j,i) = CCC.cval(10+i,j); 01534 01535 } 01536 01537 } 01538 01539 A31 = C99.BJmatrix2BJtensor_2(); 01540 01541 01542 01543 A32.val(1,1)=CCC.cval(11,10); 01544 01545 A32.val(1,2)=CCC.cval(12,10); 01546 01547 A32.val(1,3)=CCC.cval(13,10); 01548 01549 A32.val(2,1)=CCC.cval(14,10); 01550 01551 A32.val(2,2)=CCC.cval(15,10); 01552 01553 A32.val(2,3)=CCC.cval(16,10); 01554 01555 A32.val(3,1)=CCC.cval(17,10); 01556 01557 A32.val(3,2)=CCC.cval(18,10); 01558 01559 A32.val(3,3)=CCC.cval(19,10); 01560 01561 01562 01563 A23.val(1,1)=CCC.cval(10,11); 01564 01565 A23.val(1,2)=CCC.cval(10,12); 01566 01567 A23.val(1,3)=CCC.cval(10,13); 01568 01569 A23.val(2,1)=CCC.cval(10,14); 01570 01571 A23.val(2,2)=CCC.cval(10,15); 01572 01573 A23.val(2,3)=CCC.cval(10,16); 01574 01575 A23.val(3,1)=CCC.cval(10,17); 01576 01577 A23.val(3,2)=CCC.cval(10,18); 01578 01579 A23.val(3,3)=CCC.cval(10,19); 01580 01581 01582 01583 for (i =1; i <=9; i++) { 01584 01585 for (j =1; j <=9; j++) { 01586 01587 C99.val(i,j) = CCC.cval(10+i,10+j); 01588 01589 } 01590 01591 } 01592 01593 A33 = C99.BJmatrix2BJtensor_2(); 01594 01595 01596 01597 01598 01599 // Return d_gamma 01600 01601 tensorTemp6 = dyods("ij")*LM("ijkl"); 01602 01603 tensorTemp6.null_indices(); 01604 01605 if ( fdEvolutionT ) { 01606 01607 tensorTemp7 = dyoda("ij")*Ktensor("ijkl"); 01608 01609 tensorTemp7.null_indices(); 01610 01611 } 01612 01613 // !!!tensorTemp6&7 will be used later 01614 01615 01616 01617 //Begin, evaluate f^{T}*C 01618 01619 tensorTemp1 = tensorTemp6("kl")*A11("klmn"); 01620 01621 tensorTemp1.null_indices(); 01622 01623 01624 01625 tensorTemp2 = A21 *(dyodq *Kscalar); 01626 01627 tensorTemp4 = tensorTemp1 + tensorTemp2; 01628 01629 if ( fdEvolutionT ) { 01630 01631 tensorTemp3 = tensorTemp7("kl")*A31("klmn"); 01632 01633 tensorTemp3.null_indices(); 01634 01635 tensorTemp4 += tensorTemp3; 01636 01637 } 01638 01639 01640 01641 tensorTemp1 = tensorTemp6("kl")*A12("kl"); 01642 01643 tensorTemp1.null_indices(); 01644 01645 temp1 = tensorTemp1.trace(); 01646 01647 temp2 = a22 *(dyodq *Kscalar); 01648 01649 temp2 += temp1; 01650 01651 if ( fdEvolutionT ) { 01652 01653 tensorTemp3 = tensorTemp7("kl")*A32("kl"); 01654 01655 tensorTemp3.null_indices(); 01656 01657 temp3 = tensorTemp3.trace(); 01658 01659 temp2 += temp3; 01660 01661 } 01662 01663 01664 01665 tensorTemp1 = tensorTemp6("kl")*A13("klmn"); 01666 01667 tensorTemp1.null_indices(); 01668 01669 tensorTemp2 = A23 *(dyodq *Kscalar); 01670 01671 tensorTemp5 = tensorTemp1 + tensorTemp2; 01672 01673 if ( fdEvolutionT ) { 01674 01675 tensorTemp3 = tensorTemp7("kl")*A33("klmn"); 01676 01677 tensorTemp3.null_indices(); 01678 01679 tensorTemp5 += tensorTemp3; 01680 01681 } 01682 01683 //End, evaluate f^{T}*C: tensorTemp4-temp2-tensorTemp5 01684 01685 01686 01687 01688 01689 //Begin, evaluate f^{T}*C*r 01690 01691 tensorTemp1 = tensorTemp4("ij") *res_Ee("ij"); 01692 01693 tensorTemp1.null_indices(); 01694 01695 temp1 = tensorTemp1.trace(); 01696 01697 temp4 = temp1 + temp2*res_xi; 01698 01699 if ( fdEvolutionT ) { 01700 01701 tensorTemp3 = tensorTemp5("ij") *res_eta("ij"); 01702 01703 tensorTemp3.null_indices(); 01704 01705 temp3 = tensorTemp3.trace(); 01706 01707 temp4 += temp3; 01708 01709 } 01710 01711 //End, evaluate f^{T}*C*r: temp4 01712 01713 01714 01715 //Begin, evaluate f^{T}*C*m 01716 01717 tensorTemp1 = tensorTemp4("ij") *MCtensor("ij"); 01718 01719 tensorTemp1.null_indices(); 01720 01721 temp1 = tensorTemp1.trace(); 01722 01723 temp5 = temp1 + temp2*nscalar; 01724 01725 if ( fdEvolutionT ) { 01726 01727 tensorTemp3 = tensorTemp5("ij") *ntensor("ij"); 01728 01729 tensorTemp3.null_indices(); 01730 01731 temp3 = tensorTemp3.trace(); 01732 01733 temp5 += temp3; 01734 01735 } 01736 01737 //End, evaluate f^{T}*C*m: temp5 01738 01739 // !!! temp5 will be used later 01740 01741 01742 01743 d_gamma = (yieldfun - temp4) / temp5; // Here is d_gamma!!! 01744 01745 //if(d_gamma < 0.0) d_gamma = 0.0; 01746 01747 01748 01749 //Begin, Calculate incremental variables 01750 01751 tensorTemp2 = MCtensor*(-d_gamma) - res_Ee; 01752 01753 temp2 = nscalar*(-d_gamma) - res_xi; 01754 01755 if ( fdEvolutionT ) { 01756 01757 tensorTemp8 = ntensor*(-d_gamma) - res_eta; 01758 01759 } 01760 01761 01762 01763 tensorTemp1 = A11("ijkl") *tensorTemp2("kl"); 01764 01765 tensorTemp1.null_indices(); 01766 01767 D_Ee = tensorTemp1 + A12*temp2; 01768 01769 if ( fdEvolutionT ) { 01770 01771 tensorTemp3 = A13("ijkl") *tensorTemp8("kl"); 01772 01773 tensorTemp3.null_indices(); 01774 01775 D_Ee += tensorTemp3; 01776 01777 } 01778 01779 01780 01781 tensorTemp1 = A21("kl") *tensorTemp2("kl"); 01782 01783 tensorTemp1.null_indices(); 01784 01785 temp1 = tensorTemp1.trace(); 01786 01787 D_xi = temp1 + a22*temp2 ; 01788 01789 if ( fdEvolutionT ) { 01790 01791 tensorTemp3 = A23("kl") *tensorTemp8("kl"); 01792 01793 tensorTemp3.null_indices(); 01794 01795 temp3 = tensorTemp3.trace(); 01796 01797 D_xi += temp3; 01798 01799 } 01800 01801 01802 01803 if ( fdEvolutionT ) { 01804 01805 tensorTemp1 = A31("ijkl") *tensorTemp2("kl"); 01806 01807 tensorTemp1.null_indices(); 01808 01809 tensorTemp3 = A33("ijkl") *tensorTemp8("kl"); 01810 01811 tensorTemp3.null_indices(); 01812 01813 D_eta = tensorTemp1 + tensorTemp3 + (A32*temp2); 01814 01815 } 01816 01817 //End, Calculate incremental variables: D_Ee, D_xi, D_eta 01818 01819 01820 01821 // Update Variables 01822 01823 D_gamma += d_gamma; // updated D_gamma 01824 01825 01826 01827 Ee += D_Ee; 01828 01829 Ce = Ee*2.0 + tensorI2; 01830 01831 fde3d->setTrialC(Ce); // Note: It is C, not F!!! 01832 01833 B_PK2 = fde3d->getStressTensor(); // Updated B_PK2 01834 01835 B_Mandel = Ce("ik")*B_PK2("kj"); // Update Mandel Stress 01836 01837 B_Mandel.null_indices(); 01838 01839 01840 01841 xi += D_xi; 01842 01843 q += Kscalar*D_xi; 01844 01845 fdeps->setStrainLikeInVar(xi); 01846 01847 fdeps->setStressLikeInVar(q); 01848 01849 01850 01851 if ( fdEvolutionT ) { 01852 01853 eta += D_eta; 01854 01855 tensorTemp2 = Ktensor("ijkl")*D_eta("kl"); 01856 01857 tensorTemp2.null_indices(); 01858 01859 a += tensorTemp2; 01860 01861 fdeps->setStrainLikeKiVar(eta); 01862 01863 fdeps->setStressLikeKiVar(a); 01864 01865 } 01866 01867 01868 01869 //Begin, Calculate residuals 01870 01871 res_Ee = (MCtensor*D_gamma) + Ee - Ee_n; 01872 01873 res_xi = (nscalar*D_gamma) + xi - xi_n; 01874 01875 01876 01877 tensorTemp1 = res_Ee("ij")*res_Ee("ij"); 01878 01879 tensorTemp1.null_indices(); 01880 01881 res_norm_Ee = tensorTemp1.trace(); 01882 01883 residual = sqrt( res_norm_Ee + res_xi*res_xi ); 01884 01885 01886 01887 if ( fdEvolutionT ) { 01888 01889 res_eta = (ntensor*D_gamma) + eta - eta_n; 01890 01891 tensorTemp3 = res_eta("ij")*res_eta("ij"); 01892 01893 tensorTemp3.null_indices(); 01894 01895 res_norm_eta = tensorTemp3.trace(); 01896 01897 residual = sqrt( res_norm_Ee + res_xi*res_xi + res_norm_eta ); 01898 01899 } 01900 01901 //End, Calculate residuals: residual 01902 01903 01904 01905 yieldfun = fdy->Yd(B_Mandel, *fdeps); // Updated yieldfun 01906 01907 //printf("Y= %e\n ", yieldfun); 01908 01909 01910 01911 iter_counter++; 01912 01913 01914 01915 if ( iter_counter > Max_Iter ) { 01916 01917 opserr << "Warnning: Iteration More than " << Max_Iter; 01918 01919 opserr << " in return mapping algorithm of FD EP model" << "\n"; 01920 01921 //exit (-1); 01922 01923 } 01924 01925 01926 01927 } while ( yieldfun > fdy->getTolerance() || residual > tolerance && iter_counter < Max_Iter ); // end of do - while 01928 01929 01930 01931 // For Numerical stability 01932 01933 D_gamma *= (1.0 - tolerance); 01934 01935 if ( D_gamma < 0.0 ) 01936 01937 D_gamma = 0.0; 01938 01939 01940 01941 // Update Fp 01942 01943 tensorTemp2 = Mtensor("ij")*Fp_n("jk"); 01944 01945 tensorTemp2.null_indices(); 01946 01947 Fp = Fp_n + tensorTemp2 *D_gamma; 01948 01949 fdeps->setFpInVar(Fp); 01950 01951 Fpinv = Fp.inverse(); //Using the iterative FP 01952 01953 01954 01955 Fe = F("ij")*Fpinv("jk"); Fe.null_indices(); 01956 01957 Ce = Fe("ki")*Fe("kj"); Ce.null_indices(); 01958 01959 fde3d->setTrialC(Ce); // Note: It is C, not F!!! 01960 01961 B_PK2 = fde3d->getStressTensor(); 01962 01963 01964 01965 // iniPK2 01966 01967 iniPK2 = Fpinv("ip")*Fpinv("jq")*B_PK2("pq"); 01968 01969 iniPK2.null_indices(); 01970 01971 01972 01973 //Begin, evaluate the first term of D^{1} 01974 01975 tensorTemp8 = tensorTemp4 * (1.0/temp5); 01976 01977 tensorTemp9 = Mtensor("kj")*tensorTemp8("mn"); 01978 01979 tensorTemp9.null_indices(); 01980 01981 //End, evaluate the first term of D^{1}: tensorTemp9 01982 01983 01984 01985 //Begin, Evaluate \hat{C}^{11} 01986 01987 tensorTemp6 = A11("klmn") *MCtensor("mn"); 01988 01989 tensorTemp6.null_indices(); 01990 01991 tensorTemp7 = A12 *nscalar; 01992 01993 tensorTemp0 = tensorTemp6 + tensorTemp7; 01994 01995 if ( fdEvolutionT ) { 01996 01997 tensorTemp8 = A13("klmn") *ntensor("mn"); 01998 01999 tensorTemp8.null_indices(); 02000 02001 tensorTemp0 += tensorTemp8; 02002 02003 } 02004 02005 tensorTemp0 = tensorTemp0("ij")*tensorTemp4("kl"); 02006 02007 tensorTemp0.null_indices(); 02008 02009 tensorTemp1 = A11 - tensorTemp0*(1.0/temp5); 02010 02011 //End, Evaluate \hat{C}^{11} 02012 02013 02014 02015 //Here using the results of tensorTemp1 02016 02017 LATStensor = Ltensor("ijkl")*tensorTemp1("klmn"); 02018 02019 LATStensor.null_indices(); 02020 02021 02022 02023 //Begin, Evaluate \hat{C}^{21} 02024 02025 tensorTemp6 = A21("mn") *MCtensor("mn"); 02026 02027 tensorTemp6.null_indices(); 02028 02029 temp1 = tensorTemp6.trace(); 02030 02031 temp4 = temp1 + a22 *nscalar; 02032 02033 if ( fdEvolutionT ) { 02034 02035 tensorTemp8 = A23("mn") *ntensor("mn"); 02036 02037 tensorTemp8.null_indices(); 02038 02039 temp3 = tensorTemp8.trace(); 02040 02041 temp4 += temp3; 02042 02043 } 02044 02045 tensorTemp0 = tensorTemp4 *temp4; 02046 02047 tensorTemp2 = A21 - tensorTemp0*(1.0/temp5); 02048 02049 //End, Evaluate \hat{C}^{21} 02050 02051 02052 02053 //Begin, Evaluate \hat{C}^{31} 02054 02055 if ( fdEvolutionT ) { 02056 02057 tensorTemp6 = A31("klmn") *MCtensor("mn"); 02058 02059 tensorTemp6.null_indices(); 02060 02061 tensorTemp7 = A32 *nscalar; 02062 02063 tensorTemp8 = A33("klmn") *ntensor("mn"); 02064 02065 tensorTemp8.null_indices(); 02066 02067 tensorTemp0 = tensorTemp6 + tensorTemp7 + tensorTemp8; 02068 02069 tensorTemp0 = tensorTemp0("ij")*tensorTemp4("kl"); 02070 02071 tensorTemp0.null_indices(); 02072 02073 tensorTemp3 = A31 - tensorTemp0*(1.0/temp5); 02074 02075 } 02076 02077 //End, Evaluate \hat{C}^{31} 02078 02079 02080 02081 // Begin... 02082 02083 tensorTemp6 = (fdf->d2Fodsds(B_Mandel, *fdeps))("ijkl") *LM("klmn"); 02084 02085 tensorTemp6.null_indices(); 02086 02087 tensorTemp7 = (fdf->d2Fodsdq(B_Mandel, *fdeps)) *Kscalar; 02088 02089 if ( fdEvolutionT ) { 02090 02091 tensorTemp8 = (fdf->d2Fodada(B_Mandel, *fdeps))("ijkl")*Ktensor("klmn"); 02092 02093 tensorTemp8.null_indices(); 02094 02095 } 02096 02097 02098 02099 tensorTemp5 = tensorTemp6("ijkl")*tensorTemp1("klmn"); 02100 02101 tensorTemp5.null_indices(); 02102 02103 tensorTemp1 = tensorTemp7("ij")*tensorTemp2("mn"); 02104 02105 tensorTemp1.null_indices(); 02106 02107 tensorTemp5 += tensorTemp1; 02108 02109 if ( fdEvolutionT ) { 02110 02111 tensorTemp2 = tensorTemp8("ijkl")*tensorTemp3("klmn"); 02112 02113 tensorTemp2.null_indices(); 02114 02115 tensorTemp5 += tensorTemp2; 02116 02117 } 02118 02119 02120 02121 tensorTemp9 += (tensorTemp5 *D_gamma); 02122 02123 02124 02125 tensorTemp9 = Fp_ninv("it") * tensorTemp9("tjmn"); 02126 02127 tensorTemp9.null_indices(); 02128 02129 // ...End 02130 02131 02132 02133 tensorTemp5 = tensorI2("il")*Fpinv("nj")*B_PK2("nk")*tensorTemp9("klpq"); 02134 02135 tensorTemp5.null_indices(); 02136 02137 tensorTemp6 = Fpinv("ni")*tensorI2("jl")*B_PK2("nk")*tensorTemp9("klpq"); 02138 02139 tensorTemp6.null_indices(); 02140 02141 02142 02143 tensorTemp1 = Fpinv("im") * Fpinv("jn"); 02144 02145 tensorTemp1.null_indices(); 02146 02147 02148 02149 tensorTemp7 = tensorTemp1("imjn") * LATStensor("mnrs"); 02150 02151 tensorTemp8.null_indices(); 02152 02153 02154 02155 tensorTemp8 = tensorTemp7 - tensorTemp5 - tensorTemp6; 02156 02157 02158 02159 tensorTemp2 = Fp_ninv("kr") * Fp_ninv("ls"); 02160 02161 tensorTemp2.null_indices(); 02162 02163 02164 02165 iniTangent = tensorTemp8("ijrs") * tensorTemp2("krls"); 02166 02167 iniTangent.null_indices(); 02168 02169 02170 02171 } // end of plastic part 02172 02173 else { // start of elastic part 02174 02175 02176 02177 Fe = F("ij") * Fp_ninv("jk"); 02178 02179 Fe.null_indices(); 02180 02181 Ce = Fe("ki") * Fe("kj"); 02182 02183 Ce.null_indices(); 02184 02185 02186 02187 fde3d->setTrialC(Ce); // Note: It is C, not F!!! 02188 02189 B_PK2 = fde3d->getStressTensor(); 02190 02191 02192 02193 iniPK2 = Fp_ninv("ip")*Fp_ninv("jq")*B_PK2("pq"); 02194 02195 iniPK2.null_indices(); 02196 02197 02198 02199 LATStensor = fde3d->getTangentTensor(); 02200 02201 02202 02203 tensorTemp1 = Fp_ninv("im")*Fp_ninv("jn"); 02204 02205 tensorTemp1.null_indices(); 02206 02207 tensorTemp2 = Fp_ninv("kr")*Fp_ninv("ls"); 02208 02209 tensorTemp2.null_indices(); 02210 02211 tensorTemp3 = tensorTemp1("imjn")*LATStensor("mnrs"); 02212 02213 tensorTemp3.null_indices(); 02214 02215 iniTangent = tensorTemp3("ijrs")*tensorTemp2("krls"); 02216 02217 iniTangent.null_indices(); 02218 02219 } 02220 02221 02222 02223 iniGreen = F("ki")*F("kj"); 02224 02225 iniGreen.null_indices(); 02226 02227 iniGreen = (iniGreen - tensorI2) * 0.5; 02228 02229 02230 02231 //cauchystress = Fe("ip")*B_PK2("pq"); 02232 02233 // cauchystress.null_indices(); 02234 02235 //cauchystress = cauchystress("iq")*Fe("jq"); 02236 02237 // cauchystress.null_indices(); 02238 02239 //cauchystress = cauchystress*(1.0/F.determinant()); 02240 02241 02242 02243 return 0; 02244 02245 } 02246 02247 //*/ 02248 02249 02250 02251 //---------------------------------------------------------------------- 02252 02253 //Mandel Version 02254 02255 int FiniteDeformationEP3D::SemiImplicitAlgorithm() 02256 02257 { 02258 02259 // *** This is the key function *** 02260 02261 02262 02263 02264 02265 // Initializing 02266 02267 double yieldfun = 0.0; // trial value of yield function 02268 02269 double D_gamma = 0.0; // consistency parameter /Delta_{gamma} 02270 02271 double d_gamma = 0.0; // increment of consistency parameter /delta_{gamma} 02272 02273 int iter_counter = 0; 02274 02275 02276 02277 straintensor res_Ee = tensorZ2; // residual of intermediate Ee 02278 02279 straintensor res_eta = tensorZ2; // norm of residual of eta 02280 02281 02282 02283 straintensor Fp = tensorI2;; // Plastic deformation gradient 02284 02285 double xi = 0.0; // strain like internal isotropic variable 02286 02287 double q = 0.0; // stress like internal isotropic variable 02288 02289 straintensor eta; // strain like internal kinematic variable 02290 02291 stresstensor a; // stress like internal kinematic variable 02292 02293 02294 02295 straintensor Fpinv = tensorI2; // inverse of Fp 02296 02297 straintensor Ce = tensorI2; // intermediate C 02298 02299 straintensor Ee = tensorZ2; // intermediate Ee 02300 02301 02302 02303 straintensor Fp_n = tensorI2; // Plastic deformation gradient at time n 02304 02305 straintensor Fp_ninv = tensorI2; // Plastic deformation gradient at time n 02306 02307 straintensor Ee_n = tensorZ2; // Ee at the incremental step n, calculated from Fp_n 02308 02309 double xi_n; // xi at the incremental step n, known 02310 02311 straintensor eta_n = tensorZ2; // eta at the incremental step n, known 02312 02313 02314 02315 stresstensor Mtensor = tensorZ2; // --> dFl/dT 02316 02317 stresstensor MCtensor = tensorZ2; // --> dFl/dS 02318 02319 tensor Ltensor = tensorZ4 ; // Tangent tensor in the intermediate configuration 02320 02321 tensor LATStensor = tensorZ4; // Consistent tangent tensor in the intermediate configuration 02322 02323 02324 02325 double nscalar = 0.0; // --> dFl/d(xi) 02326 02327 double Kscalar = 0.0; // Isotropic hardening modoulus 02328 02329 straintensor ntensor; // --> dFl/d(eta) 02330 02331 tensor Ktensor = tensorZ4; // Kinematic hardening modoulus 02332 02333 02334 02335 stresstensor dyods = tensorZ2 ; // --> dY/d(stress) 02336 02337 double dyodq = 0.0; // --> dY/d(xi) 02338 02339 stresstensor dyoda = tensorZ2 ; // --> dY/d(eta) 02340 02341 02342 02343 straintensor D_Ee = tensorZ2; 02344 02345 double D_xi = 0.0; 02346 02347 straintensor D_eta = tensorZ2; 02348 02349 02350 02351 tensor tensorTemp0; 02352 02353 tensor tensorTemp1; 02354 02355 tensor tensorTemp2; 02356 02357 tensor tensorTemp3; 02358 02359 tensor tensorTemp4; 02360 02361 tensor tensorTemp5; 02362 02363 double temp0 = 0.0; 02364 02365 double lowerPart = 0.0; 02366 02367 02368 02369 tensor LM = tensorZ4; // For Mandel Tangent Stiffness 02370 02371 stresstensor B_Mandel = tensorZ2; // For Mandel stress 02372 02373 02374 02375 tensorTemp1 = tensorI2("ij")*tensorI2("kl"); 02376 02377 tensorTemp1.null_indices(); 02378 02379 tensor tensorI4 = tensorTemp1.transpose0110(); 02380 02381 02382 02383 // Read the previous incremental step history variables 02384 02385 Fp = fdeps->getCommitedFpInVar(); 02386 02387 Fp_n = Fp; 02388 02389 Fp_ninv = Fp_n.inverse(); 02390 02391 Fpinv = Fp.inverse(); 02392 02393 Fe = F("ij")*Fpinv("jk"); Fe.null_indices(); 02394 02395 Ce = Fe("ki")*Fe("kj"); Ce.null_indices(); 02396 02397 Ee = (Ce - tensorI2) * 0.5; 02398 02399 Ee_n = Ee; 02400 02401 02402 02403 if ( fdEvolutionS ) { 02404 02405 xi = fdeps->getCommitedStrainLikeInVar(); 02406 02407 xi_n = xi; 02408 02409 q = fdeps->getCommitedStressLikeInVar(); 02410 02411 } 02412 02413 02414 02415 if ( fdEvolutionT ) { 02416 02417 eta = fdeps->getCommitedStrainLikeKiVar(); 02418 02419 eta_n = eta; 02420 02421 a = fdeps->getCommitedStressLikeKiVar(); 02422 02423 } 02424 02425 02426 02427 // Return stress from finite deformation elastic model 02428 02429 fde3d->setTrialC(Ce); // Note: It is C, not F!!! 02430 02431 B_PK2 = fde3d->getStressTensor(); 02432 02433 02434 02435 B_Mandel = Ce("ik")*B_PK2("kj"); // Mandel Stress 02436 02437 B_Mandel.null_indices(); 02438 02439 02440 02441 // Evaluate the value of yield function 02442 02443 yieldfun = fdy->Yd(B_Mandel, *fdeps); 02444 02445 //printf("\nY0= %e\n", yieldfun); 02446 02447 02448 02449 if ( yieldfun > (fdy->getTolerance()) ) { // start of plastic part 02450 02451 D_gamma = 0.0; 02452 02453 d_gamma = 0.0; 02454 02455 iter_counter = 0; 02456 02457 02458 02459 Mtensor = fdf->dFods(B_Mandel, *fdeps); 02460 02461 //MCtensor = Ce("ik")*Mtensor("kj"); 02462 02463 //MCtensor.null_indices(); 02464 02465 //MCtensor = (MCtensor + MCtensor.transpose11()) * 0.5; 02466 02467 if ( fdEvolutionS) 02468 02469 nscalar = fdf->dFodq(B_Mandel, *fdeps); 02470 02471 if ( fdEvolutionT) 02472 02473 ntensor = fdf->dFoda(B_Mandel, *fdeps); 02474 02475 02476 02477 do { // beginning of do - while 02478 02479 MCtensor = Ce("ik")*Mtensor("kj"); 02480 02481 MCtensor.null_indices(); 02482 02483 MCtensor = (MCtensor + MCtensor.transpose11()) * 0.5; 02484 02485 02486 02487 // Return tangent variables 02488 02489 Ltensor = fde3d->getTangentTensor(); 02490 02491 tensorTemp1 = tensorI2("ij") *B_PK2("mn"); 02492 02493 tensorTemp1.null_indices(); 02494 02495 tensorTemp2 = Ce("ik") *Ltensor("kjmn"); 02496 02497 tensorTemp2.null_indices(); 02498 02499 LM = tensorTemp2 + tensorTemp1.transpose0110() + tensorTemp1.transpose0111(); 02500 02501 dyods = fdy->dYods(B_Mandel, *fdeps); 02502 02503 02504 02505 if ( fdEvolutionS) { 02506 02507 Kscalar = fdEvolutionS->HModulus(B_Mandel, *fdeps); 02508 02509 dyodq = fdy->dYodq(B_Mandel, *fdeps); 02510 02511 } 02512 02513 02514 02515 if ( fdEvolutionT) { 02516 02517 Ktensor = fdEvolutionT->HModulus(B_Mandel, *fdeps); 02518 02519 dyoda = fdy->dYoda(B_Mandel, *fdeps); 02520 02521 } 02522 02523 02524 02525 // Return d_gamma 02526 02527 tensorTemp4 = dyods("ij")*LM("ijkl"); 02528 02529 tensorTemp4.null_indices(); 02530 02531 tensorTemp1 = tensorTemp4("kl")*MCtensor("kl"); 02532 02533 tensorTemp1.null_indices(); 02534 02535 lowerPart = tensorTemp1.trace(); 02536 02537 02538 02539 if ( fdEvolutionS ) 02540 02541 lowerPart += dyodq * (Kscalar*nscalar); 02542 02543 02544 02545 if ( fdEvolutionT ) { 02546 02547 tensorTemp5 = dyoda("ij")*Ktensor("ijkl"); 02548 02549 tensorTemp5.null_indices(); 02550 02551 tensorTemp2 = tensorTemp5("kl")*ntensor("kl"); 02552 02553 tensorTemp2.null_indices(); 02554 02555 temp0 = tensorTemp2.trace(); 02556 02557 lowerPart += temp0; 02558 02559 } 02560 02561 02562 02563 if (lowerPart != 0.0) 02564 02565 d_gamma = yieldfun / lowerPart; 02566 02567 //if(d_gamma < 0.0) d_gamma = 0.0; 02568 02569 //printf("d_gamma= %e\n ", d_gamma); 02570 02571 02572 02573 //Begin, Calculate incremental variables 02574 02575 D_Ee = MCtensor * (-d_gamma); 02576 02577 if ( fdEvolutionS ) 02578 02579 D_xi = nscalar * (-d_gamma); 02580 02581 if ( fdEvolutionT ) 02582 02583 D_eta = ntensor * (-d_gamma); 02584 02585 02586 02587 // Update Variable 02588 02589 D_gamma += d_gamma; // updated D_gamma 02590 02591 02592 02593 Ee += D_Ee; 02594 02595 Ce = Ee*2.0 + tensorI2; 02596 02597 fde3d->setTrialC(Ce); // Note: It is C, not F!!! 02598 02599 B_PK2 = fde3d->getStressTensor(); // Updated B_PK2 02600 02601 B_Mandel = Ce("ik")*B_PK2("kj"); // Update Mandel Stress 02602 02603 B_Mandel.null_indices(); 02604 02605 02606 02607 if ( fdEvolutionS ) { 02608 02609 xi += D_xi; 02610 02611 q += (Kscalar * D_xi); 02612 02613 fdeps->setStrainLikeInVar(xi); 02614 02615 fdeps->setStressLikeInVar(q); 02616 02617 } 02618 02619 02620 02621 if ( fdEvolutionT ) { 02622 02623 eta += D_eta; 02624 02625 tensorTemp1 = Ktensor("ijkl") *D_eta("kl"); 02626 02627 tensorTemp1.null_indices(); 02628 02629 a += tensorTemp1; 02630 02631 fdeps->setStrainLikeKiVar(eta); 02632 02633 fdeps->setStressLikeKiVar(a); 02634 02635 } 02636 02637 02638 02639 yieldfun = fdy->Yd(B_Mandel, *fdeps); // Updated yieldfun 02640 02641 //printf("Y= %e\n ", yieldfun); 02642 02643 02644 02645 iter_counter++; 02646 02647 02648 02649 if ( iter_counter > Max_Iter ) { 02650 02651 opserr << "Warnning: Iteration More than " << Max_Iter; 02652 02653 opserr << " in return mapping algorithm of FD EP model" << "\n"; 02654 02655 //exit (-1); 02656 02657 } 02658 02659 02660 02661 } while ( yieldfun > fdy->getTolerance() && iter_counter < Max_Iter); // end of do - while 02662 02663 02664 02666 02667 //D_gamma *= (1.0 - tolerance); 02668 02669 //if ( D_gamma < 0.0 ) 02670 02671 // D_gamma = 0.0; 02672 02673 02674 02675 // Update Fp 02676 02677 tensorTemp2 = Mtensor("ij")*Fp_n("jk"); 02678 02679 tensorTemp2.null_indices(); 02680 02681 Fp = Fp_n + (tensorTemp2 *D_gamma); 02682 02683 fdeps->setFpInVar(Fp); 02684 02685 02686 02687 // Return iniTangent and iniPK2 02688 02689 Fpinv = Fp.inverse(); // Using the iterative FP 02690 02691 Fe = F("ij")*Fpinv("jk"); 02692 02693 Fe.null_indices(); 02694 02695 Ce = Fe("ki")*Fe("kj"); 02696 02697 Ce.null_indices(); 02698 02699 02700 02701 fde3d->setTrialC(Ce); // Note: It is C, not F!!! 02702 02703 B_PK2 = fde3d->getStressTensor(); 02704 02705 02706 02707 //iniPK2 = B_PK2; 02708 02709 iniPK2 = Fpinv("ip")*Fpinv("jq")*B_PK2("pq"); 02710 02711 iniPK2.null_indices(); 02712 02713 02714 02715 tensorTemp5 = Ltensor("ijkl") * MCtensor("kl"); 02716 02717 tensorTemp5.null_indices(); 02718 02719 02720 02721 tensorTemp3 = tensorTemp5("ij") * tensorTemp4("mn"); 02722 02723 tensorTemp3.null_indices(); 02724 02725 02726 02727 if (lowerPart != 0.0) 02728 02729 LATStensor = Ltensor - ( tensorTemp3 * (1.0/lowerPart) ); 02730 02731 02732 02733 tensorTemp1 = Mtensor("ri")*Fpinv("lj")*Fp_ninv("kr")*B_PK2("kl"); 02734 02735 tensorTemp1.null_indices(); 02736 02737 tensorTemp2 = Fpinv("ki")*Mtensor("rj")*Fp_ninv("lr")*B_PK2("kl"); 02738 02739 tensorTemp2.null_indices(); 02740 02741 tensorTemp3 = tensorTemp1 + tensorTemp2; 02742 02743 tensorTemp5 = tensorTemp3("ij") * tensorTemp4("mn"); 02744 02745 tensorTemp5.null_indices(); 02746 02747 02748 02749 tensorTemp1 = Fpinv("im") * Fpinv("jn"); 02750 02751 tensorTemp1.null_indices(); 02752 02753 tensorTemp3 = tensorTemp1("imjn") * LATStensor("mnrs"); 02754 02755 tensorTemp3.null_indices(); 02756 02757 02758 02759 tensorTemp3 = tensorTemp3 - ( tensorTemp5 * (1.0/lowerPart) ); 02760 02761 02762 02763 tensorTemp2 = Fpinv("kr") * Fpinv("ls"); 02764 02765 tensorTemp2.null_indices(); 02766 02767 iniTangent = tensorTemp3("ijrs") * tensorTemp2("krls"); 02768 02769 iniTangent.null_indices(); 02770 02771 02772 02773 } // end of plastic part 02774 02775 else { // start of elastic part 02776 02777 02778 02779 // Return iniTangent and iniPK2 02780 02781 02782 02783 Fe = F("ij") * Fp_ninv("jk"); 02784 02785 Fe.null_indices(); 02786 02787 Ce = Fe("ki") * Fe("kj"); 02788 02789 Ce.null_indices(); 02790 02791 02792 02793 fde3d->setTrialC(Ce); // Note: It is C, not F!!! 02794 02795 B_PK2 = fde3d->getStressTensor(); 02796 02797 02798 02799 //iniPK2 = B_PK2; 02800 02801 iniPK2 = Fp_ninv("ip")*Fp_ninv("jq")*B_PK2("pq"); 02802 02803 iniPK2.null_indices(); 02804 02805 02806 02807 LATStensor = fde3d->getTangentTensor(); 02808 02809 02810 02811 tensorTemp1 = Fp_ninv("im") * Fp_ninv("jn"); 02812 02813 tensorTemp1.null_indices(); 02814 02815 tensorTemp2 = Fp_ninv("kr") * Fp_ninv("ls"); 02816 02817 tensorTemp2.null_indices(); 02818 02819 tensorTemp3 = tensorTemp1("imjn") * LATStensor("mnrs"); 02820 02821 tensorTemp3.null_indices(); 02822 02823 iniTangent = tensorTemp3("ijrs") * tensorTemp2("krls"); 02824 02825 iniTangent.null_indices(); 02826 02827 02828 02829 } // end of elastic part 02830 02831 02832 02833 iniGreen = F("ki") * F("kj"); 02834 02835 iniGreen.null_indices(); 02836 02837 iniGreen = (iniGreen - tensorI2) * 0.5; 02838 02839 02840 02841 //cauchystress = Fe("ip") * B_PK2("pq"); 02842 02843 // cauchystress.null_indices(); 02844 02845 //cauchystress = cauchystress("iq") * Fe("jq"); 02846 02847 // cauchystress.null_indices(); 02848 02849 //cauchystress = cauchystress * (1.0/F.determinant()); 02850 02851 02852 02853 return 0; 02854 02855 } 02856 02857 02858 02859 02860 02861 #endif 02862
Generated on Mon Oct 23 15:05:14 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents