exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Solid models update - boundary traction

Browse source
This commit is contained in:
Hrvoje Jasak 2014-02-23 10:53:50 +00:00
parent 47c824e967
commit 3810b29f37
20 changed files with 1491 additions and 1249 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

986
src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C
View file

File diff suppressed because it is too large Load diff

25
src/solidModels/constitutiveModel/constitutiveModel.H
View file

@ -81,8 +81,8 @@ class constitutiveModel
//- Plane stress
Switch planeStress_;
//- Run-time selectable solidInterface method for correct discretisation
// on bi-material interfaces
//- Run-time selectable solidInterface method for correct
// discretisation on bi-material interfaces
autoPtr<solidInterface> solidInterfacePtr_;
// we use IOReferencer to allow lookup of solidInterface object in
@ -137,17 +137,18 @@ public:
//- If plasticity is active
bool plasticityActive() const
{
if (plasticityStressReturnPtr_.valid())
{
return plasticityStressReturnPtr_->plasticityActive();
}
return false;
if (plasticityStressReturnPtr_.valid())
{
return plasticityStressReturnPtr_->plasticityActive();
}
return false;
}
//- If visco effects are active
bool viscoActive() const
{
return rheologyLawPtr_->viscoActive();
return rheologyLawPtr_->viscoActive();
}
//- Return reference to stress field
@ -241,14 +242,14 @@ public:
//- Return reference to solidInterface
virtual solidInterface& solInterface()
{
return solidInterfacePtr_();
}
return solidInterfacePtr_();
}
//- Return true if solidInterface is active
virtual bool solidInterfaceActive()
{
return solidInterfaceActive_;
}
return solidInterfaceActive_;
}
//- Update yield stress
void updateYieldStress();

885
src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.C
View file

