Chaste Release::3.1
|
00001 /* 00002 00003 Copyright (c) 2005-2012, University of Oxford. 00004 All rights reserved. 00005 00006 University of Oxford means the Chancellor, Masters and Scholars of the 00007 University of Oxford, having an administrative office at Wellington 00008 Square, Oxford OX1 2JD, UK. 00009 00010 This file is part of Chaste. 00011 00012 Redistribution and use in source and binary forms, with or without 00013 modification, are permitted provided that the following conditions are met: 00014 * Redistributions of source code must retain the above copyright notice, 00015 this list of conditions and the following disclaimer. 00016 * Redistributions in binary form must reproduce the above copyright notice, 00017 this list of conditions and the following disclaimer in the documentation 00018 and/or other materials provided with the distribution. 00019 * Neither the name of the University of Oxford nor the names of its 00020 contributors may be used to endorse or promote products derived from this 00021 software without specific prior written permission. 00022 00023 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 00024 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00025 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 00026 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 00027 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 00028 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE 00029 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 00030 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 00031 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 00032 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00033 00034 */ 00035 00036 #include "SchmidCostaExponentialLaw2d.hpp" 00037 00038 SchmidCostaExponentialLaw2d::SchmidCostaExponentialLaw2d() 00039 { 00040 mA = 0.221; // kiloPascals, presumably, although the paper doesn't say. 00041 // gives results matching Pole-zero anyway. 00042 // Obtained from Table 1 of Schmid reference (see class doxygen), the mu (mean) value. 00043 00044 double bff = 42.5; // dimensionless 00045 double bfs = 11.0; // dimensionless 00046 double bss = 18.6; // dimensionless 00047 00048 mB.resize(2); 00049 mB[0].resize(2); 00050 mB[1].resize(2); 00051 00052 mB[0][0] = bff; 00053 mB[0][1] = bfs; 00054 mB[1][0] = bfs; 00055 mB[1][1] = bss; 00056 00057 for (unsigned M=0; M<2; M++) 00058 { 00059 for (unsigned N=0; N<2; N++) 00060 { 00061 mIdentity(M,N) = M==N ? 1.0 : 0.0; 00062 } 00063 } 00064 } 00065 00066 void SchmidCostaExponentialLaw2d::ComputeStressAndStressDerivative(c_matrix<double,2,2>& rC, 00067 c_matrix<double,2,2>& rInvC, 00068 double pressure, 00069 c_matrix<double,2,2>& rT, 00070 FourthOrderTensor<2,2,2,2>& rDTdE, 00071 bool computeDTdE) 00072 { 00073 static c_matrix<double,2,2> C_transformed; 00074 static c_matrix<double,2,2> invC_transformed; 00075 00076 // The material law parameters are set up assuming the fibre direction is (1,0,0) 00077 // and sheet direction is (0,1,0), so we have to transform C,inv(C),and T. 00078 // Let P be the change-of-basis matrix P = (\mathbf{m}_f, \mathbf{m}_s, \mathbf{m}_n). 00079 // The transformed C for the fibre/sheet basis is C* = P^T C P. 00080 // We then compute T* = T*(C*), and then compute T = P T* P^T. 00081 00082 ComputeTransformedDeformationTensor(rC, rInvC, C_transformed, invC_transformed); 00083 00084 // Compute T* 00085 00086 c_matrix<double,2,2> E = 0.5*(C_transformed - mIdentity); 00087 00088 double Q = 0; 00089 for (unsigned M=0; M<2; M++) 00090 { 00091 for (unsigned N=0; N<2; N++) 00092 { 00093 Q += mB[M][N]*E(M,N)*E(M,N); 00094 } 00095 } 00096 00097 double multiplier = mA*exp(Q)/2; 00098 rDTdE.Zero(); 00099 00100 for (unsigned M=0; M<2; M++) 00101 { 00102 for (unsigned N=0; N<2; N++) 00103 { 00104 rT(M,N) = multiplier*mB[M][N]*E(M,N) - pressure*invC_transformed(M,N); 00105 00106 if (computeDTdE) 00107 { 00108 for (unsigned P=0; P<2; P++) 00109 { 00110 for (unsigned Q=0; Q<2; Q++) 00111 { 00112 rDTdE(M,N,P,Q) = multiplier * mB[M][N] * (M==P)*(N==Q) 00113 + 2*multiplier*mB[M][N]*mB[P][Q]*E(M,N)*E(P,Q) 00114 + 2*pressure*invC_transformed(M,P)*invC_transformed(Q,N); 00115 } 00116 } 00117 } 00118 } 00119 } 00120 00121 // Now do: T = P T* P^T and dTdE_{MNPQ} = P_{Mm}P_{Nn}P_{Pp}P_{Qq} dT*dE*_{mnpq} 00122 this->TransformStressAndStressDerivative(rT, rDTdE, computeDTdE); 00123 } 00124 00125 double SchmidCostaExponentialLaw2d::GetA() 00126 { 00127 return mA; 00128 } 00129 00130 std::vector<std::vector<double> > SchmidCostaExponentialLaw2d::GetB() 00131 { 00132 return mB; 00133 } 00134 00135 double SchmidCostaExponentialLaw2d::GetZeroStrainPressure() 00136 { 00137 return 0.0; 00138 }