32 #include "libmesh/tensor_value.h"
36 template <
typename StrainEnergy>
43 template <
typename StrainEnergy>
49 template <
typename StrainEnergy>
58 libMesh::Real lambda_sq = A_cov(2,2);
60 libMesh::Real A_over_a = 1.0/lambda_sq;
63 this->compute_I1_I2(a_contra,a_cov,A_contra,A_cov,lambda_sq,A_over_a,I1,I2);
65 libMesh::Real a_term, A_term;
66 this->compute_stress_terms( lambda_sq, A_over_a, I1, I2, a_term, A_term );
70 for(
unsigned int alpha = 0; alpha < 2; alpha++ )
72 for(
unsigned int beta = 0; beta < 2; beta++ )
74 stress(alpha,beta) = a_contra(alpha,beta)*a_term + A_contra(alpha,beta)*A_term;
81 template <
typename StrainEnergy>
91 libMesh::Real lambda_sq = A_cov(2,2);
93 libMesh::Real A_over_a = 1.0/lambda_sq;
96 this->compute_I1_I2(a_contra,a_cov,A_contra,A_cov,lambda_sq,A_over_a,I1,I2);
98 libMesh::Real a_term, A_term;
99 this->compute_stress_terms( lambda_sq, A_over_a, I1, I2, a_term, A_term );
102 this->compute_stress_deriv_terms( lambda_sq, A_over_a, I1, I2, a_contra, A_contra, daterm_dstrain, dAterm_dstrain );
105 this->compute_Acontra_deriv( A_contra, dAcontra_dstrain );
109 for(
unsigned int alpha = 0; alpha < 2; alpha++ )
111 for(
unsigned int beta = 0; beta < 2; beta++ )
113 stress(alpha,beta) = a_contra(alpha,beta)*a_term + A_contra(alpha,beta)*A_term;
115 for(
unsigned int lambda = 0; lambda < 2; lambda++ )
117 for(
unsigned int mu = 0; mu < 2; mu++ )
119 C(alpha,beta,lambda,mu) = a_contra(alpha,beta)*daterm_dstrain(lambda,mu)
120 + dAcontra_dstrain(alpha,beta,lambda,mu)*A_term
121 + A_contra(alpha,beta)*dAterm_dstrain(lambda,mu);
130 template <
typename StrainEnergy>
135 libMesh::Real lambda_sq, libMesh::Real A_over_a,
136 libMesh::Real& I1, libMesh::Real& I2 )
const
140 for(
unsigned int alpha = 0; alpha < 2; alpha++ )
142 for(
unsigned int beta = 0; beta < 2; beta++ )
144 I1 += a_contra(alpha,beta)*A_cov(alpha,beta);
145 I2 += a_cov(alpha,beta)*A_contra(alpha,beta);
155 template <
typename StrainEnergy>
157 libMesh::Real I1, libMesh::Real I2,
158 libMesh::Real& a_term, libMesh::Real& A_term )
const
160 libMesh::Real dWdI1 = _W.dI1(I1,I2,1.0);
161 libMesh::Real dWdI2 = _W.dI2(I1,I2,1.0);
165 libMesh::Real p = -2.0*lambda_sq*( dWdI1 + dWdI2*(I1-lambda_sq) );
167 a_term = 2.0*(dWdI1 + dWdI2*lambda_sq);
168 A_term = 2.0*dWdI2*A_over_a + p;
173 template <
typename StrainEnergy>
175 libMesh::Real I1, libMesh::Real I2,
181 daterm_dstrain.zero();
182 dAterm_dstrain.zero();
184 libMesh::Real dWdI1 = _W.dI1(I1,I2,1.0);
185 libMesh::Real dWdI2 = _W.dI2(I1,I2,1.0);
190 for(
unsigned int alpha = 0; alpha < 2; alpha++ )
192 for(
unsigned int beta = 0; beta < 2; beta++ )
197 const libMesh::Real dlamsq_dstrain = -2.0*lambda_sq*A_contra(alpha,beta);
199 daterm_dstrain(alpha,beta) = 2.0*dWdI2*dlamsq_dstrain;
204 const libMesh::Real dI1_dstrain = 2.0*a_contra(alpha,beta) + dlamsq_dstrain;
206 const libMesh::Real dp_dstrain = -2.0*dlamsq_dstrain*( dWdI1 + dWdI2*(I1-lambda_sq) )
207 -2.0*lambda_sq*dWdI2*(dI1_dstrain - dlamsq_dstrain);
209 dAterm_dstrain(alpha,beta) = 2.0*dWdI2*A_over_a*(2.0*A_contra(alpha,beta)) + dp_dstrain;
216 template <
typename StrainEnergy>
220 for(
unsigned int alpha = 0; alpha < 2; alpha++ )
222 for(
unsigned int beta = 0; beta < 2; beta++ )
224 for(
unsigned int lambda = 0; lambda < 2; lambda++ )
226 for(
unsigned int mu = 0; mu < 2; mu++ )
228 dAcontra_dstrain(alpha,beta,lambda,mu) = -( A_contra(alpha,lambda)*A_contra(beta,mu)
229 + A_contra(alpha,mu)*A_contra(beta,lambda) );
238 template <
typename StrainEnergy>
244 std::cerr <<
"Error: compute_33_stress shouldn't be called for incompressible materials." << std::endl;
void compute_Acontra_deriv(const libMesh::TensorValue< libMesh::Real > &A_contra, ElasticityTensor &dAcontra_dstrain) const
virtual ~IncompressiblePlaneStressHyperelasticity()
void compute_stress_deriv_terms(libMesh::Real lambda_sq, libMesh::Real A_over_a, libMesh::Real I1, libMesh::Real I2, const libMesh::TensorValue< libMesh::Real > &a_contra, const libMesh::TensorValue< libMesh::Real > &A_contra, libMesh::TensorValue< libMesh::Real > &daterm_dstrain, libMesh::TensorValue< libMesh::Real > &dAterm_dstrain) const
void compute_stress_imp(unsigned int dim, const libMesh::TensorValue< libMesh::Real > &a_contra, const libMesh::TensorValue< libMesh::Real > &a_cov, const libMesh::TensorValue< libMesh::Real > &A_contra, const libMesh::TensorValue< libMesh::Real > &A_cov, libMesh::TensorValue< libMesh::Real > &stress)
void compute_stress_terms(libMesh::Real lambda_sq, libMesh::Real A_over_a, libMesh::Real I1, libMesh::Real I2, libMesh::Real &a_term, libMesh::Real &A_term) const
libMesh::Real compute_33_stress_imp(const libMesh::TensorValue< libMesh::Real > &g_contra, const libMesh::TensorValue< libMesh::Real > &g_cov, const libMesh::TensorValue< libMesh::Real > &G_contra, const libMesh::TensorValue< libMesh::Real > &G_cov)
void compute_I1_I2(const libMesh::TensorValue< libMesh::Real > &a_contra, const libMesh::TensorValue< libMesh::Real > &a_cov, const libMesh::TensorValue< libMesh::Real > &A_contra, const libMesh::TensorValue< libMesh::Real > &A_cov, libMesh::Real lambda_sq, libMesh::Real A_over_a, libMesh::Real &I1, libMesh::Real &I2) const
IncompressiblePlaneStressHyperelasticity(const GetPot &input)
void compute_stress_and_elasticity_imp(unsigned int dim, const libMesh::TensorValue< libMesh::Real > &a_contra, const libMesh::TensorValue< libMesh::Real > &a_cov, const libMesh::TensorValue< libMesh::Real > &A_contra, const libMesh::TensorValue< libMesh::Real > &A_cov, libMesh::TensorValue< libMesh::Real > &stress, ElasticityTensor &C)