@ -36,438 +36,587 @@ Author
#include "constitutiveModel.H"
#include "thermalModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(tractionBoundaryGradient, 0);
// * * * * * * * * * * * * * * Member functions * * * * * * * * * * * * * * //
tmp<vectorField> tractionBoundaryGradient::traction
Foam::tmp<Foam::vectorField> Foam::tractionBoundaryGradient::traction
(
const tensorField& gradField,
const word fieldName,
const word& workingFieldName,
const word& integralFieldName,
const fvPatch& patch,
Switch orthotropic,
word nonLinear
) const
const bool orthotropic,
const nonLinearGeometry::nonLinearType& nonLinear,
const bool incremental
)
{
// create tmp
// Create result
tmp<vectorField> ttraction
(
new vectorField(gradField.size(), vector::zero)
);
vectorField& traction = ttraction();
// Orthotropic material
if (orthotropic)
{
// for now, turn off orthotropic
FatalError
<< "tractionBoundaryGradient::traction is not written for"
<< " orthotropic yet" << nl
// For now, turn off orthotropic
FatalErrorIn
(
"tmp<vectorField> tractionBoundaryGradient::traction\n"
"(\n"
" const tensorField& gradField,\n"
" const word& workingFieldName,\n"
" const word& integralFieldName,\n"
" const fvPatch& patch,\n"
" const bool orthotropic,\n"
" const nonLinearGeometry::nonLinearType& nonLinear\n"
") const"
) << "tractionBoundaryGradient::traction is not written for"
<< " orthotropic materials yet" << nl
<< "you are probably trying to use solidContact boundaries "
<< "with orthotropic solver" << nl
<< "it should be straight-forward but I have not done it yet!"
<< "with the orthotropic solver" << nl
<< exit(FatalError);
}
else
{
// lookup material properties from the solver
const fvPatchField<scalar>& mu =
patch.lookupPatchField<volScalarField, scalar>("mu");
const fvPatchField<scalar>& lambda =
patch.lookupPatchField<volScalarField, scalar>("lambda");
{
// Lookup material properties from the solver
const fvPatchScalarField& mu =
patch.lookupPatchField<volScalarField, scalar>("mu");
// only for nonlinear elastic properties
// if (rheology.type() == plasticityModel::typeName)
// {
// const plasticityModel& plasticity =
// refCast<const plasticityModel>(rheology);
const fvPatchScalarField& lambda =
patch.lookupPatchField<volScalarField, scalar>("lambda");
// mu = plasticity.newMu().boundaryField()[patch.index()];
// lambda = plasticity.newLambda().boundaryField()[patch.index()];
// }
// only for nonlinear elastic properties
// if (rheology.type() == plasticityModel::typeName)
// {
// const plasticityModel& plasticity =
// refCast<const plasticityModel>(rheology);
// required fields
vectorField n = patch.nf();
// mu = plasticity.newMu().boundaryField()[patch.index()];
// lambda = plasticity.newLambda().boundaryField()[patch.index()];
// }
// Get patch normal
vectorField n = patch.nf();
// Calculate traction
traction = 2*mu*(n & symm(gradField)) + lambda*tr(gradField)*n;
// Plasticity effects
const constitutiveModel& rheology =
patch.boundaryMesh().mesh().objectRegistry::
lookupObject<constitutiveModel>("rheologyProperties");
if (rheology.plasticityActive())
{
traction -=
2*mu*(n & rheology.DEpsilonP().boundaryField()[patch.index()]);
}
// Thermal effects
if
(
patch.boundaryMesh().mesh().objectRegistry::
foundObject<thermalModel>("thermalProperties")
)
{
const thermalModel& thermo =
patch.boundaryMesh().mesh().objectRegistry::
lookupObject<thermalModel>("thermalProperties");
const fvPatchScalarField& threeKalpha =
patch.lookupPatchField<volScalarField, scalar>
("((threeK*rho)*alpha)");
// Incremental thermal contribution
if (incremental)
{
const fvPatchScalarField& DT =
patch.lookupPatchField<volScalarField, scalar>("DT");
traction -= (n*threeKalpha*DT);
}
else
{
const fvPatchScalarField& T =
patch.lookupPatchField<volScalarField, scalar>("T");
const scalarField T0 =
thermo.T0()().boundaryField()[patch.index()];
traction -= (n*threeKalpha*(T - T0));
}
}
// Non-linear terms
if
(
nonLinear == nonLinearGeometry::UPDATED_LAGRANGIAN
|| nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
)
{
traction +=
(n & (mu*(gradField & gradField.T())))
+ 0.5*n*lambda*(gradField && gradField);
if
(
incremental
&& nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
)
{
// Incremental total Lagrangian
const fvPatchTensorField& gradU =
patch.lookupPatchField<volTensorField, tensor>
(
"grad(" + integralFieldName + ")"
);
traction +=
(
n
& (
mu*
(
(gradU & gradField.T())
+ (gradField & gradU.T())
)
)
)
+ 0.5*n*lambda*tr
(
(gradU & gradField.T())
+ (gradField & gradU.T())
);
}
}
// calculate traction
traction = 2*mu*(n&symm(gradField)) + lambda*tr(gradField)*n;
// Add old stress for incremental approaches
if (incremental)
{
// add on old traction
const fvPatchSymmTensorField& sigma =
patch.lookupPatchField<volSymmTensorField, symmTensor>
(
"sigma"
);
traction += (n & sigma);
}
//- if there is plasticity
const constitutiveModel& rheology =
patch.boundaryMesh().mesh().objectRegistry::
lookupObject<constitutiveModel>("rheologyProperties");
if (rheology.plasticityActive())
{
traction -=
2*mu*(n & rheology.DEpsilonP().boundaryField()[patch.index()]);
}
// Visco-elastic effects
if (rheology.viscoActive())
{
const fvPatchSymmTensorField& DSigmaCorr =
patch.lookupPatchField<volSymmTensorField, symmTensor>
(
"DSigmaCorr"
);
//- if there are thermal effects
if (patch.boundaryMesh().mesh().objectRegistry::
foundObject<thermalModel>("thermalProperties"))
{
const thermalModel& thermo =
patch.boundaryMesh().mesh().objectRegistry::
lookupObject<thermalModel>("thermalProperties");
traction += (n & DSigmaCorr);
}
const fvPatchField<scalar>& threeKalpha =
patch.lookupPatchField<volScalarField, scalar>
("((threeK*rho)*alpha)");
// Updated Lagrangian or total Lagrangian large strain
if
(
nonLinear == nonLinearGeometry::UPDATED_LAGRANGIAN
|| nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
)
{
tensorField F = I + gradField;
if (fieldName == "DU")
{
const fvPatchField<scalar>& DT =
patch.lookupPatchField<volScalarField, scalar>("DT");
if
(
incremental
&& nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
)
{
// Incremental total Lagrangian
traction -= (n*threeKalpha*DT);
}
else
{
const fvPatchField<scalar>& T =
patch.lookupPatchField<volScalarField, scalar>("T");
const fvPatchTensorField& gradU =
patch.lookupPatchField<volTensorField, tensor>
(
"grad(" + integralFieldName + ")"
);
const scalarField T0 = thermo.T0()().boundaryField()[patch.index()];
F += gradU;
}
traction -= (n*threeKalpha*(T - T0));
}
}
tensorField Finv = inv(F);
scalarField J = det(F);
// non linear terms
if (nonLinear == "updatedLagrangian" || nonLinear == "totalLagrangian")
{
traction +=
(n & (mu*(gradField & gradField.T())))
+ 0.5*n*lambda*(gradField && gradField);
if (fieldName == "DU" && nonLinear == "totalLagrangian")
{
// incremental total Lagrangian
const fvPatchField<tensor>& gradU =
patch.lookupPatchField<volTensorField, tensor>("grad(U)");
traction +=
(n & (mu*( (gradU & gradField.T()) + (gradField & gradU.T()))))
+ 0.5*n*lambda*tr((gradU & gradField.T()) + (gradField & gradU.T()));
}
}
//- add old stress for incremental approaches
if (fieldName == "DU")
{
// add on old traction
const fvPatchField<symmTensor>& sigma =
patch.lookupPatchField<volSymmTensorField, symmTensor>("sigma");
traction += (n & sigma);
}
//- visco-elastic
if (rheology.viscoActive())
{
const fvPatchField<symmTensor>& DSigmaCorr =
patch.lookupPatchField<volSymmTensorField, symmTensor>("DSigmaCorr");
traction += (n & DSigmaCorr);
}
//- updated Lagrangian or total Lagrangian large strain
if (nonLinear == "updatedLagrangian" || nonLinear == "totalLagrangian")
{
tensorField F = I + gradField;
if (fieldName == "DU" && nonLinear == "totalLagrangian")
{
const fvPatchField<tensor>& gradU =
patch.lookupPatchField<volTensorField, tensor>("grad(U)");
F += gradU;
}
tensorField Finv = inv(F);
scalarField J = det(F);
//- rotate second Piola Kirchhoff traction to be Cauchy traction
traction = (traction & F)/(mag(J * Finv & n));
}
}
// Rotate second Piola Kirchhoff traction to be Cauchy traction
traction = (traction & F)/(mag(J * Finv & n));
}
}
return ttraction;
}
}
// * * * * * * * * * * * * * * * * Operators * * * * * * * * * * * * * * //
tmp<vectorField> tractionBoundaryGradient::operator()
(
const vectorField& traction,
const scalarField& pressure,
const word fieldName,
const fvPatch& patch,
Switch orthotropic,
word nonLinear
) const
{
// create tmp
Foam::tmp<Foam::vectorField> Foam::tractionBoundaryGradient::snGrad
(
const vectorField& traction,
const scalarField& pressure,
const word& workingFieldName,
const word& integralFieldName,
const fvPatch& patch,
const bool orthotropic,
const nonLinearGeometry::nonLinearType& nonLinear,
const bool incremental
)
{
// Create result
tmp<vectorField> tgradient(new vectorField(traction.size(), vector::zero));
vectorField& gradient = tgradient();
// lookup switches
// orthotropic solvers
// Orthotropic material
if (orthotropic)
{
// mechanical properties
const constitutiveModel& rheology =
patch.boundaryMesh().mesh().objectRegistry::
lookupObject<constitutiveModel>("rheologyProperties");
const diagTensorField K = rheology.K()().boundaryField()[patch.index()];
const symmTensor4thOrderField C =
rheology.C()().boundaryField()[patch.index()];
{
// Get mechanical properties
const constitutiveModel& rheology =
patch.boundaryMesh().mesh().objectRegistry::
lookupObject<constitutiveModel>("rheologyProperties");
// required fields
vectorField n = patch.nf();
const diagTensorField Kinv = inv(K);
const fvPatchField<tensor>& gradField =
patch.lookupPatchField<volTensorField, tensor>("grad(" + fieldName + ")");
const diagTensorField K =
rheology.K()().boundaryField()[patch.index()];
const symmTensor4thOrderField C =
rheology.C()().boundaryField()[patch.index()];
// calculate the traction to apply
vectorField Traction(n.size(), vector::zero);
tensorField sigmaExp(n.size(), tensor::zero);
// Required fields
vectorField n = patch.nf();
const diagTensorField Kinv = inv(K);
const fvPatchTensorField& gradField =
patch.lookupPatchField<volTensorField, tensor>
(
"grad(" + workingFieldName + ")"
);
//- total Lagrangian small strain
if (fieldName == "U" && nonLinear == "off")
{
//- total traction
Traction = (traction - n*pressure);
// Calculate the traction to apply
vectorField Traction(n.size(), vector::zero);
tensorField sigmaExp(n.size(), tensor::zero);
sigmaExp = (n*(n&(C && symm(gradField)))) - (K & gradField);
}
//- incremental total Lagrangian small strain
else if (fieldName == "DU" && nonLinear == "off")
{
const fvPatchField<symmTensor>& sigma =
patch.lookupPatchField<volSymmTensorField, symmTensor>("sigma");
// Total Lagrangian, small strain
if
(
!incremental
&& nonLinear == nonLinearGeometry::OFF
)
{
// Use total traction
Traction = (traction - n*pressure);
//- increment of traction
Traction = (traction - n*pressure) - (n & sigma);
sigmaExp = (n*(n&(C && symm(gradField)))) - (K & gradField);
}
// Incremental total Lagrangian small strain
else if
(
incremental
&& nonLinear == nonLinearGeometry::OFF
)
{
const fvPatchSymmTensorField& sigma =
patch.lookupPatchField<volSymmTensorField, symmTensor>
(
"sigma"
);
sigmaExp = (n*(n&(C && symm(gradField)))) - (K & gradField);
}
//- updated Lagrangian large strain
else if (nonLinear == "updatedLagrangian")
{
const fvPatchField<symmTensor>& sigma =
patch.lookupPatchField<volSymmTensorField, symmTensor>("sigma");
// Increment of traction
Traction = (traction - n*pressure) - (n & 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;
sigmaExp = (n*(n&(C && symm(gradField)))) - (K & gradField);
}
// Updated Lagrangian large strain
else if
(
nonLinear == nonLinearGeometry::UPDATED_LAGRANGIAN
)
{
const fvPatchSymmTensorField& sigma =
patch.lookupPatchField<volSymmTensorField, symmTensor>
(
"sigma"
);
//- increment of 2nd Piola-Kirchhoff traction
Traction = (mag(J * Finv & n) * tractionCauchy & Finv) - (n & sigma);
tensorField F = I + gradField;
tensorField Finv = inv(F);
scalarField J = det(F);
vectorField nCurrent = Finv & n;
nCurrent /= mag(nCurrent) + SMALL;
vectorField tractionCauchy = traction - nCurrent*pressure;
sigmaExp = (n*(n&(C && symm(gradField)))) - (K & gradField);
}
else
{
FatalError
<< "solidTractionOrtho boundary condition not suitable for "
<< " field " << fieldName << " and " << nonLinear
<< abort(FatalError);
}
// Increment of 2nd Piola-Kirchhoff traction
Traction = mag(J*(Finv & n))*(tractionCauchy & Finv) - (n & sigma);
gradient =
n &
(
Kinv & ( n*(Traction) - sigmaExp )
);
}
// standard isotropic solvers
else
{
// lookup material properties from the solver
const fvPatchField<scalar>& mu =
patch.lookupPatchField<volScalarField, scalar>("mu");
const fvPatchField<scalar>& lambda =
patch.lookupPatchField<volScalarField, scalar>("lambda");
// only for nonlinear elastic properties
// 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 actual 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")
{
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 == "updatedLagrangian" || nonLinear == "totalLagrangian")
{
const fvPatchField<symmTensor>& sigma =
patch.lookupPatchField<volSymmTensorField, symmTensor>("sigma");
tensorField F = I + gradField;
if (fieldName == "DU" && nonLinear == "totalLagrangian") // for incr TL
{
const fvPatchField<tensor>& gradU =
patch.lookupPatchField<volTensorField, tensor>("grad(U)");
F += gradU;
}
tensorField Finv = hinv(F);
scalarField J = det(F);
vectorField nCurrent = Finv & n;
nCurrent /= mag(nCurrent);
vectorField tractionCauchy = traction - nCurrent*pressure;
Traction = (mag(J * Finv & n) * tractionCauchy & Finv);
if (fieldName == "DU")
{
//- increment of 2nd Piola-Kirchhoff traction
Traction -= (n & sigma);
}
}
else
{
FatalError
<< "Field " << fieldName << " and " << nonLinear
<< " nonLinear are not compatible!"
<< exit(FatalError);
}
//- visco-elastic
const constitutiveModel& rheology =
patch.boundaryMesh().mesh().objectRegistry::
lookupObject<constitutiveModel>("rheologyProperties");
if (rheology.viscoActive())
{
//Info << "visco active" << endl;
const fvPatchField<symmTensor>& DSigmaCorr =
patch.lookupPatchField<volSymmTensorField, symmTensor>("DSigmaCorr");
Traction -= (n & DSigmaCorr);
}
//---------------------------//
//- calculate the normal gradient based on the traction
//---------------------------//
gradient =
Traction
- (n & (mu*gradField.T() - (mu + lambda)*gradField))
- n*lambda*tr(gradField);
//- if there is plasticity
if (rheology.plasticityActive())
{
//Info << "plasticity is active" << endl;
gradient +=
2*mu*(n & rheology.DEpsilonP().boundaryField()[patch.index()]);
}
//- if there are thermal effects
if (patch.boundaryMesh().mesh().objectRegistry::
foundObject<thermalModel>("thermalProperties"))
{
const thermalModel& thermo =
patch.boundaryMesh().mesh().objectRegistry::
lookupObject<thermalModel>("thermalProperties");
const fvPatchField<scalar>& threeKalpha =
patch.lookupPatchField<volScalarField, scalar>
("((threeK*rho)*alpha)");
if (fieldName == "DU")
{
const fvPatchField<scalar>& DT =
patch.lookupPatchField<volScalarField, scalar>("DT");
gradient += (n*threeKalpha*DT);
}
sigmaExp = n*(n &(C && symm(gradField))) - (K & gradField);
}
else
{
const fvPatchField<scalar>& T =
patch.lookupPatchField<volScalarField, scalar>("T");
{
FatalErrorIn
(
"tmp<vectorField> tractionBoundaryGradient::snGrad\n"
"(\n"
" const vectorField& traction,\n"
" const scalarField& pressure,\n"
" const word& workingFieldName,\n"
" const word& integralFieldName,\n"
" const fvPatch& patch,\n"
" const bool orthotropic,\n"
" const nonLinearGeometry::nonLinearType& nonLinear,\n"
" const bool incremental\n"
") const"
) << "solidTractionOrtho boundary condition not suitable for "
<< " field " << workingFieldName << " and "
<< nonLinearGeometry::nonLinearNames_[nonLinear]
<< abort(FatalError);
}
const scalarField T0 = thermo.T0()().boundaryField()[patch.index()];
gradient = n & (Kinv & ( n*(Traction) - sigmaExp ));
}
else
{
// Standard isotropic solvers
gradient += (n*threeKalpha*(T - T0));
}
// Lookup material properties from the solver
const fvPatchScalarField& mu =
patch.lookupPatchField<volScalarField, scalar>("mu");
const fvPatchScalarField& lambda =
patch.lookupPatchField<volScalarField, scalar>("lambda");
// only for nonlinear elastic properties
// if (rheology.type() == plasticityModel::typeName)
// {
// const plasticityModel& plasticity =
// refCast<const plasticityModel>(rheology);
// mu = plasticity.newMu().boundaryField()[patch.index()];
// lambda = plasticity.newLambda().boundaryField()[patch.index()];
// }
vectorField n = patch.nf();
// gradient of the field
const fvPatchTensorField& gradField =
patch.lookupPatchField<volTensorField, tensor>
(
"grad(" + workingFieldName + ")"
);
vectorField Traction(n.size(), vector::zero);
// Total Lagrangian, small strain
if
(
!incremental
&& nonLinear == nonLinearGeometry::OFF
)
{
// Use total traction
Traction = (traction - n*pressure);
}
// Incremental total Lagrangian small strain
else if
(
incremental
&& nonLinear == nonLinearGeometry::OFF
)
{
const fvPatchSymmTensorField& sigma =
patch.lookupPatchField<volSymmTensorField, symmTensor>
(
"sigma"
);
// Increment of traction
Traction = (traction - n*pressure) - (n & sigma);
}
else if
(
nonLinear == nonLinearGeometry::UPDATED_LAGRANGIAN
|| nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
)
{
const fvPatchSymmTensorField& sigma =
patch.lookupPatchField<volSymmTensorField, symmTensor>
(
"sigma"
);
tensorField F = I + gradField;
if
(
incremental
&& nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
)
{
// Incremental total Lagrangian
const fvPatchTensorField& gradU =
patch.lookupPatchField<volTensorField, tensor>
(
"grad(" + integralFieldName + ")"
);
F += gradU;
}
tensorField Finv = hinv(F);
scalarField J = det(F);
vectorField nCurrent = Finv & n;
nCurrent /= mag(nCurrent) + SMALL;
vectorField tractionCauchy = traction - nCurrent*pressure;
Traction = mag(J*(Finv & n))*(tractionCauchy & Finv);
if (incremental)
{
// Increment of 2nd Piola-Kirchhoff traction
Traction -= (n & sigma);
}
}
else
{
FatalErrorIn
(
"tmp<vectorField> tractionBoundaryGradient::snGrad\n"
"(\n"
" const vectorField& traction,\n"
" const scalarField& pressure,\n"
" const word& workingFieldName,\n"
" const word& integralFieldName,\n"
" const fvPatch& patch,\n"
" const bool orthotropic,\n"
" const nonLinearGeometry::nonLinearType& nonLinear,\n"
" const bool incremental\n"
") const"
) << " field " << workingFieldName << " and "
<< nonLinearGeometry::nonLinearNames_[nonLinear]
<< " nonLinear are not compatible!"
<< abort(FatalError);
}
//- higher order non-linear terms
if (nonLinear == "updatedLagrangian" || nonLinear == "totalLagrangian")
{
// no extra relaxation
gradient -=
(n & (mu*(gradField & gradField.T())))
// + 0.5*n*lambda*(gradField && gradField);
//- visco-elastic
const constitutiveModel& rheology =
patch.boundaryMesh().mesh().objectRegistry::
lookupObject<constitutiveModel>("rheologyProperties");
if (rheology.viscoActive())
{
const fvPatchSymmTensorField& DSigmaCorr =
patch.lookupPatchField<volSymmTensorField, symmTensor>
(
"DSigmaCorr"
);
Traction -= (n & DSigmaCorr);
}
// Calculate the normal gradient based on the traction
gradient =
Traction
- (n & (mu*gradField.T() - (mu + lambda)*gradField))
- n*lambda*tr(gradField);
//- Plasticity contribution
if (rheology.plasticityActive())
{
gradient +=
2*mu*(n & rheology.DEpsilonP().boundaryField()[patch.index()]);
}
// Thermal effects
if
(
patch.boundaryMesh().mesh().objectRegistry::
foundObject<thermalModel>("thermalProperties")
)
{
const thermalModel& thermo =
patch.boundaryMesh().mesh().objectRegistry::
lookupObject<thermalModel>("thermalProperties");
const fvPatchScalarField& threeKalpha =
patch.lookupPatchField<volScalarField, scalar>
(
"((threeK*rho)*alpha)"
);
if (!incremental)
{
const fvPatchScalarField& DT =
patch.lookupPatchField<volScalarField, scalar>("DT");
gradient += n*threeKalpha*DT;
}
else
{
const fvPatchScalarField& T =
patch.lookupPatchField<volScalarField, scalar>("T");
const scalarField T0 =
thermo.T0()().boundaryField()[patch.index()];
gradient += n*threeKalpha*(T - T0);
}
}
// Higher order non-linear terms
if
(
nonLinear == nonLinearGeometry::UPDATED_LAGRANGIAN
|| nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
)
{
// no extra relaxation
gradient -=
(n & (mu*(gradField & gradField.T())))
// + 0.5*n*lambda*(gradField && gradField);
+ 0.5*n*lambda*tr(gradField & gradField.T());
//- tensorial identity
//- tr(gradField & gradField.T())*I == (gradField && gradField)*I
if (fieldName == "DU" && nonLinear == "totalLagrangian")
{
// gradU is const in a time step
const fvPatchField<tensor>& gradU =
patch.lookupPatchField<volTensorField, tensor>("grad(U)");
if
(
incremental
&& nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN
)
{
// gradU is const in a time step
const fvPatchTensorField& gradU =
patch.lookupPatchField<volTensorField, tensor>("grad(U)");
gradient -=
(n &
(mu*( (gradField & gradU.T())
+ (gradU & gradField.T()) ))
)
+ 0.5*n*lambda*tr( (gradField & gradU.T())
+ (gradU & gradField.T()) );
}
}
gradient -=
(
n &
(
mu*
(
(gradField & gradU.T())
+ (gradU & gradField.T())
)
)
)
+ 0.5*n*lambda*tr
(
(gradField & gradU.T())
+ (gradU & gradField.T())
);
}
}
gradient /= (2.0*mu + lambda);
}
gradient /= (2.0*mu + lambda);
}
return tgradient;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

45
src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.H
View file

@ -35,6 +35,7 @@ SourceFiles
Author
Philip Cardiff UCD
Clean-up and re-factoring Hrvoje Jasak
\*---------------------------------------------------------------------------*/
@ -47,6 +48,7 @@ Author
#include "volFields.H"
#include "tmp.H"
#include "rheologyLaw.H"
#include "nonLinearGeometry.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -59,49 +61,38 @@ namespace Foam
class tractionBoundaryGradient
{
public:
//- Runtime type information
TypeName("tractionBoundaryGradient");
// Constructors
tractionBoundaryGradient()
{}
// Destructor
virtual ~tractionBoundaryGradient()
{}
// Member functions
// Static member functions
//- Return the boundary Cauchy traction corresponding to
// the given gradient
tmp<vectorField> traction
static tmp<vectorField> traction
(
const tensorField& gradField,
const word fieldName,
const word& workingFieldName, // Working variable
const word& integralFieldName, // Integrated displacement
const fvPatch& patch,
Switch orthotropic,
word nonLinear
) const;
const bool orthotropic,
const nonLinearGeometry::nonLinearType& nonLinear,
const bool incremental
);
// Operators
tmp<vectorField> operator()
//- Return surface-normal gradient given traction and pressure
static tmp<vectorField> snGrad
(
const vectorField& traction,
const scalarField& pressure,
const word fieldName,
const word& workingFieldName, // Working variable
const word& integralFieldName, // Integrated displacement
const fvPatch& patch,
Switch orthotropic,
word nonLinear
) const;
const bool orthotropic,
const nonLinearGeometry::nonLinearType& nonLinear,
const bool incremental
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

30
src/solidModels/contactModels/frictionContactModels/dirichletNeumannFriction/dirichletNeumannFriction.C
View file

@ -269,18 +269,24 @@ Foam::dirichletNeumannFriction::dirichletNeumannFriction
const fvPatch& slavePatch = mesh.boundary()[slavePatchIndex];
const fvPatchField<tensor>& gradField =
slavePatch.lookupPatchField<volTensorField, tensor>
("grad("+fieldName+")");
vectorField slaveShearTraction =
(I - sqr(slaveFaceNormals))
&
tractionBoundaryGradient().traction
(
gradField,
fieldName,
slavePatch,
orthotropic,
nonLinear
);
("grad(" + fieldName + ")");
bool incremental(fieldName == "DU");
vectorField slaveShearTraction
(
(I - sqr(slaveFaceNormals))
& tractionBoundaryGradient::traction
(
gradField,
fieldName,
"U",
slavePatch,
orthotropic,
nonLinearGeometry::nonLinearNames_[nonLinear],
incremental
)
);
// algorithm

336
src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C
View file

@ -117,57 +117,59 @@ dirichletNeumann::dirichletNeumann
contactFilePtr_(NULL)
{
// master proc open contact info file
if (Pstream::master())
if (Pstream::master())
{
word masterName = mesh_.boundary()[masterPatchID].name();
word slaveName = mesh_.boundary()[slavePatchID].name();
contactFilePtr_ =
new OFstream
(fileName("normalContact_"+masterName+"_"+slaveName+".txt"));
OFstream& contactFile = *contactFilePtr_;
int width = 20;
contactFile << "time";
contactFile.width(width);
contactFile << "iter";
contactFile.width(width);
contactFile << "contactFaces";
contactFile.width(width);
contactFile << "settledFaces";
contactFile.width(width);
contactFile << "minPen";
contactFile.width(width);
contactFile << "maxPress";
contactFile.width(width);
contactFile << "corrPoints" << endl;
}
word masterName = mesh_.boundary()[masterPatchID].name();
word slaveName = mesh_.boundary()[slavePatchID].name();
contactFilePtr_ =
new OFstream
(fileName("normalContact_"+masterName+"_"+slaveName+".txt"));
OFstream& contactFile = *contactFilePtr_;
int width = 20;
contactFile << "time";
contactFile.width(width);
contactFile << "iter";
contactFile.width(width);
contactFile << "contactFaces";
contactFile.width(width);
contactFile << "settledFaces";
contactFile.width(width);
contactFile << "minPen";
contactFile.width(width);
contactFile << "maxPress";
contactFile.width(width);
contactFile << "corrPoints" << endl;
}
// calc point points for missed vertices
if (correctMissedVertices_)
// calc point points for missed vertices
if (correctMissedVertices_)
{
calcSlavePointPoints(slavePatchID);
calcSlavePointPoints(slavePatchID);
}
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void dirichletNeumann::correct
(
const PrimitivePatch<face, List, pointField>& masterFaceZonePatch,
const PrimitivePatch<face, List, pointField>& slaveFaceZonePatch,
const intersection::algorithm alg,
const intersection::direction dir,
word fieldName,
const Switch orthotropic,
const word nonLinear,
vectorField& slaveFaceNormals,
GGIInterpolation< PrimitivePatch< face, List, pointField >,
PrimitivePatch< face, List, pointField >
>* ggiInterpolatorPtr
)
{
void dirichletNeumann::correct
(
const PrimitivePatch<face, List, pointField>& masterFaceZonePatch,
const PrimitivePatch<face, List, pointField>& slaveFaceZonePatch,
const intersection::algorithm alg,
const intersection::direction dir,
word fieldName,
const Switch orthotropic,
const word nonLinear,
vectorField& slaveFaceNormals,
GGIInterpolation
<
PrimitivePatch< face, List, pointField >,
PrimitivePatch< face, List, pointField >
>* ggiInterpolatorPtr
)
{
if (!settleContact_ || (iCorr_ < settleIterationNumber_))
{
{
//Info << "Correcting contact..." << flush;
const fvMesh& mesh = mesh_;
@ -186,11 +188,11 @@ dirichletNeumann::dirichletNeumann
>, PrimitivePatch<face, List, pointField>
> masterToSlavePatchToPatchInterpolator
(
masterFaceZonePatch, // from zone
slaveFaceZonePatch, // to zone
alg,
dir
);
masterFaceZonePatch, // from zone
slaveFaceZonePatch, // to zone
alg,
dir
);
//- calculate intersection distances
//- this is the slowest part of the contact correction
@ -205,43 +207,40 @@ dirichletNeumann::dirichletNeumann
(mesh.boundaryMesh()[slavePatchID()].nPoints(), 0.0);
label numCorrectedPoints = 0;
if (distanceMethod_ == "point")
{
{
// pointDistanceToIntersection() sometimes gives a seg fault when the
// momentum equation diverges
globalSlavePointPenetration
= masterToSlavePatchToPatchInterpolator.pointDistanceToIntersection();
globalSlavePointPenetration =
masterToSlavePatchToPatchInterpolator.pointDistanceToIntersection();
// get local point values from global values
forAll(slavePointPenetration, pointI)
{
// the local point values seem to be kept at the start
// of the global field
slavePointPenetration[pointI] =
globalSlavePointPenetration
[
pointI
];
{
// the local point values seem to be kept at the start
// of the global field
slavePointPenetration[pointI] =
globalSlavePointPenetration[pointI];
//- when the master surface surrounds the slave (like the pelvis and
// femur head) then
//- the slave penetration can sometimes calculate the distance through
// the femur head
//- to the pelvis which is wrong so I limit slavePenetration here
if (limitPenetration_)
{
if (slavePointPenetration[pointI] < penetrationLimit_)
{
slavePointPenetration[pointI] = 0.0;
}
}
}
//- when the master surface surrounds the slave (like the pelvis and
// femur head) then
//- the slave penetration can sometimes calculate the distance through
// the femur head
//- to the pelvis which is wrong so I limit slavePenetration here
if (limitPenetration_)
{
if (slavePointPenetration[pointI] < penetrationLimit_)
{
slavePointPenetration[pointI] = 0.0;
}
}
}
if (correctMissedVertices_)
{
scalarField& slavePointPenetration_ = slavePointPenetration;
{
scalarField& slavePointPenetration_ = slavePointPenetration;
# include "iterativePenaltyCorrectMissedVertices.H"
}
}
minSlavePointPenetration = gMin(globalSlavePointPenetration);
@ -485,7 +484,7 @@ dirichletNeumann::dirichletNeumann
// remove tengential component
slaveDisp_ = slaveFaceNormals*(slaveFaceNormals & slaveDisp_);
//--------Try to remove any oscillations----------//
//--------Try to remove any oscillations----------//
if (oscillationCorr_)
{
if (smoothingSteps_ < 1)
@ -502,15 +501,16 @@ dirichletNeumann::dirichletNeumann
(mesh.boundaryMesh()[slavePatchIndex].nPoints(), vector::zero);
for (int i=0; i<smoothingSteps_; i++)
{
slaveDispPoints =
localSlaveInterpolator.faceToPointInterpolate<vector>(slaveDisp_);
slaveDisp_ =
localSlaveInterpolator.pointToFaceInterpolate<vector>(slaveDispPoints);
{
slaveDispPoints =
localSlaveInterpolator.faceToPointInterpolate<vector>(slaveDisp_);
// make sure no tangential component
slaveDisp_ = slaveFaceNormals*(slaveFaceNormals & slaveDisp_);
}
slaveDisp_ = localSlaveInterpolator.pointToFaceInterpolate<vector>
(slaveDispPoints);
// make sure no tangential component
slaveDisp_ = slaveFaceNormals*(slaveFaceNormals & slaveDisp_);
}
}
@ -530,66 +530,80 @@ dirichletNeumann::dirichletNeumann
// calculate current slave traction
{
const fvPatch& slavePatch = mesh.boundary()[slavePatchIndex];
const fvPatchField<tensor>& gradField =
slavePatch.lookupPatchField<volTensorField, tensor>
("grad("+fieldName+")");
slavePressure_ = tractionBoundaryGradient().traction
const fvPatch& slavePatch = mesh.boundary()[slavePatchIndex];
const fvPatchField<tensor>& gradField =
slavePatch.lookupPatchField<volTensorField, tensor>
(
"grad(" + fieldName + ")"
);
bool incremental(fieldName == "DU");
slavePressure_ = tractionBoundaryGradient::traction
(
gradField,
fieldName,
slavePatch,
orthotropic,
nonLinear
);
gradField,
fieldName,
"U",
slavePatch,
orthotropic,
nonLinearGeometry::nonLinearNames_[nonLinear],
incremental
);
// set traction to zero on faces not in contact
// and set tensile stresses to zero
//const scalar maxMagSlavePressure = gMax(mag(slavePressure_));
scalar maxMagSlavePressure = 0.0;
if (slavePressure_.size() > 0)
{
maxMagSlavePressure = max(mag(slavePressure_));
}
reduce(maxMagSlavePressure, maxOp<scalar>());
forAll(touchFraction_, facei)
// set traction to zero on faces not in contact
// and set tensile stresses to zero
//const scalar maxMagSlavePressure = gMax(mag(slavePressure_));
scalar maxMagSlavePressure = 0.0;
if (slavePressure_.size() > 0)
{
if (touchFraction_[facei] < SMALL
||
( (slaveFaceNormals[facei] & slavePressure_[facei])
> 1e-3*maxMagSlavePressure)
)
{
slavePressure_[facei] = vector::zero;
slaveValueFrac_[facei] = symmTensor::zero;
slaveDisp_[facei] =
slaveFaceNormals[facei]
*(slaveFaceNormals[facei]&oldSlaveDisp[facei]);
}
maxMagSlavePressure = max(mag(slavePressure_));
}
reduce(maxMagSlavePressure, maxOp<scalar>());
// relax traction
slavePressure_ =
relaxFactor_*slavePressure_ + (1-relaxFactor_)*oldSlavePressure_;
// remove any shears
slavePressure_ = slaveFaceNormals*(slaveFaceNormals & slavePressure_);
// limit pressure to help convergence
if (limitPressure_)
forAll(touchFraction_, facei)
{
forAll(slavePressure_, facei)
{
if ( (slaveFaceNormals[facei]&slavePressure_[facei]) < -pressureLimit_)
if
(
touchFraction_[facei] < SMALL
|| (
(slaveFaceNormals[facei] & slavePressure_[facei])
> 1e-3*maxMagSlavePressure
)
)
{
slavePressure_[facei] = -pressureLimit_*slaveFaceNormals[facei];
slavePressure_[facei] = vector::zero;
slaveValueFrac_[facei] = symmTensor::zero;
slaveDisp_[facei] = slaveFaceNormals[facei]
*(slaveFaceNormals[facei]&oldSlaveDisp[facei]);
}
}
// relax traction
slavePressure_ =
relaxFactor_*slavePressure_ + (1-relaxFactor_)*oldSlavePressure_;
// remove any shears
slavePressure_ = slaveFaceNormals*(slaveFaceNormals & slavePressure_);
// limit pressure to help convergence
if (limitPressure_)
{
forAll(slavePressure_, facei)
{
if
(
(slaveFaceNormals[facei] & slavePressure_[facei])
< -pressureLimit_
)
{
slavePressure_[facei] =
-pressureLimit_*slaveFaceNormals[facei];
}
}
}
// update old slave traction
oldSlavePressure_ = slavePressure_;
// update old slave traction
oldSlavePressure_ = slavePressure_;
}
// in parallel, the log is poluted with warnings that
@ -604,36 +618,36 @@ dirichletNeumann::dirichletNeumann
}
reduce(maxMagMasterTraction, maxOp<scalar>());
// under-relax value fraction
slaveValueFrac_ =
relaxFactor_*slaveValueFrac_ + (1-relaxFactor_)*oldSlaveValueFrac_;
oldSlaveValueFrac_ = slaveValueFrac_;
// under-relax value fraction
slaveValueFrac_ =
relaxFactor_*slaveValueFrac_ + (1-relaxFactor_)*oldSlaveValueFrac_;
oldSlaveValueFrac_ = slaveValueFrac_;
//--------MASTER PROCS WRITES CONTACT INFO FILE----------//
if (Pstream::master() && (contactIterNum_ % infoFreq_ == 0))
{
OFstream& contactFile = *contactFilePtr_;
int width = 20;
contactFile << mesh.time().value();
contactFile.width(width);
contactFile << contactIterNum_;
contactFile.width(width);
contactFile << numSlaveContactFaces;
contactFile.width(width);
contactFile << numSlaveSettledFaces;
contactFile.width(width);
contactFile << minSlavePointPenetration;
contactFile.width(width);
contactFile << maxMagMasterTraction;
contactFile.width(width);
contactFile << numCorrectedPoints;
if (aitkenRelaxation_)
{
contactFile.width(width);
contactFile << aitkenTheta_;
}
contactFile << endl;
}
//--------MASTER PROCS WRITES CONTACT INFO FILE----------//
if (Pstream::master() && (contactIterNum_ % infoFreq_ == 0))
{
OFstream& contactFile = *contactFilePtr_;
int width = 20;
contactFile << mesh.time().value();
contactFile.width(width);
contactFile << contactIterNum_;
contactFile.width(width);
contactFile << numSlaveContactFaces;
contactFile.width(width);
contactFile << numSlaveSettledFaces;
contactFile.width(width);
contactFile << minSlavePointPenetration;
contactFile.width(width);
contactFile << maxMagMasterTraction;
contactFile.width(width);
contactFile << numCorrectedPoints;
if (aitkenRelaxation_)
{
contactFile.width(width);
contactFile << aitkenTheta_;
}
contactFile << endl;
}
//Info << "\tdone" << endl;
}

27
src/solidModels/fvPatchFields/fixedDisplacementOrSolidTraction/fixedDisplacementOrSolidTractionFvPatchVectorField.C
View file

@ -200,21 +200,22 @@ void fixedDisplacementOrSolidTractionFvPatchVectorField::updateCoeffs()
if ( mag(timeSeries_(this->db().time().timeOutputValue())) < SMALL)
{
// traction boundary
// traction boundary
// set valueFraction to zero
this->valueFraction() = symmTensor::zero;
// set valueFraction to zero
this->valueFraction() = symmTensor::zero;
// set gradient to enfore specified traction
refGrad() = tractionBoundaryGradient()
(
traction_,
pressure_,
word(fieldName_),
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
)();
// set gradient to enfore specified traction
refGrad() = tractionBoundaryGradient::snGrad
(
traction_,
pressure_,
fieldName_,
"U",
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
);
}
else
{

53
src/solidModels/fvPatchFields/fixedDisplacementZeroShear/fixedDisplacementZeroShearFvPatchVectorField.C
View file

@ -130,31 +130,32 @@ fixedDisplacementZeroShearFvPatchVectorField
<< "non-orthogonal correction to be right" << endl;
}
this->refGrad() = vector::zero;
this->refGrad() = vector::zero;
vectorField n = patch().nf();
this->valueFraction() = sqr(n);
vectorField n = patch().nf();
this->valueFraction() = sqr(n);
if (dict.found("value"))
if (dict.found("value"))
{
Field<vector>::operator=(vectorField("value", dict, p.size()));
Field<vector>::operator=(vectorField("value", dict, p.size()));
}
else
else
{
FatalError << "value entry not found for patch " << patch().name()
<< exit(FatalError);
FatalError << "value entry not found for patch " << patch().name()
<< exit(FatalError);
}
this->refValue() = *this;
Field<vector> normalValue = transform(valueFraction(), refValue());
this->refValue() = *this;
Field<vector> gradValue =
this->patchInternalField() + refGrad()/this->patch().deltaCoeffs();
Field<vector> normalValue = transform(valueFraction(), refValue());
Field<vector> transformGradValue =
transform(I - valueFraction(), gradValue);
Field<vector> gradValue =
this->patchInternalField() + refGrad()/this->patch().deltaCoeffs();
Field<vector>::operator=(normalValue + transformGradValue);
Field<vector> transformGradValue =
transform(I - valueFraction(), gradValue);
Field<vector>::operator=(normalValue + transformGradValue);
}
@ -223,15 +224,19 @@ void fixedDisplacementZeroShearFvPatchVectorField::updateCoeffs()
// this->valueFraction() = sqr(nCurrent);
// }
refGrad() = tractionBoundaryGradient()
(
vectorField(patch().size(), vector::zero),
scalarField(patch().size(), 0.0),
word(fieldName_),
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
)();
bool incremental(fieldName_ == "DU");
refGrad() = tractionBoundaryGradient::snGrad
(
vectorField(patch().size(), vector::zero),
scalarField(patch().size(), 0.0),
fieldName_,
"U",
patch(),
orthotropic_,
nonLinear_,
incremental
);
directionMixedFvPatchVectorField::updateCoeffs();
}

63
src/solidModels/fvPatchFields/fixedDisplacementZeroShearOrSolidTraction/fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField.C
View file

@ -201,48 +201,51 @@ void fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField::rmap
}
void fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField::updateCoeffs()
void
fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField::updateCoeffs()
{
if (this->updated())
{
return;
}
if ( mag(timeSeries_(this->db().time().timeOutputValue())) < SMALL)
if (mag(timeSeries_(this->db().time().timeOutputValue())) < SMALL)
{
// traction boundary
// traction boundary
// set valueFraction to zero
this->valueFraction() = symmTensor::zero;
// set valueFraction to zero
this->valueFraction() = symmTensor::zero;
// set gradient to enfore specified traction
refGrad() = tractionBoundaryGradient()
(
traction_,
pressure_,
word(fieldName_),
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
)();
// set gradient to enfore specified traction
refGrad() = tractionBoundaryGradient::snGrad
(
traction_,
pressure_,
fieldName_,
"U",
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
);
}
else
{
// fixed displacement zero shear
{
// fixed displacement zero shear
// set valueFraction to fix normal
this->valueFraction() = sqr(fixedNormal_);
// set valueFraction to fix normal
this->valueFraction() = sqr(fixedNormal_);
// force zero shear stresses
refGrad() = tractionBoundaryGradient()
(
vectorField(traction_.size(), vector::zero),
scalarField(traction_.size(), 0.0),
word(fieldName_),
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
)();
// force zero shear stresses
refGrad() = tractionBoundaryGradient::snGrad
(
vectorField(traction_.size(), vector::zero),
scalarField(traction_.size(), 0.0),
fieldName_,
"U",
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
);
// set displacement
refValue() = displacement_;
@ -256,7 +259,7 @@ void fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField::updateCoeffs()
void fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField::write
(
Ostream& os
) const
) const
{
directionMixedFvPatchVectorField::write(os);
os.writeKeyword("nonLinear")

205
src/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C
View file

@ -767,6 +767,8 @@ void solidContactFvPatchVectorField::updateCoeffs()
}
// set boundary conditions
bool incremental(fieldName_ == "DU");
if (!master_)
{
// set refValue, refGrad and valueFraction
@ -778,16 +780,18 @@ void solidContactFvPatchVectorField::updateCoeffs()
//refGrad - set traction
refGrad() =
tractionBoundaryGradient()
tractionBoundaryGradient::snGrad
(
frictionContactModelPtr_->slaveTraction()
+ normalContactModelPtr_->slavePressure(),
scalarField(patch().size(),0.0),
word(fieldName_),
fieldName_,
"U",
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
)();
nonLinear_,
incremental
);
//valueFraction
valueFraction() =
@ -805,32 +809,35 @@ void solidContactFvPatchVectorField::updateCoeffs()
{
// set to master to traction free if it is rigid
refGrad() =
tractionBoundaryGradient()
tractionBoundaryGradient::snGrad
(
vectorField(patch().size(),vector::zero),
scalarField(patch().size(),0.0),
word(fieldName_),
fieldName_,
"U",
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
)();
nonLinear_,
incremental
);
}
else
{
refGrad() =
tractionBoundaryGradient()
refGrad() = tractionBoundaryGradient::snGrad
(
interpolateSlaveToMaster
(
interpolateSlaveToMaster
(
-frictionContactModelPtr_->slaveTraction()
-normalContactModelPtr_->slavePressure()
),
scalarField(patch().size(),0.0),
word(fieldName_),
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
)();
-frictionContactModelPtr_->slaveTraction()
-normalContactModelPtr_->slavePressure()
),
scalarField(patch().size(),0.0),
fieldName_,
"U",
patch(),
orthotropic_,
nonLinear_,
incremental
);
}
}
} // if correction freqeuncy
@ -840,28 +847,33 @@ void solidContactFvPatchVectorField::updateCoeffs()
else
{
// set refGrad to traction free for master and slave
bool incremental(fieldName_ == "DU");
refGrad() =
tractionBoundaryGradient()
tractionBoundaryGradient::snGrad
(
vectorField(patch().size(),vector::zero),
scalarField(patch().size(),0.0),
word(fieldName_),
fieldName_,
"U",
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
)();
nonLinear_,
incremental
);
}
directionMixedFvPatchVectorField::updateCoeffs();
}
// Interpolate traction from slave to master
tmp<vectorField> solidContactFvPatchVectorField::interpolateSlaveToMaster
(
const vectorField slaveField
)
{
if (!master_)
if (!master_)
{
FatalError
<< "only the master may call the function "
@ -1315,85 +1327,96 @@ snGrad() const
//- Increment of dissipated energy due to friction
tmp<scalarField> solidContactFvPatchVectorField::Qc() const
{
tmp<scalarField> tQc(new scalarField(patch().size(), 0.0));
scalarField& Qc = tQc();
tmp<scalarField> tQc(new scalarField(patch().size(), 0.0));
scalarField& Qc = tQc();
// Integrate energy using trapezoidal rule
// 0.5*averageForce*incrementOfDisplacement
// Integrate energy using trapezoidal rule
// 0.5*averageForce*incrementOfDisplacement
// displacement increment
const vectorField& curDU = *this;
// displacement increment
const vectorField& curDU = *this;
// average force
const fvPatchField<tensor>& gradField =
patch().lookupPatchField<volTensorField, tensor>("grad("+fieldName_+")");
// current Cauchy traction for nonlinear
const vectorField curTrac =
tractionBoundaryGradient().traction
(
gradField,
fieldName_,
patch(),
orthotropic_,
word(nonLinearGeometry::nonLinearNames_[nonLinear_])
);
// average force
const fvPatchField<tensor>& gradField =
patch().lookupPatchField<volTensorField, tensor>
(
"grad(" + fieldName_ + ")"
);
vectorField Sf = patch().Sf();
if (nonLinear_ != nonLinearGeometry::OFF
&& nonLinear_ == nonLinearGeometry::TOTAL_LAGRANGIAN)
// current Cauchy traction for nonlinear
bool incremental(fieldName_ == "DU");
const vectorField curTrac =
tractionBoundaryGradient::traction
(
gradField,
fieldName_,
"U",
patch(),
orthotropic_,
nonLinear_,
incremental
);
vectorField Sf = patch().Sf();
if
(
nonLinear_ != nonLinearGeometry::OFF
&& nonLinear_ == nonLinearGeometry::TOTAL_LAGRANGIAN
)
{
// current areas
const fvPatchField<tensor>& gradU =
patch().lookupPatchField<volTensorField, tensor>("grad(U)");
tensorField F = I + gradField + gradU;
tensorField Finv = inv(F);
scalarField J = det(F);
Sf = J*(Finv & Sf);
}
const scalarField magSf = mag(Sf);
const vectorField curForce = magSf * curTrac;
const fvPatchField<symmTensor>& oldSigma =
patch().lookupPatchField<volSymmTensorField, symmTensor>("sigma_0");
vectorField oldSf = patch().Sf();
if (nonLinear_ != nonLinearGeometry::OFF)
{
// current areas
//tensorField F = I + gradField;
if (nonLinear_ == nonLinearGeometry::TOTAL_LAGRANGIAN)
{
const fvPatchField<tensor>& gradU =
// current areas
const fvPatchField<tensor>& gradU =
patch().lookupPatchField<volTensorField, tensor>("grad(U)");
tensorField F = I + gradU;
tensorField Finv = inv(F);
scalarField J = det(F);
// rotate initial reference area to area of previous time-step
oldSf = J*(Finv & Sf);
}
else if (nonLinear_ == nonLinearGeometry::UPDATED_LAGRANGIAN)
tensorField F = I + gradField + gradU;
tensorField Finv = inv(F);
scalarField J = det(F);
Sf = J*(Finv & Sf);
}
const scalarField magSf = mag(Sf);
const vectorField curForce = magSf*curTrac;
const fvPatchField<symmTensor>& oldSigma =
patch().lookupPatchField<volSymmTensorField, symmTensor>("sigma_0");
vectorField oldSf = patch().Sf();
if (nonLinear_ != nonLinearGeometry::OFF)
{
// current areas
//tensorField F = I + gradField;
if (nonLinear_ == nonLinearGeometry::TOTAL_LAGRANGIAN)
{
tensorField F = I + gradField;
tensorField Finv = inv(F);
scalarField J = det(F);
// rotate current area to area of previous time-step
oldSf = (F & Sf)/J;
const fvPatchField<tensor>& gradU =
patch().lookupPatchField<volTensorField, tensor>("grad(U)");
tensorField F = I + gradU;
tensorField Finv = inv(F);
scalarField J = det(F);
// rotate initial reference area to area of previous time-step
oldSf = J*(Finv & Sf);
}
else
else if (nonLinear_ == nonLinearGeometry::UPDATED_LAGRANGIAN)
{
FatalError << "solidContact::Qc() nonLinear type not known!"
<< exit(FatalError);
tensorField F = I + gradField;
tensorField Finv = inv(F);
scalarField J = det(F);
// rotate current area to area of previous time-step
oldSf = (F & Sf)/J;
}
else
{
FatalError << "solidContact::Qc() nonLinear type not known!"
<< exit(FatalError);
}
}
const vectorField oldForce = oldSf & oldSigma;
const vectorField avForce = 0.5*(oldForce + curForce);
const vectorField oldForce = oldSf & oldSigma;
const vectorField avForce = 0.5*(oldForce + curForce);
// Increment of dissiptaed frictional energy for this timestep
Qc = mag(0.5*(avForce & curDU));
// Increment of dissiptaed frictional energy for this timestep
Qc = mag(0.5*(avForce & curDU));
return tQc;
return tQc;
}

26
src/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.C
View file

@ -210,15 +210,19 @@ void solidTractionFvPatchVectorField::updateCoeffs()
return;
}
gradient() = tractionBoundaryGradient()
(
traction_,
pressure_,
word(fieldName_),
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
)();
bool incremental(fieldName_ == "DU");
gradient() = tractionBoundaryGradient::snGrad
(
traction_,
pressure_,
fieldName_,
"U",
patch(),
orthotropic_,
nonLinear_,
incremental
);
fixedGradientFvPatchVectorField::updateCoeffs();
}
@ -243,14 +247,14 @@ void solidTractionFvPatchVectorField::evaluate(const Pstream::commsTypes)
//- non-orthogonal correction vectors
vectorField k = delta - n*(n&delta);
Field<vector>::operator=
vectorField::operator=
(
this->patchInternalField()
+ (k&gradField.patchInternalField())
+ gradient()/this->patch().deltaCoeffs()
);
fvPatchField<vector>::evaluate();
fvPatchVectorField::evaluate();
}
// Write

2
src/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.H
View file

@ -53,7 +53,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class solidTractionFvPatchVectorField Declaration
Class solidTractionFvPatchVectorField Declaration
\*---------------------------------------------------------------------------*/
class solidTractionFvPatchVectorField

20
src/solidModels/fvPatchFields/solidTractionFree/solidTractionFreeFvPatchVectorField.C
View file

@ -75,9 +75,9 @@ solidTractionFreeFvPatchVectorField
if (dict.found("nonLinear"))
{
nonLinear_ = nonLinearGeometry::nonLinearNames_.read
(
dict.lookup("nonLinear")
);
(
dict.lookup("nonLinear")
);
if (nonLinear_ == nonLinearGeometry::UPDATED_LAGRANGIAN)
{
@ -191,15 +191,19 @@ void solidTractionFreeFvPatchVectorField::updateCoeffs()
return;
}
gradient() = tractionBoundaryGradient()
bool incremental(fieldName_ == "DU");
gradient() = tractionBoundaryGradient::snGrad
(
vectorField(patch().size(), vector::zero),
scalarField(patch().size(), 0.0),
word(fieldName_),
fieldName_,
"U",
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
)();
nonLinear_,
incremental
);
fixedGradientFvPatchVectorField::updateCoeffs();
}
@ -226,7 +230,7 @@ void solidTractionFreeFvPatchVectorField::evaluate(const Pstream::commsTypes)
Field<vector>::operator=
(
this->patchInternalField()
+ (k&gradField.patchInternalField())
+ (k & gradField.patchInternalField())
+ gradient()/this->patch().deltaCoeffs()
);

2
src/solidModels/fvPatchFields/solidTractionFree/solidTractionFreeFvPatchVectorField.H
View file

@ -66,7 +66,7 @@ class solidTractionFreeFvPatchVectorField
//- Name of the displacement field
const word fieldName_;
//- if it is a non linear solver
//- Is it a non linear solver
nonLinearGeometry::nonLinearType nonLinear_;
//- Is it an orthropic solver

10
src/solidModels/fvPatchFields/timeVaryingFixedDisplacementZeroShear/timeVaryingFixedDisplacementZeroShearFvPatchVectorField.C
View file

@ -233,14 +233,18 @@ void timeVaryingFixedDisplacementZeroShearFvPatchVectorField::updateCoeffs()
//vectorField n = patch().nf();
//this->valueFraction() = sqr(n);
refGrad() = tractionBoundaryGradient()
bool incremental(fieldName_ == "DU");
refGrad() = tractionBoundaryGradient::snGrad
(
vectorField(patch().size(), vector::zero),
scalarField(patch().size(), 0),
word(fieldName_),
fieldName_,
"U",
patch(),
orthotropic_,
nonLinearGeometry::nonLinearNames_[nonLinear_]
nonLinear_,
incremental
);
directionMixedFvPatchVectorField::updateCoeffs();

2
src/solidModels/nonLinearGeometry/nonLinearGeometry.H
View file

@ -53,7 +53,7 @@ class nonLinearGeometry
{
public:
//- non-linear solver options
//- Non-linear solver options
enum nonLinearType
{
OFF,

8
tutorials/solidMechanics/elasticIncrSolidFoam/incrPlateHole/0/DU
View file

@ -28,10 +28,10 @@ boundaryField
right
{
type timeVaryingSolidTraction;
nonLinear off;
outOfBounds clamp;
fileName "$FOAM_CASE/constant/timeVsRightTraction";
type timeVaryingSolidTraction;
nonLinear off;
outOfBounds clamp;
fileName "$FOAM_CASE/constant/timeVsRightTraction";
}
down

2
tutorials/solidMechanics/elasticIncrSolidFoam/incrPlateHole/constant/polyMesh/boundary
View file

@ -1,7 +1,7 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.0 |
| \\ / O peration | Version: 3.0 |
| \\ / A nd | Web: http://www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/

11
tutorials/solidMechanics/elasticSolidFoam/plateHole/0/U
View file

@ -27,19 +27,20 @@ boundaryField
}
right
{
type analyticalPlateHoleTraction;
// type analyticalPlateHoleTraction;
type solidTraction;
traction uniform (0 0 0);
pressure uniform 1e7;
}
down
{
type symmetryPlane;
}
up
{
type analyticalPlateHoleTraction;
// type analyticalPlateHoleTraction;
type solidTractionFree;
}
hole
{
type solidTractionFree;

2
tutorials/solidMechanics/elasticSolidFoam/plateHole/constant/polyMesh/boundary
View file

@ -1,7 +1,7 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.0 |
| \\ / O peration | Version: 3.0 |
| \\ / A nd | Web: http://www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/