Skip to content
Snippets Groups Projects
Commit 898b8bbb authored by wenqing's avatar wenqing
Browse files

[BGRa] Added some comments

parent 86dbee6c
No related branches found
No related tags found
No related merge requests found
...@@ -48,19 +48,29 @@ CreepBGRa<DisplacementDim>::integrateStress( ...@@ -48,19 +48,29 @@ CreepBGRa<DisplacementDim>::integrateStress(
std::exp(-_q / (MaterialLib::PhysicalConstant::IdealGasConstant * T)); std::exp(-_q / (MaterialLib::PhysicalConstant::IdealGasConstant * T));
auto const& deviatoric_matrix = Invariants::deviatoric_projection; auto const& deviatoric_matrix = Invariants::deviatoric_projection;
// In newton_solver.solve(), the Jacobian is calculated first, and then
// then comes the assembly of the residue vector. In order to save
// computation time, the following two variables are set visible in this
// function. The two variables keep the computation results of the
// norm of the deviatoric stress and 2bG||s_{n+1}||^{n-1} in
// update_jacobian, and then they are directly used in update_residual
// without repeating the computation. Of course, update_jacobian is not a
// pure function anymore, but has side effects.
double pow_norm_s_n1_n_minus_one_2b_G = 0.; double pow_norm_s_n1_n_minus_one_2b_G = 0.;
KelvinVector s_n1; KelvinVector s_n1;
auto const update_jacobian = [&](JacobianMatrix& jacobian) { auto const update_jacobian = [&](JacobianMatrix& jacobian) {
// side effect
s_n1 = deviatoric_matrix * solution; s_n1 = deviatoric_matrix * solution;
double const norm_s_n1 = Invariants::FrobeniusNorm(s_n1); double const norm_s_n1 = Invariants::FrobeniusNorm(s_n1);
// side effect
pow_norm_s_n1_n_minus_one_2b_G = pow_norm_s_n1_n_minus_one_2b_G =
2.0 * b * this->_mp.mu(t, x) * std::pow(norm_s_n1, _n - 1); 2.0 * b * this->_mp.mu(t, x) * std::pow(norm_s_n1, _n - 1);
jacobian = jacobian = KelvinMatrix::Identity() +
KelvinMatrix::Identity() + pow_norm_s_n1_n_minus_one_2b_G *
pow_norm_s_n1_n_minus_one_2b_G * (deviatoric_matrix +
(deviatoric_matrix + ((_n - 1) / (norm_s_n1 * norm_s_n1)) * ((_n - 1) / (norm_s_n1 * norm_s_n1)) * s_n1 *
s_n1 * s_n1.transpose()); s_n1.transpose());
}; };
auto const update_residual = [&](ResidualVectorType& r) { auto const update_residual = [&](ResidualVectorType& r) {
......
...@@ -93,7 +93,7 @@ public: ...@@ -93,7 +93,7 @@ public:
private: private:
NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters; NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters;
const double _n; /// parameter determined by experiments. const double _n; /// Creep rate exponent n.
/// $fA\left(\frac{3}{2}\right)^{n/2+1}/\sigma_{eff}^n $f /// $fA\left(\frac{3}{2}\right)^{n/2+1}/\sigma_{eff}^n $f
const double _coef; const double _coef;
const double _q; /// Activation energy const double _q; /// Activation energy
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment