CompressibleExponentialLaw.cpp
00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "CompressibleExponentialLaw.hpp"
00030
00031 template<unsigned DIM>
00032 CompressibleExponentialLaw<DIM>::CompressibleExponentialLaw()
00033 {
00034 mA = 0.88;
00035
00036 double bff = 18.5;
00037 double bss = 3.58;
00038 double bnn = 3.58;
00039 double bfn = 2.8;
00040 double bfs = 2.8;
00041 double bsn = 2.8;
00042
00043 mCompressibilityParam = 100.0;
00044
00045 mB.resize(DIM);
00046 for (unsigned i=0; i<DIM; i++)
00047 {
00048 mB[i].resize(DIM);
00049 }
00050
00051 mB[0][0] = bff;
00052 mB[0][1] = mB[1][0] = bfs;
00053 mB[1][1] = bss;
00054
00055 if (DIM > 2)
00056 {
00057 mB[2][2] = bnn;
00058 mB[0][2] = mB[2][0] = bfn;
00059 mB[2][1] = mB[1][2] = bsn;
00060 }
00061
00062 for (unsigned M=0; M<DIM; M++)
00063 {
00064 for (unsigned N=0; N<DIM; N++)
00065 {
00066 mIdentity(M,N) = M==N ? 1.0 : 0.0;
00067 }
00068 }
00069 }
00070
00071 template<unsigned DIM>
00072 void CompressibleExponentialLaw<DIM>::ComputeStressAndStressDerivative(c_matrix<double,DIM,DIM>& rC,
00073 c_matrix<double,DIM,DIM>& rInvC,
00074 double pressure ,
00075 c_matrix<double,DIM,DIM>& rT,
00076 FourthOrderTensor<DIM,DIM,DIM,DIM>& rDTdE,
00077 bool computeDTdE)
00078 {
00079 static c_matrix<double,DIM,DIM> C_transformed;
00080 static c_matrix<double,DIM,DIM> invC_transformed;
00081
00082
00083
00084
00085
00086
00087
00088 ComputeTransformedDeformationTensor(rC, rInvC, C_transformed, invC_transformed);
00089
00090
00091
00092 c_matrix<double,DIM,DIM> E = 0.5*(C_transformed - mIdentity);
00093
00094 double QQ = 0;
00095 for (unsigned M=0; M<DIM; M++)
00096 {
00097 for (unsigned N=0; N<DIM; N++)
00098 {
00099 QQ += mB[M][N]*E(M,N)*E(M,N);
00100 }
00101 }
00102
00103 double multiplier = mA*exp(QQ)/2;
00104 rDTdE.Zero();
00105
00106 double J = sqrt(Determinant(rC));
00107
00108 for (unsigned M=0; M<DIM; M++)
00109 {
00110 for (unsigned N=0; N<DIM; N++)
00111 {
00112 rT(M,N) = multiplier*mB[M][N]*E(M,N) + mCompressibilityParam * J*log(J)*invC_transformed(M,N);
00113
00114 if (computeDTdE)
00115 {
00116 for (unsigned P=0; P<DIM; P++)
00117 {
00118 for (unsigned Q=0; Q<DIM; Q++)
00119 {
00120 rDTdE(M,N,P,Q) = multiplier * mB[M][N] * (M==P)*(N==Q)
00121 + 2*multiplier*mB[M][N]*mB[P][Q]*E(M,N)*E(P,Q)
00122 + mCompressibilityParam * (J*log(J) + J) * invC_transformed(M,N) * invC_transformed(P,Q)
00123 - mCompressibilityParam * 2*J*log(J) * invC_transformed(M,P) * invC_transformed(Q,N);
00124 }
00125 }
00126 }
00127 }
00128 }
00129
00130
00131 this->TransformStressAndStressDerivative(rT, rDTdE, computeDTdE);
00132 }
00133
00135
00137
00138 template class CompressibleExponentialLaw<2>;
00139 template class CompressibleExponentialLaw<3>;