Dear All,
I am solving an elastic problem coupled with heat equation transfer
(thermal-elasticity). I should solve the problem in a curved domain (circle
in 2d and sphere in 3d) . I am developing step 44 in this case. I have to
define some symmetric tensors. I have attached some parts of my code that i
have problem with. I have two unknown variables. one for displacement
(vector) and another for other parameter(scalar). I do really appreciate
your kindness if you let me know how I can define symmetric tensor or in
general tensor for this case. In addition, should I define DofHandler<dim,
spacedim> for both scalar and vector valued unknown? My main problem is how
to define tensor in different dim and spacedim problem.
template <int spacedim>
class PointHistory
{
public:
PointHistory()
:
material(NULL),
F_inv(StandardTensors<spacedim>::I),
elasStress(SymmetricTensor<2, spacedim>()),
totalStress(SymmetricTensor<2, spacedim>()),
Jc(SymmetricTensor<4, spacedim>())
{}
virtual ~PointHistory()
{
delete material;
material = NULL;
}
void setup_lqp (const Parameters::AllParameters ¶meters)
{
material = new Material_Variables<spacedim>(parameters.lambda0,
parameters.mu0, parameters.lambda1, parameters.mu1,
parameters.chi, parameters.H,
parameters.thetain, parameters.thetacr, parameters.thetaeq,
parameters.theta, parameters.gammasv);
update_values(Tensor<2, dim>(), double ());
}
void update_values (const Tensor<2, dim> &Grad_u_n, const double eta
)
{
const Tensor<2, dim> eps = 0.5 * (Grad_u_n +
transpose(Grad_u_n)); // Total
Strain
const double phi = std::pow(eta, 2) * (3 - 2 * eta);
// Function
of order parameter
dphi = 6 * eta * (1 - eta);
//
Derivative of the interpolation function
const double doubleWellBarrier = 2*eta + 4 * std::pow(eta, 3) -
6*std::pow(eta, 2);
const double etaBVG = -6 * chi1 * delg1 * (1+eps_V) * dphi;
// Boundary value for
G-L equation
const double etaBVM = 2 * gam1 / 30e-9; // Boundary value for
Mechanics
// Thermal Strain Tensor
const SymmetricTensor<2, dim> epsThermal = alphas * (thetaeq1 -
thetain1) * StandardTensors<dim>::I +
(alpham + (alphas - alpham) * phi) * (theta1 - thetaeq1) *
StandardTensors<dim>::I;
// Transformation Strain Tensor
const SymmetricTensor<2, dim> epsTransformation = (1.0/3.0) * 0.02 *
StandardTensors<dim>::I *
(1 - phi);
//The Initial Strain
const SymmetricTensor<2, dim> epsInitial = SymmetricTensor<2,
dim>(epsThermal) +
(1-phi) * 0.02 * StandardTensors<dim>::I;
// Elastic Strain Tensor
const SymmetricTensor<2, dim> epsElastic = symmetrize(Tensor<2,
dim>(eps)-
SymmetricTensor<2,
dim>(epsTransformation)-
SymmetricTensor<2,
dim>(epsInitial) -
SymmetricTensor<2, dim>(epsThermal));
// The Deviatoric Strain
const SymmetricTensor<2, dim> epsDeviatoric = symmetrize(Tensor<2,
dim>(eps)-(1.0/3.0) * eps_V*
StandardTensors<dim>::I);
material->update_material_data(F, eps, phi, dphi, doubleWellBarrier,
etaBVG, etaBVM
);
F_inv = invert(F);
eps_V = material->get_eps_V();
elasStress = material->get_elasStress();
// surfTenStress = material->get_surfTenStress(); // It is exactly
equal to initial stress...
totalStress = material->get_totalStress(); // Extracting total
Cauchy stress
Jc = material->get_Jc();
// Extracting elasticity stiffnes tensor
driving_force_noStress = material->get_driving_force_noStress();
// Extracting driving force with no stress
eta_boundary_valueG = material->get_eta_boundary_valueG();
// Extracting the boundary value for G-L equation
eta_boundary_valueM = material->get_eta_boundary_valueM();
// Extracting the boundary value for
mechanics
const double temp = -(1.0/(1+eps_V))*0.02*mean_elas_stress;
const double temp1 = temp + 3*(1.0/(1+eps_V))*mean_elas_stress*
(alphas - alpham)*(theta1 - thetaeq1);
const double temp2 = temp1 - ((1.0/(1+eps_V))*(0.5*(lambda1-lambda0)*
SymmetricTensor<2, dim>(epsElastic) *
SymmetricTensor<2, dim>(epsElastic)));
const double temp3 = -(1.0/(1+eps_V))*mu*SymmetricTensor<2,
dim>(epsDeviatoric)*
SymmetricTensor<2, dim>(epsDeviatoric);
// Computing total driving force
get_driv_forc = temp2 * temp3 * dphi + driving_force_noStress;
// computing the total boundary value for the G-L equation
get_eta_bou_valG = eta_boundary_valueG;
// computing the total boundary value for the Mechanics
get_eta_bou_valM = eta_boundary_valueM;
Assert(determinant(F_inv) > 0, ExcInternalError());
}
double get_det_F() const
{
return material->get_det_F();
}
double get_eps_V() const
{
return material->get_eps_V();
}
const Tensor<2, dim> &get_F_inv() const
{
return F_inv;
}
const SymmetricTensor<2, dim> &get_elasStress() const
{
return elasStress;
}
const SymmetricTensor<2, dim> &get_totalStress() const
{
return totalStress;
}
// const SymmetricTensor<2, dim> &get_surfTenStress() const
// {
// return surfTenStress;
// }
const SymmetricTensor<4, dim> &get_Jc() const
{
return Jc;
}
double get_driving_force() const
{
return get_driv_forc;
}
double get_meanElStress() const
{
return mean_elas_stress;
}
// The boundary value for the G-L equation
double get_eta_orderG() const
{
return get_eta_bou_valG;
}
// The boundary value for the Mechanics
double get_eta_orderM() const
{
return get_eta_bou_valM;
}
private:
Material_Variables<dim> *material;
Tensor<2, dim> eps;
Tensor<2, dim> F_inv;
SymmetricTensor<2, dim> epsTransformation;
SymmetricTensor<2, dim> epsThermal;
SymmetricTensor<2, dim> epsInitial;
SymmetricTensor<2, dim> epsElastic;
SymmetricTensor<2, dim> epsDeviatoric;
double eps_V;
Tensor<2, dim> elasStress;
Tensor<2, dim> totalStress;
// SymmetricTensor<2, dim> initialStress;
SymmetricTensor<4, dim> Jc;
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
