swh:1:snp:f521c49ab17ef7db6ec70b2430e1ed203f50383f
Tip revision: 588e391cde1c9be3cb8ae9670e10fff235d3c4ca authored by Dmitry Yu. Naumov on 01 November 2021, 09:31:34 UTC
Merge branch 'UpdateXDMFLibraryForNewerHDF5' into 'master'
Merge branch 'UpdateXDMFLibraryForNewerHDF5' into 'master'
Tip revision: 588e391
LinearElasticIsotropicPhaseField.h
/**
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#pragma once
#include "MathLib/KelvinVector.h"
namespace MaterialLib
{
namespace Solids
{
namespace Phasefield
{
/** Decompose the stiffness into tensile and compressive part.
* Judging by the physical observations, compression perpendicular
* to a crack does not cause crack propagation. Thus,
* the phase-field parameter is only involved into the tensile part
* to degrade the elastic strain energy.
*/
template <int DisplacementDim>
std::tuple<MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_real */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_compressive */,
double /* strain_energy_tensile */, double /* elastic_energy */
>
calculateVolDevDegradedStress(
double const degradation,
double const bulk_modulus,
double const mu,
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps)
{
static int const KelvinVectorSize =
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
using KelvinMatrix =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
// calculation of deviatoric parts
auto const& P_dev = Invariants::deviatoric_projection;
KelvinVector const epsd_curr = P_dev * eps;
// Hydrostatic part for the stress and the tangent.
double const eps_curr_trace = Invariants::trace(eps);
KelvinMatrix C_tensile = KelvinMatrix::Zero();
KelvinMatrix C_compressive = KelvinMatrix::Zero();
if (eps_curr_trace >= 0)
{
double const strain_energy_tensile =
bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
mu * epsd_curr.transpose() * epsd_curr;
KelvinVector const sigma_tensile =
bulk_modulus * eps_curr_trace * Invariants::identity2 +
2 * mu * epsd_curr;
C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
double const elastic_energy = degradation * strain_energy_tensile;
KelvinVector const sigma_real =
degradation * sigma_tensile; // + sigma_compressive, which is zero;
return std::make_tuple(sigma_real, sigma_tensile, C_tensile,
C_compressive, strain_energy_tensile,
elastic_energy);
}
double const strain_energy_tensile = mu * epsd_curr.transpose() * epsd_curr;
KelvinVector const sigma_tensile = 2 * mu * epsd_curr;
KelvinVector const sigma_compressive =
bulk_modulus * eps_curr_trace * Invariants::identity2;
C_tensile.noalias() = 2 * mu * P_dev * KelvinMatrix::Identity();
C_compressive.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
double const elastic_energy =
bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
mu * epsd_curr.transpose() * epsd_curr;
KelvinVector const sigma_real =
degradation * sigma_tensile + sigma_compressive;
return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy);
}
template <int DisplacementDim>
std::tuple<MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_real */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_compressive */,
double /* strain_energy_tensile */, double /* elastic_energy */
>
calculateIsotropicDegradedStress(
double const degradation,
double const bulk_modulus,
double const mu,
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps)
{
static int const KelvinVectorSize =
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
using KelvinMatrix =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
// calculation of deviatoric parts
auto const& P_dev = Invariants::deviatoric_projection;
KelvinVector const epsd_curr = P_dev * eps;
// Hydrostatic part for the stress and the tangent.
double const eps_curr_trace = Invariants::trace(eps);
KelvinMatrix C_tensile = KelvinMatrix::Zero();
KelvinMatrix C_compressive = KelvinMatrix::Zero();
double const strain_energy_tensile =
bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
mu * epsd_curr.transpose() * epsd_curr;
double const elastic_energy = degradation * strain_energy_tensile;
KelvinVector const sigma_tensile =
bulk_modulus * eps_curr_trace * Invariants::identity2 +
2 * mu * epsd_curr;
C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
KelvinVector const sigma_real = degradation * sigma_tensile;
return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy);
}
} // namespace Phasefield
} // namespace Solids
} // namespace MaterialLib