This repository has been archived on 2023-11-20. You can view files and clone it, but cannot push or open issues or pull requests.
foam-extend4.1-coherent-io/applications/solvers/solidMechanics/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.C
2013-07-09 11:06:50 +01:00

396 lines
11 KiB
C

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
\*---------------------------------------------------------------------------*/
#include "solidTractionFvPatchVectorField.H"
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"
#include "rheologyModel.H"
#include "plasticityModel.H"
#include "thermalModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
solidTractionFvPatchVectorField::
solidTractionFvPatchVectorField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF
)
:
fixedGradientFvPatchVectorField(p, iF),
fieldName_("undefined"),
traction_(p.size(), vector::zero),
pressure_(p.size(), 0.0),
nonLinear_(OFF)
{
fvPatchVectorField::operator=(patchInternalField());
gradient() = vector::zero;
}
solidTractionFvPatchVectorField::
solidTractionFvPatchVectorField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const dictionary& dict
)
:
fixedGradientFvPatchVectorField(p, iF),
fieldName_(dimensionedInternalField().name()),
traction_("traction", dict, p.size()),
pressure_("pressure", dict, p.size()),
nonLinear_(OFF)
{
fvPatchVectorField::operator=(patchInternalField());
gradient() = vector::zero;
Info << "Patch " << patch().name()
<< "\tTraction boundary field: " << fieldName_ << endl;
//- check if traction boundary is for non linear solver
if (dict.found("nonLinear"))
{
nonLinear_ = nonLinearNames_.read(dict.lookup("nonLinear"));
if (nonLinear_ == UPDATED_LAGRANGIAN)
{
Info << "\tnonLinear set to updated Lagrangian"
<< endl;
}
else if (nonLinear_ == TOTAL_LAGRANGIAN)
{
Info << "\tnonLinear set to total Lagrangian"
<< endl;
}
}
//- the leastSquares has zero non-orthogonal correction
//- on the boundary
//- so the gradient scheme should be extendedLeastSquares
if
(
Foam::word
(
dimensionedInternalField().mesh().schemesDict().gradScheme
(
"grad(" + fieldName_ + ")"
)
) != "extendedLeastSquares"
)
{
Warning << "The gradScheme for " << fieldName_
<< " should be \"extendedLeastSquares 0\" for the boundary "
<< "non-orthogonal correction to be right" << endl;
}
}
solidTractionFvPatchVectorField::
solidTractionFvPatchVectorField
(
const solidTractionFvPatchVectorField& stpvf,
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedGradientFvPatchVectorField(stpvf, p, iF, mapper),
fieldName_(stpvf.fieldName_),
traction_(stpvf.traction_, mapper),
pressure_(stpvf.pressure_, mapper),
nonLinear_(stpvf.nonLinear_)
{}
solidTractionFvPatchVectorField::
solidTractionFvPatchVectorField
(
const solidTractionFvPatchVectorField& stpvf
)
:
fixedGradientFvPatchVectorField(stpvf),
fieldName_(stpvf.fieldName_),
traction_(stpvf.traction_),
pressure_(stpvf.pressure_),
nonLinear_(stpvf.nonLinear_)
{}
solidTractionFvPatchVectorField::
solidTractionFvPatchVectorField
(
const solidTractionFvPatchVectorField& stpvf,
const DimensionedField<vector, volMesh>& iF
)
:
fixedGradientFvPatchVectorField(stpvf, iF),
fieldName_(stpvf.fieldName_),
traction_(stpvf.traction_),
pressure_(stpvf.pressure_),
nonLinear_(stpvf.nonLinear_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void solidTractionFvPatchVectorField::autoMap
(
const fvPatchFieldMapper& m
)
{
fixedGradientFvPatchVectorField::autoMap(m);
traction_.autoMap(m);
pressure_.autoMap(m);
}
// Reverse-map the given fvPatchField onto this fvPatchField
void solidTractionFvPatchVectorField::rmap
(
const fvPatchVectorField& ptf,
const labelList& addr
)
{
fixedGradientFvPatchVectorField::rmap(ptf, addr);
const solidTractionFvPatchVectorField& dmptf =
refCast<const solidTractionFvPatchVectorField>(ptf);
traction_.rmap(dmptf.traction_, addr);
pressure_.rmap(dmptf.pressure_, addr);
}
// Update the coefficients associated with the patch field
void solidTractionFvPatchVectorField::updateCoeffs()
{
if (updated())
{
return;
}
//---------------------------//
//- material properties
//---------------------------//
const rheologyModel& rheology =
this->db().objectRegistry::lookupObject<rheologyModel>("rheologyProperties");
scalarField mu =
rheology.mu()().boundaryField()[patch().index()];
scalarField lambda =
rheology.lambda()().boundaryField()[patch().index()];
if(rheology.type() == plasticityModel::typeName)
{
const plasticityModel& plasticity =
refCast<const plasticityModel>(rheology);
mu = plasticity.newMu().boundaryField()[patch().index()];
lambda = plasticity.newLambda().boundaryField()[patch().index()];
}
//---------------------------//
//- required fields
//---------------------------//
vectorField n = patch().nf();
//- gradient of the field
const fvPatchField<tensor>& gradField =
patch().lookupPatchField<volTensorField, tensor>("grad(" + fieldName_ + ")");
//---------------------------//
//- calculate the traction to apply
//---------------------------//
vectorField Traction(n.size(),vector::zero);
//- total Lagrangian small strain
if(fieldName_ == "U" && nonLinear_ == OFF)
{
//- total traction
Traction = (traction_ - n*pressure_);
}
//- incremental total Lagrangian small strain
else if(fieldName_ == "DU" && nonLinear_ == OFF) //- incremental small strain
{
const fvPatchField<symmTensor>& sigma =
patch().lookupPatchField<volSymmTensorField, symmTensor>("sigma");
//- increment of traction
Traction = (traction_ - n*pressure_) - (n & sigma);
}
//- updated Lagrangian or total Lagrangian large strain
else if (nonLinear_ == UPDATED_LAGRANGIAN || nonLinear_ == TOTAL_LAGRANGIAN)
{
const fvPatchField<symmTensor>& sigma =
patch().lookupPatchField<volSymmTensorField, symmTensor>("sigma");
tensorField F = I + gradField;
tensorField Finv = inv(F);
scalarField J = det(F);
vectorField nCurrent = Finv & n;
nCurrent /= mag(nCurrent);
vectorField tractionCauchy = traction_ - nCurrent*pressure_;
if(nonLinear_ == UPDATED_LAGRANGIAN)
{
//- increment of 2nd Piola-Kirchhoff traction
Traction = (mag(J * Finv & n) * tractionCauchy & Finv) - (n & sigma);
}
else if(nonLinear_ == TOTAL_LAGRANGIAN)
{
//- total 2nd Piola-Kirchhoff traction
Traction = mag(J * Finv & n) * tractionCauchy & Finv;
}
}
else
{
FatalError << "Field " << fieldName_ << " and " << nonLinear_ << " nonLinear are not compatible!"
<< exit(FatalError);
}
//---------------------------//
//- calculate the normal gradient based on the traction
//---------------------------//
vectorField newGradient =
Traction
- (n & (mu*gradField.T() - (mu + lambda)*gradField))
- n*lambda*tr(gradField);
//- if there is plasticity
if(rheology.type() == plasticityModel::typeName)
{
const plasticityModel& plasticity =
refCast<const plasticityModel>(rheology);
newGradient +=
2*mu*(n & plasticity.DEpsilonP().boundaryField()[patch().index()]);
}
//- if there are thermal effects
if(this->db().objectRegistry::foundObject<thermalModel>("thermalProperties"))
{
const thermalModel& thermo =
this->db().objectRegistry::lookupObject<thermalModel>("thermalProperties");
const fvPatchField<scalar>& T =
patch().lookupPatchField<volScalarField, scalar>("T");
const fvPatchField<scalar>& threeKalpha =
patch().lookupPatchField<volScalarField, scalar>("((threeK*rho)*alpha)");
const scalarField T0 = thermo.T0()().boundaryField()[patch().index()];
newGradient += (n*threeKalpha*(T - T0));
}
//- higher order non-linear terms
if(nonLinear_ == UPDATED_LAGRANGIAN || nonLinear_ == TOTAL_LAGRANGIAN)
{
newGradient -=
(n & (mu*(gradField & gradField.T())))
+ 0.5*n*lambda*(gradField && gradField);
//- tensorial identity
//- tr(gradField & gradField.T())*I == (gradField && gradField)*I
}
newGradient /= (2.0*mu + lambda);
gradient() = newGradient;
fixedGradientFvPatchVectorField::updateCoeffs();
}
void solidTractionFvPatchVectorField::evaluate(const Pstream::commsTypes)
{
if (!this->updated())
{
this->updateCoeffs();
}
const fvPatchField<tensor>& gradField =
patch().lookupPatchField<volTensorField, tensor>
(
"grad(" + fieldName_ + ")"
);
vectorField n = patch().nf();
vectorField delta = patch().delta();
//- non-orthogonal correction vectors
vectorField k = delta - n*(n&delta);
Field<vector>::operator=
(
this->patchInternalField()
+ (k&gradField.patchInternalField())
+ gradient()/this->patch().deltaCoeffs()
);
fvPatchField<vector>::evaluate();
}
// Write
void solidTractionFvPatchVectorField::write(Ostream& os) const
{
fvPatchVectorField::write(os);
os.writeKeyword("nonLinear") << nonLinearNames_[nonLinear_] << token::END_STATEMENT << nl;
traction_.writeEntry("traction", os);
pressure_.writeEntry("pressure", os);
writeEntry("value", os);
}
template<>
const char* Foam::NamedEnum<Foam::solidTractionFvPatchVectorField::nonLinearType, 3>::names[] =
{
"off",
"updatedLagrangian",
"totalLagrangian"
};
const Foam::NamedEnum<Foam::solidTractionFvPatchVectorField::nonLinearType, 3>
Foam::solidTractionFvPatchVectorField::nonLinearNames_;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField(fvPatchVectorField, solidTractionFvPatchVectorField);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //