diff --git a/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C b/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C index 7a0accb32..a7a78cb44 100644 --- a/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C +++ b/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C @@ -116,10 +116,10 @@ solidCohesiveFvPatchVectorField::solidCohesiveFvPatchVectorField nonLinear_(nonLinearGeometry::OFF), orthotropic_(false) { - Info << "Creating solidCohesive patch" << nl - << "\tOnly Dugdale law currently available!" << endl; + Info<< "Creating solidCohesive patch" << nl + << "\tOnly Dugdale law currently available!" << endl; - //- check if traction boundary is for non linear solver + // Check if traction boundary is for non linear solver if (dict.found("nonLinear")) { nonLinear_ = nonLinearGeometry::nonLinearNames_.read @@ -127,23 +127,23 @@ solidCohesiveFvPatchVectorField::solidCohesiveFvPatchVectorField dict.lookup("nonLinear") ); - if (nonLinear_ == nonLinearGeometry::UPDATED_LAGRANGIAN) - { - Info << "\tnonLinear set to updated Lagrangian" - << endl; - } - else if (nonLinear_ == nonLinearGeometry::TOTAL_LAGRANGIAN) - { - Info << "\tnonLinear set to total Lagrangian" - << endl; - } - } + if (nonLinear_ == nonLinearGeometry::UPDATED_LAGRANGIAN) + { + Info << "\tnonLinear set to updated Lagrangian" + << endl; + } + else if (nonLinear_ == nonLinearGeometry::TOTAL_LAGRANGIAN) + { + Info << "\tnonLinear set to total Lagrangian" + << endl; + } + } if (dict.found("orthotropic")) - { - orthotropic_ = Switch(dict.lookup("orthotropic")); - Info << "\t\torthotropic set to " << orthotropic_ << endl; - } + { + orthotropic_ = Switch(dict.lookup("orthotropic")); + Info << "\t\torthotropic set to " << orthotropic_ << endl; + } // if (minUnloadingSeparationDistance_ < 0.001) // { @@ -216,73 +216,73 @@ solidCohesiveFvPatchVectorField::solidCohesiveFvPatchVectorField // } if (dict.found("cracked")) - { + { // cracked_ = Field("cracked", dict, p.size()); - cracked_ = Field("cracked", dict, p.size()); - } + cracked_ = Field("cracked", dict, p.size()); + } if (dict.found("curTractionN")) - { - curTractionN_ = scalarField("curTractionN", dict, p.size()); - } + { + curTractionN_ = scalarField("curTractionN", dict, p.size()); + } if (dict.found("curTractionS")) - { - curTractionS_ = scalarField("curTractionS", dict, p.size()); - } + { + curTractionS_ = scalarField("curTractionS", dict, p.size()); + } if (dict.found("oldTractionN")) - { - oldTractionN_ = scalarField("oldTractionN", dict, p.size()); - } + { + oldTractionN_ = scalarField("oldTractionN", dict, p.size()); + } if (dict.found("oldTractionS")) - { - oldTractionS_ = scalarField("oldTractionS", dict, p.size()); - } + { + oldTractionS_ = scalarField("oldTractionS", dict, p.size()); + } if (dict.found("deltaN")) - { - deltaN_ = scalarField("deltaN", dict, p.size()); - } + { + deltaN_ = scalarField("deltaN", dict, p.size()); + } if (dict.found("deltaS")) - { - deltaS_ = scalarField("deltaS", dict, p.size()); - } + { + deltaS_ = scalarField("deltaS", dict, p.size()); + } if (dict.found("oldDeltaN")) - { - oldDeltaN_ = scalarField("oldDeltaN", dict, p.size()); - } + { + oldDeltaN_ = scalarField("oldDeltaN", dict, p.size()); + } if (dict.found("oldDeltaS")) - { - oldDeltaS_ = scalarField("oldDeltaS", dict, p.size()); - } + { + oldDeltaS_ = scalarField("oldDeltaS", dict, p.size()); + } if (dict.found("curDeltaN")) - { - curDeltaN_ = scalarField("curDeltaN", dict, p.size()); - } + { + curDeltaN_ = scalarField("curDeltaN", dict, p.size()); + } if (dict.found("curDeltaS")) - { - curDeltaS_ = scalarField("curDeltaS", dict, p.size()); - } + { + curDeltaS_ = scalarField("curDeltaS", dict, p.size()); + } if (dict.found("unloadingDeltaEff")) - { - unloadingDeltaEff_ = scalarField("unloadingDeltaEff", dict, p.size()); - } + { + unloadingDeltaEff_ = scalarField("unloadingDeltaEff", dict, p.size()); + } if (dict.found("currentGI")) - { - currentGI_ = scalarField("currentGI", dict, p.size()); - } + { + currentGI_ = scalarField("currentGI", dict, p.size()); + } if (dict.found("currentGII")) - { - currentGII_ = scalarField("currentGII", dict, p.size()); - } + { + currentGII_ = scalarField("currentGII", dict, p.size()); + } if (dict.found("oldGI")) - { - oldGI_ = scalarField("oldGI", dict, p.size()); - } + { + oldGI_ = scalarField("oldGI", dict, p.size()); + } if (dict.found("oldGII")) - { - oldGII_ = scalarField("oldGII", dict, p.size()); - } + { + oldGII_ = scalarField("oldGII", dict, p.size()); + } } @@ -415,14 +415,15 @@ tmp solidCohesiveFvPatchVectorField::crackingAndDamage() const ( new scalarField(size(), 1.0) ); + scalarField& cad = tcrackingAndDamage(); - forAll(tcrackingAndDamage(), facei) - { - if (cracked_[facei]) - { - tcrackingAndDamage()[facei] = 2.0; - } - } + forAll(cad, facei) + { + if (cracked_[facei]) + { + cad[facei] = 2.0; + } + } return tcrackingAndDamage; } @@ -434,11 +435,12 @@ tmp solidCohesiveFvPatchVectorField::GI() const ( new scalarField(size(), 0.0) ); + scalarField& GI = tGI(); - forAll(tGI(), facei) - { - tGI()[facei] = currentGI_[facei]; - } + forAll(GI, facei) + { + GI[facei] = currentGI_[facei]; + } return tGI; } @@ -450,41 +452,44 @@ tmp solidCohesiveFvPatchVectorField::GII() const ( new scalarField(size(), 0.0) ); + scalarField& GII = tGII(); - forAll(tGII(), facei) - { - tGII()[facei] = currentGII_[facei]; - } + forAll(GII, facei) + { + GII[facei] = currentGII_[facei]; + } return tGII; } + bool solidCohesiveFvPatchVectorField::cracking() { - const fvMesh& mesh = patch().boundaryMesh().mesh(); - const crackerFvMesh& crackerMesh = refCast(mesh); + const fvMesh& mesh = patch().boundaryMesh().mesh(); + const crackerFvMesh& crackerMesh = refCast(mesh); - // global fields - Field globalCracked = - crackerMesh.globalCrackField(cracked_); + // global fields + Field globalCracked = + crackerMesh.globalCrackField(cracked_); - label sumDamaged = 0; - label sumCracked = 0; - forAll(globalCracked, facei) - { - if (globalCracked[facei] > 0.0) - { - sumCracked++; - } - else - { - sumDamaged++; - } - } - Info << "\t\tThere are " << sumDamaged << " damaged faces" << endl; - Info << "\t\tThere are " << sumCracked << " cracked faces" << endl; + label sumDamaged = 0; + label sumCracked = 0; - return false; + forAll(globalCracked, facei) + { + if (globalCracked[facei] > 0.0) + { + sumCracked++; + } + else + { + sumDamaged++; + } + } + Info<< "\t\tThere are " << sumDamaged << " damaged faces" << nl + << "\t\tThere are " << sumCracked << " cracked faces" << endl; + + return false; } @@ -534,84 +539,13 @@ void solidCohesiveFvPatchVectorField::autoMap // the traction_ field after mapping and call updateCoeffs. // For here, we will just set traction to zero. - if ( (patch().size()==1) && (nNewFaces == 1) ) + if ( (patch().size() == 1) && (nNewFaces == 1) ) { label i=0; this->valueFraction()[i] = symmTensor::zero; - traction_[i] = vector::zero; - cracked_[i] = false; - curTractionN_[i] = 0.0; - oldTractionN_[i] = 0.0; - curTractionS_[i] = 0.0; - oldTractionS_[i] = 0.0; - deltaN_[i] = 0.0; - oldDeltaN_[i] = 0.0; - deltaS_[i] = 0.0; - oldDeltaS_[i] = 0.0; - curDeltaN_[i] = 0.0; - curDeltaS_[i] = 0.0; - unloadingDeltaEff_[i] = 0.0; - currentGI_[i] = 0.0; - oldGI_[i] = 0.0; - currentGII_[i] = 0.0; - oldGII_[i] = 0.0; - } - else if ( (patch().size()==2) && (nNewFaces == 1) ) - { - label i=1; - this->valueFraction()[i] = symmTensor::zero; + traction_[i] = vector::zero; cracked_[i] = false; - curTractionN_[i] = 0.0; - oldTractionN_[i] = 0.0; - curTractionS_[i] = 0.0; - oldTractionS_[i] = 0.0; - deltaN_[i] = 0.0; - oldDeltaN_[i] = 0.0; - deltaS_[i] = 0.0; - oldDeltaS_[i] = 0.0; - curDeltaN_[i] = 0.0; - curDeltaS_[i] = 0.0; - unloadingDeltaEff_[i] = 0.0; - currentGI_[i] = 0.0; - oldGI_[i] = 0.0; - currentGII_[i] = 0.0; - oldGII_[i] = 0.0; - } - else if ( (patch().size()==2) && (nNewFaces == 2) ) - { - label i=0; - this->valueFraction()[i] = symmTensor::zero; - traction_[i] = vector::zero; - i=1; - this->valueFraction()[i] = symmTensor::zero; - traction_[i] = vector::zero; - cracked_[i] = false; - curTractionN_[i] = 0.0; - oldTractionN_[i] = 0.0; - curTractionS_[i] = 0.0; - oldTractionS_[i] = 0.0; - deltaN_[i] = 0.0; - oldDeltaN_[i] = 0.0; - deltaS_[i] = 0.0; - oldDeltaS_[i] = 0.0; - curDeltaN_[i] = 0.0; - curDeltaS_[i] = 0.0; - unloadingDeltaEff_[i] = 0.0; - currentGI_[i] = 0.0; - oldGI_[i] = 0.0; - currentGII_[i] = 0.0; - oldGII_[i] = 0.0; - } - else - { - for (label i = 1; i < patch().size(); i++) - { - if (addressing[i] == 0) - { - this->valueFraction()[i] = symmTensor::zero; - traction_[i] = vector::zero; - cracked_[i] = false; curTractionN_[i] = 0.0; oldTractionN_[i] = 0.0; curTractionS_[i] = 0.0; @@ -627,6 +561,81 @@ void solidCohesiveFvPatchVectorField::autoMap oldGI_[i] = 0.0; currentGII_[i] = 0.0; oldGII_[i] = 0.0; + } + else if ( (patch().size() == 2) && (nNewFaces == 1) ) + { + label i=1; + this->valueFraction()[i] = symmTensor::zero; + traction_[i] = vector::zero; + cracked_[i] = false; + + curTractionN_[i] = 0.0; + oldTractionN_[i] = 0.0; + curTractionS_[i] = 0.0; + oldTractionS_[i] = 0.0; + deltaN_[i] = 0.0; + oldDeltaN_[i] = 0.0; + deltaS_[i] = 0.0; + oldDeltaS_[i] = 0.0; + curDeltaN_[i] = 0.0; + curDeltaS_[i] = 0.0; + unloadingDeltaEff_[i] = 0.0; + currentGI_[i] = 0.0; + oldGI_[i] = 0.0; + currentGII_[i] = 0.0; + oldGII_[i] = 0.0; + } + else if ( (patch().size()==2) && (nNewFaces == 2) ) + { + label i=0; + this->valueFraction()[i] = symmTensor::zero; + traction_[i] = vector::zero; + i=1; + this->valueFraction()[i] = symmTensor::zero; + traction_[i] = vector::zero; + + cracked_[i] = false; + curTractionN_[i] = 0.0; + oldTractionN_[i] = 0.0; + curTractionS_[i] = 0.0; + oldTractionS_[i] = 0.0; + deltaN_[i] = 0.0; + oldDeltaN_[i] = 0.0; + deltaS_[i] = 0.0; + oldDeltaS_[i] = 0.0; + curDeltaN_[i] = 0.0; + curDeltaS_[i] = 0.0; + unloadingDeltaEff_[i] = 0.0; + currentGI_[i] = 0.0; + oldGI_[i] = 0.0; + currentGII_[i] = 0.0; + oldGII_[i] = 0.0; + } + else + { + for (label i = 1; i < patch().size(); i++) + { + if (addressing[i] == 0) + { + this->valueFraction()[i] = symmTensor::zero; + traction_[i] = vector::zero; + + cracked_[i] = false; + curTractionN_[i] = 0.0; + oldTractionN_[i] = 0.0; + curTractionS_[i] = 0.0; + oldTractionS_[i] = 0.0; + deltaN_[i] = 0.0; + oldDeltaN_[i] = 0.0; + deltaS_[i] = 0.0; + oldDeltaS_[i] = 0.0; + curDeltaN_[i] = 0.0; + curDeltaS_[i] = 0.0; + unloadingDeltaEff_[i] = 0.0; + currentGI_[i] = 0.0; + oldGI_[i] = 0.0; + currentGII_[i] = 0.0; + oldGII_[i] = 0.0; } } } @@ -640,10 +649,10 @@ void solidCohesiveFvPatchVectorField::rmap const labelList& addr ) { - directionMixedFvPatchVectorField::rmap(ptf, addr); + directionMixedFvPatchVectorField::rmap(ptf, addr); - const solidCohesiveFvPatchVectorField& dmptf = - refCast(ptf); + const solidCohesiveFvPatchVectorField& dmptf = + refCast(ptf); relaxationFactor_ = dmptf.relaxationFactor_; traction_.rmap(dmptf.traction_, addr); @@ -681,7 +690,8 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() this->db().objectRegistry::lookupObject ( "rheologyProperties" - ); + ); + label patchID = patch().index(); const scalarField sigmaMax = rheology.cohLaw().sigmaMax()().boundaryField()[patchID]; @@ -695,32 +705,35 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // update old values if (curTimeIndex_ != this->db().time().timeIndex()) { - // we force the penalty factor to be calculated here for the first time - // as all processors must call this at the same time - if (contact_ && patch().size() > 0) - { - // force calculation of penalty factor here - penaltyFactor(); - } + // we force the penalty factor to be calculated here for the first time + // as all processors must call this at the same time + if (contact_ && patch().size() > 0) + { + // force calculation of penalty factor here + penaltyFactor(); + } - oldDeltaN_ = deltaN_; - oldDeltaS_ = deltaS_; - // curDelta is always latest delta - // whereas delta is only latest when explicitDist is false - //oldDeltaN_ = curDeltaN_; - //oldDeltaS_ = curDeltaS_; - oldTractionN_ = curTractionN_; - oldTractionS_ = curTractionS_; - oldGI_ = currentGI_; - oldGII_ = currentGII_; - unloadingDeltaEff_ = - max - ( - unloadingDeltaEff_, - Foam::sqrt(max(deltaN_, 0.0)*max(deltaN_, 0.0) + deltaS_*deltaS_) - ); + oldDeltaN_ = deltaN_; + oldDeltaS_ = deltaS_; + // curDelta is always latest delta + // whereas delta is only latest when explicitDist is false + //oldDeltaN_ = curDeltaN_; + //oldDeltaS_ = curDeltaS_; + oldTractionN_ = curTractionN_; + oldTractionS_ = curTractionS_; + oldGI_ = currentGI_; + oldGII_ = currentGII_; + unloadingDeltaEff_ = + max + ( + unloadingDeltaEff_, + Foam::sqrt + ( + max(deltaN_, 0.0)*max(deltaN_, 0.0) + deltaS_*deltaS_ + ) + ); - curTimeIndex_ = this->db().time().timeIndex(); + curTimeIndex_ = this->db().time().timeIndex(); } @@ -741,40 +754,40 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // global patch material field is needed for multimaterial cohesive laws if (updateGlobalPatchMaterials_) - { - updateGlobalPatchMaterials_ = false; + { + updateGlobalPatchMaterials_ = false; - if (mesh.objectRegistry::foundObject("materials")) - { - scalarField localPatchMaterials = - patch().lookupPatchField - ("materials").patchInternalField(); - scalarField newGlobalPatchMaterials = - crackerMesh.globalCrackField(localPatchMaterials); + if (mesh.objectRegistry::foundObject("materials")) + { + scalarField localPatchMaterials = + patch().lookupPatchField + ("materials").patchInternalField(); + scalarField newGlobalPatchMaterials = + crackerMesh.globalCrackField(localPatchMaterials); - globalPatchMaterials_.setSize(newGlobalPatchMaterials.size()); - globalPatchMaterials_ = newGlobalPatchMaterials; - } - } + globalPatchMaterials_.setSize(newGlobalPatchMaterials.size()); + globalPatchMaterials_ = newGlobalPatchMaterials; + } + } const regionSplit& regions = crackerMesh.regions(); labelField faceCellRegion(size(), -1); + forAll(faceCellRegion, faceI) { faceCellRegion[faceI] = regions[faceCells[faceI]]; } + labelField globalFaceCellRegion = crackerMesh.globalCrackField(faceCellRegion); - // Patch displacement vectorField UPatch = *this; if (fieldName_ == "DU") - { - UPatch += - patch().lookupPatchField("U"); - } + { + UPatch += patch().lookupPatchField("U"); + } // Global displacement vectorField globalUPatch = crackerMesh.globalCrackField(UPatch); @@ -796,285 +809,304 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() //globalIndex = crackerMesh.localCrackStart(); for (label i = 0; i < patch().size(); i++) - { - // update deltas - curDeltaN_[i] = n[i] & delta[i]; - curDeltaS_[i] = mag( (I - sqr(n[i])) & delta[i]); // shearing - if (explicitSeparationDistance_) - { - deltaN_[i] = oldDeltaN_[i]; - deltaS_[i] = oldDeltaS_[i]; - - if (mag(deltaN_[i]) < SMALL) - { - deltaN_[i] = 2*SMALL; - } - FatalError - << "explicitSeparationDistance_ is true!" << abort(FatalError); - } - else - { - // relax - deltaN_[i] = - relaxationFactor_*curDeltaN_[i] - + (1.0 - relaxationFactor_)*deltaN_[i]; - deltaS_[i] = - relaxationFactor_*curDeltaS_[i] - + (1.0 - relaxationFactor_)*deltaS_[i]; - } - - // update tractions - curTractionN_[i] = n[i] & traction_[i]; - curTractionS_[i] = mag( (I-sqr(n[i])) & traction_[i] ); - - // effective delta - only positive deltaN is considered - const scalar deltaEff = - ::sqrt(max(deltaN_[i], 0.0)*max(deltaN_[i], 0.0) - + deltaS_[i]*deltaS_[i]); - - // update energies - // stop calculating after cracking for convergence - // because crack might jump in and out of damaged/failed - // only update energies if there is loading - if ( !cracked_[i] && (deltaEff > (unloadingDeltaEff_[i]-SMALL)) ) { - // if the average normal stress is tensile - if ((curTractionN_[i]+oldTractionN_[i]) > 0.0) + // update deltas + curDeltaN_[i] = n[i] & delta[i]; + curDeltaS_[i] = mag( (I - sqr(n[i])) & delta[i]); // shearing + if (explicitSeparationDistance_) { - // trapezoidal rule - currentGI_[i] = - oldGI_[i] - + ((0.5*(curTractionN_[i]+oldTractionN_[i])) - *(deltaN_[i]-oldDeltaN_[i])); + deltaN_[i] = oldDeltaN_[i]; + deltaS_[i] = oldDeltaS_[i]; + + if (mag(deltaN_[i]) < SMALL) + { + deltaN_[i] = 2*SMALL; + } + FatalError + << "explicitSeparationDistance_ is true!" << abort(FatalError); } else { - // no modeI energy lost if in compression - currentGI_[i] = oldGI_[i]; + // relax + deltaN_[i] = + relaxationFactor_*curDeltaN_[i] + + (1.0 - relaxationFactor_)*deltaN_[i]; + + deltaS_[i] = + relaxationFactor_*curDeltaS_[i] + + (1.0 - relaxationFactor_)*deltaS_[i]; } - // mode II - trapezoidal rule - currentGII_[i] = - oldGII_[i] - + ((0.5*(curTractionS_[i]+oldTractionS_[i])) - *(deltaS_[i]-oldDeltaS_[i])); - } - //scalar currentG = currentGI_[i] + currentGII_[i]; - // if - // ( - // ( globalFaceCellRegion[globalIndex] - // == globalFaceCellRegion[gcfa[globalIndex]] ) - // || - // !cracked_[i] - // ) + // update tractions + curTractionN_[i] = n[i] & traction_[i]; + curTractionS_[i] = mag( (I-sqr(n[i])) & traction_[i] ); - // new tractions to be calculated - scalar curNormalTraction = 0; - vector curTangentialTraction = vector::zero; + // effective delta - only positive deltaN is considered + const scalar deltaEff = + Foam::sqrt(max(deltaN_[i], 0.0)*max(deltaN_[i], 0.0) + + deltaS_[i]*deltaS_[i]); - // loading - // if the effective delta is greater than unloadingDeltaEff then - // there is loading - // unloadingDeltaEff is the maximum previsouly reached deltaEff - if (deltaEff > (unloadingDeltaEff_[i]-SMALL) ) - { - // at the moment loading and unloading are the same - // if total energy is dissipated, then fully crack face - //if ( currentG > GIc[i] || cracked_[i] ) //Gc ) - // propagation - if ( ((currentGI_[i]/GIc[i]) + (currentGII_[i]/GIIc[i])) >= 1 ) - { - //Pout << "GIc[i] is " << GIc[i] << ", curG is " << currentG << endl; - if (!cracked_[i]) Pout << "Face " << i << " is fully cracked" << endl; + // update energies + // stop calculating after cracking for convergence + // because crack might jump in and out of damaged/failed + // only update energies if there is loading + if ( !cracked_[i] && (deltaEff > (unloadingDeltaEff_[i]-SMALL)) ) + { + // if the average normal stress is tensile + if ((curTractionN_[i]+oldTractionN_[i]) > 0.0) + { + // trapezoidal rule + currentGI_[i] = + oldGI_[i] + + ((0.5*(curTractionN_[i]+oldTractionN_[i]))* + (deltaN_[i]-oldDeltaN_[i])); + } + else + { + // no modeI energy lost if in compression + currentGI_[i] = oldGI_[i]; + } - cracked_[i] = true; + // mode II - trapezoidal rule + currentGII_[i] = oldGII_[i] + + ((0.5*(curTractionS_[i]+oldTractionS_[i]))* + (deltaS_[i]-oldDeltaS_[i])); + } + //scalar currentG = currentGI_[i] + currentGII_[i]; + + // if + // ( + // ( globalFaceCellRegion[globalIndex] + // == globalFaceCellRegion[gcfa[globalIndex]] ) + // || + // !cracked_[i] + // ) + + // new tractions to be calculated + scalar curNormalTraction = 0; + vector curTangentialTraction = vector::zero; + + // loading + // if the effective delta is greater than unloadingDeltaEff then + // there is loading + // unloadingDeltaEff is the maximum previsouly reached deltaEff + if (deltaEff > (unloadingDeltaEff_[i]-SMALL) ) + { + // at the moment loading and unloading are the same + // if total energy is dissipated, then fully crack face + //if ( currentG > GIc[i] || cracked_[i] ) //Gc ) + // propagation + if ( ((currentGI_[i]/GIc[i]) + (currentGII_[i]/GIIc[i])) >= 1 ) + { + //Pout << "GIc[i] is " << GIc[i] << ", curG is " << currentG << endl; + if (!cracked_[i]) + { + Pout << "Face " << i << " is fully cracked" << endl; + } + + cracked_[i] = true; - // failed faces might come in to contact so we need to deal with them - // here we could use the face delta values but that assumes that the - // opposing crack faces are aligned direclty opposite one another, which - // in general they are not. So we actually need to calculate the actual - // face penetration distances and then a standard penalty approach - // contact is probably fine. We will use simple assumption for now which - // is fine is the relative displacement - // of the opposing faces is not too large - // To-do: we must calculate actual distances - curNormalTraction = 0.0; - curTangentialTraction = vector::zero; - if ( contact_ && deltaN_[i] <= 0.0 ) - { - curNormalTraction = deltaN_[i]*penaltyFactor(); - //Info << "penaltyFactor() is " << penaltyFactor() << endl; + // failed faces might come in to contact so we need to deal with them + // here we could use the face delta values but that assumes that the + // opposing crack faces are aligned direclty opposite one another, which + // in general they are not. So we actually need to calculate the actual + // face penetration distances and then a standard penalty approach + // contact is probably fine. We will use simple assumption for now which + // is fine is the relative displacement + // of the opposing faces is not too large + // To-do: we must calculate actual distances + curNormalTraction = 0.0; + curTangentialTraction = vector::zero; + if ( contact_ && deltaN_[i] <= 0.0 ) + { + curNormalTraction = deltaN_[i]*penaltyFactor(); + //Info << "penaltyFactor() is " << penaltyFactor() << endl; - // friction - vector sDir = (I - sqr(n[i]))&delta[i]; - sDir /= mag(sDir+vector(SMALL,SMALL,SMALL)); - scalar slip = mag(deltaS_[gcfa[i]] - deltaS_[i]); - scalar fricMag = - min + // friction + vector sDir = (I - sqr(n[i]))&delta[i]; + sDir /= mag(sDir+vector(SMALL,SMALL,SMALL)); + scalar slip = mag(deltaS_[gcfa[i]] - deltaS_[i]); + scalar fricMag = + min + ( + slip*penaltyFactor(), + frictionCoeff_*curNormalTraction // stick/slip + ); + curTangentialTraction = fricMag*sDir; + } + + // relax tractions + curNormalTraction = + relaxationFactor_*curNormalTraction + + (1-relaxationFactor_)*(n[i] & traction_[i]); + + curTangentialTraction = + relaxationFactor_*curTangentialTraction + + (1-relaxationFactor_)*( (I -sqr(n[i])) & traction_[i]); + + traction_[i] = curNormalTraction*n[i] + curTangentialTraction; + } + // damging face with positive normal delta + else if ( deltaN_[i] > 0.0 ) + { + if (cracked_[i]) + { + Pout << "Face " << i << " is un-cracked" << endl; + } + + cracked_[i] = false; + + // set traction in a fixed point iteration manner to force + // (tN/sigmaMax)^2 + (tS/tauMax)^2 = 1 + // while also allowing varying mode mix by assuming colinear traction + // and delta + // Dugdale version + // scalar tN = + // sigmaMax[i] * deltaN_[i] / + // (SMALL + ::sqrt( (deltaN_[i]*deltaN_[i]) + // + (deltaS_[i]*deltaS_[i])*(sigmaMax[i]*sigmaMax[i] + // /(tauMax[i]*tauMax[i])) )); + // scalar tS = + // tauMax[i] * deltaS_[i] / + // (SMALL + ::sqrt( (deltaS_[i]*deltaS_[i]) + // + (deltaN_[i]*deltaN_[i])*(tauMax[i]*tauMax[i] + // /(sigmaMax[i]*sigmaMax[i])) )); + + // Given current deltaN and deltaS, the cohesive law + // (Dugdale, linear, etc)gives back current traction: + scalar tN = n[i] & traction_[i]; + scalar tS = mag( (I - sqr(n[i])) & traction_[i] ); + // Update tN and tS + rheology.cohLaw().damageTractions ( - slip*penaltyFactor(), - frictionCoeff_*curNormalTraction // stick/slip - ); - curTangentialTraction = fricMag*sDir; - } - - // relax tractions - curNormalTraction = - relaxationFactor_*curNormalTraction - + (1-relaxationFactor_)*(n[i] & traction_[i]); - curTangentialTraction = - relaxationFactor_*curTangentialTraction - + (1-relaxationFactor_)*( (I -sqr(n[i])) & traction_[i]); - traction_[i] = curNormalTraction*n[i] + curTangentialTraction; - } - - // damging face with positive normal delta - else if ( deltaN_[i] > 0.0 ) - { - if (cracked_[i]) Pout << "Face " << i << " is un-cracked" << endl; - - cracked_[i] = false; - - // set traction in a fixed point iteration manner to force - // (tN/sigmaMax)^2 + (tS/tauMax)^2 = 1 - // while also allowing varying mode mix by assuming colinear traction - // and delta - // Dugdale version - // scalar tN = - // sigmaMax[i] * deltaN_[i] / - // (SMALL + ::sqrt( (deltaN_[i]*deltaN_[i]) - // + (deltaS_[i]*deltaS_[i])*(sigmaMax[i]*sigmaMax[i] - // /(tauMax[i]*tauMax[i])) )); - // scalar tS = - // tauMax[i] * deltaS_[i] / - // (SMALL + ::sqrt( (deltaS_[i]*deltaS_[i]) - // + (deltaN_[i]*deltaN_[i])*(tauMax[i]*tauMax[i] - // /(sigmaMax[i]*sigmaMax[i])) )); - - // Given current deltaN and deltaS, the cohesive law - // (Dugdale, linear, etc)gives back current traction: - scalar tN = n[i] & traction_[i]; - scalar tS = mag( (I - sqr(n[i])) & traction_[i] ); - // Update tN and tS - rheology.cohLaw().damageTractions - ( - tN, - tS, - deltaN_[i], - deltaS_[i], - currentGI_[i], - currentGII_[i], - i, - globalPatchMaterials_ + tN, + tS, + deltaN_[i], + deltaS_[i], + currentGI_[i], + currentGII_[i], + i, + globalPatchMaterials_ ); - // shear delta direction - vector sDir = (I - sqr(n[i]))&delta[i]; - sDir /= mag(sDir+vector(SMALL,SMALL,SMALL)); + // shear delta direction + vector sDir = (I - sqr(n[i]))&delta[i]; + sDir /= mag(sDir+vector(SMALL,SMALL,SMALL)); - // relax tractions - curNormalTraction = - relaxationFactor_*tN - + (1-relaxationFactor_)*(n[i] & traction_[i]); - curTangentialTraction = - relaxationFactor_*(tS*sDir) - + (1-relaxationFactor_)*( (I -sqr(n[i])) & traction_[i]); - traction_[i] = curNormalTraction*n[i] + curTangentialTraction; - } + // relax tractions + curNormalTraction = + relaxationFactor_*tN + + (1-relaxationFactor_)*(n[i] & traction_[i]); - // damaging faces with negative normal delta + curTangentialTraction = + relaxationFactor_*(tS*sDir) + + (1-relaxationFactor_)*( (I -sqr(n[i])) & traction_[i]); + + traction_[i] = curNormalTraction*n[i] + curTangentialTraction; + } + // damaging faces with negative normal delta + else + { + if (cracked_[i]) + { + Pout << "Face " << i << " is un-cracked" << endl; + } + + cracked_[i] = false; + //Pout << "Contact and shearing face " << i << endl; + + //Pout << "Face " << i << " is damaging but in contact" << endl; + // set shear traction from cohesive law + // set normal traction to contact condition + scalar tS = tauMax[i]; + + vector sDir = (I - sqr(n[i])) & delta[i]; + sDir /= mag(sDir+vector(SMALL,SMALL,SMALL)); + + // Simple penalty condition + scalar penaltyFac = 0.0; + if (contact_) + { + penaltyFac = penaltyFactor(); + } + scalar contactTN = deltaN_[i]*penaltyFac; + + // relax tractions + curNormalTraction = + relaxationFactor_*contactTN + + (1-relaxationFactor_)*(n[i] & traction_[i]); + + curTangentialTraction = + relaxationFactor_*(tS*sDir) + + (1-relaxationFactor_)*( (I -sqr(n[i])) & traction_[i]); + + traction_[i] = curNormalTraction*n[i] + curTangentialTraction; + } + } else - { - if (cracked_[i]) Pout << "Face " << i << " is un-cracked" << endl; + { + // unloading + // as the current effective delta is less than the old time + // effective delta - cracked_[i] = false; - //Pout << "Contact and shearing face " << i << endl; + // We have two choice for unloading: + // (a) ductile approach + // unload with initial stiffness (which is infinite) + // similar to ductilve metal stress-strain curve + // as we use infinite intial stiffness, this means to elastic + // energy is recovered + // and we immediately start loading in opposite direction + // (b) brittle approach + // unload straight back to the origin + // this means we recover elastic energy + // + // approach (b) is numerically "nicer"; however, for Dugdale cohesive + // zone this implys that only half the energy is dissipated just before + // failure and then the other half is dissipated at failure. This may + // not be an issue in practive but it does not really make sense. + // Obviously it is fine for a linear cohesive zone as the energy + // is smoothly dissipated up to failure. For now, we will implement + // approach (b), but this requires more thinking... - //Pout << "Face " << i << " is damaging but in contact" << endl; - // set shear traction from cohesive law - // set normal traction to contact condition - scalar tS = tauMax[i]; + // reduce traction linearly with the reduction in delta + scalar scaleFactor = deltaEff/(unloadingDeltaEff_[i]); + traction_[i] = + relaxationFactor_*scaleFactor*traction_[i] + + (1 - relaxationFactor_)*traction_[i]; + } + } - vector sDir = (I - sqr(n[i])) & delta[i]; - sDir /= mag(sDir+vector(SMALL,SMALL,SMALL)); + bool incremental(fieldName_ == "DU"); - // Simple penalty condition - scalar penaltyFac = 0.0; - if (contact_) - { - penaltyFac = penaltyFactor(); - } - scalar contactTN = deltaN_[i]*penaltyFac; - - // relax tractions - curNormalTraction = - relaxationFactor_*contactTN - + (1-relaxationFactor_)*(n[i] & traction_[i]); - curTangentialTraction = - relaxationFactor_*(tS*sDir) - + (1-relaxationFactor_)*( (I -sqr(n[i])) & traction_[i]); - traction_[i] = curNormalTraction*n[i] + curTangentialTraction; - } - } - - else - { - // unloading - // as the current effective delta is less than the old time - // effective delta - - // We have two choice for unloading: - // (a) ductile approach - // unload with initial stiffness (which is infinite) - // similar to ductilve metal stress-strain curve - // as we use infinite intial stiffness, this means to elastic - // energy is recovered - // and we immediately start loading in opposite direction - // (b) brittle approach - // unload straight back to the origin - // this means we recover elastic energy - // - // approach (b) is numerically "nicer"; however, for Dugdale cohesive - // zone this implys that only half the energy is dissipated just before - // failure and then the other half is dissipated at failure. This may - // not be an issue in practive but it does not really make sense. - // Obviously it is fine for a linear cohesive zone as the energy - // is smoothly dissipated up to failure. For now, we will implement - // approach (b), but this requires more thinking... - - // reduce traction linearly with the reduction in delta - scalar scaleFactor = deltaEff/(unloadingDeltaEff_[i]); - traction_[i] = - relaxationFactor_*scaleFactor*traction_[i] - + (1-relaxationFactor_)*traction_[i]; - } - } - - - this->refGrad() = - tractionBoundaryGradient() - ( - traction_, - scalarField(traction_.size(), 0.0), - word(fieldName_), - patch(), - orthotropic_, - nonLinearGeometry::nonLinearNames_[nonLinear_] - )(); + this->refGrad() = tractionBoundaryGradient::snGrad + ( + traction_, + scalarField(traction_.size(), 0.0), + fieldName_, + "U", + patch(), + orthotropic_, + nonLinear_, + incremental + ); directionMixedFvPatchVectorField::updateCoeffs(); } - void solidCohesiveFvPatchVectorField::calcPenaltyFactor() - { + +void solidCohesiveFvPatchVectorField::calcPenaltyFactor() +{ // calculate penalty factor similar to standardPenalty contact model // approx penaltyFactor from mechanical properties // this can then be scaled using the penaltyScale // to-do: write equivalent for orthotropic if (orthotropic_) { - FatalError << "solidCohesiveFvPatchVectorField::calcPenaltyFactor()" + FatalErrorIn + ( + "void solidCohesiveFvPatchVectorField::calcPenaltyFactor()" + ) << "solidCohesiveFvPatchVectorField::calcPenaltyFactor()" << " has yet to be written for orthotropic" << exit(FatalError); } @@ -1083,10 +1115,11 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() const fvMesh& mesh = patch().boundaryMesh().mesh(); const scalarField& mu = - mesh.objectRegistry::lookupObject + mesh.objectRegistry::lookupObject ("mu").boundaryField()[patchID]; + const scalarField& lambda = - mesh.objectRegistry::lookupObject + mesh.objectRegistry::lookupObject ("lambda").boundaryField()[patchID]; // avarage contact patch bulk modulus @@ -1097,13 +1130,16 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // average contact patch cell volume scalar cellVolume = 0.0; + const volScalarField::DimensionedInternalField & V = mesh.V(); { - const unallocLabelList& faceCells = mesh.boundary()[patchID].faceCells(); - forAll(mesh.boundary()[patchID], facei) - { - cellVolume += V[faceCells[facei]]; - } + const unallocLabelList& faceCells = + mesh.boundary()[patchID].faceCells(); + + forAll(mesh.boundary()[patchID], facei) + { + cellVolume += V[faceCells[facei]]; + } } cellVolume /= patch().size(); @@ -1111,7 +1147,7 @@ void solidCohesiveFvPatchVectorField::updateCoeffs() // we approximate penalty factor for traction instead of force penaltyFactorPtr_ = new scalar(penaltyScale_*bulkModulus*faceArea/cellVolume); - } +} // Write diff --git a/src/solidModels/constitutiveModel/constitutiveModel.H b/src/solidModels/constitutiveModel/constitutiveModel.H index f04969abb..731297324 100644 --- a/src/solidModels/constitutiveModel/constitutiveModel.H +++ b/src/solidModels/constitutiveModel/constitutiveModel.H @@ -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 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(); diff --git a/src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.C b/src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.C index d76fdaac8..59a4bf6d7 100644 --- a/src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.C +++ b/src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.C @@ -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 tractionBoundaryGradient::traction +Foam::tmp 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 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 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& mu = - patch.lookupPatchField("mu"); - const fvPatchField& lambda = - patch.lookupPatchField("lambda"); + { + // Lookup material properties from the solver + const fvPatchScalarField& mu = + patch.lookupPatchField("mu"); - // only for nonlinear elastic properties - // if (rheology.type() == plasticityModel::typeName) - // { - // const plasticityModel& plasticity = - // refCast(rheology); + const fvPatchScalarField& lambda = + patch.lookupPatchField("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(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("rheologyProperties"); + + if (rheology.plasticityActive()) + { + traction -= + 2*mu*(n & rheology.DEpsilonP().boundaryField()[patch.index()]); + } + + // Thermal effects + if + ( + patch.boundaryMesh().mesh().objectRegistry:: + foundObject("thermalProperties") + ) + { + const thermalModel& thermo = + patch.boundaryMesh().mesh().objectRegistry:: + lookupObject("thermalProperties"); + + const fvPatchScalarField& threeKalpha = + patch.lookupPatchField + ("((threeK*rho)*alpha)"); + + // Incremental thermal contribution + if (incremental) + { + const fvPatchScalarField& DT = + patch.lookupPatchField("DT"); + + traction -= (n*threeKalpha*DT); + } + else + { + const fvPatchScalarField& T = + patch.lookupPatchField("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 + ( + "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 + ( + "sigma" + ); + traction += (n & sigma); + } - //- if there is plasticity - const constitutiveModel& rheology = - patch.boundaryMesh().mesh().objectRegistry:: - lookupObject("rheologyProperties"); - if (rheology.plasticityActive()) - { - traction -= - 2*mu*(n & rheology.DEpsilonP().boundaryField()[patch.index()]); - } + // Visco-elastic effects + if (rheology.viscoActive()) + { + const fvPatchSymmTensorField& DSigmaCorr = + patch.lookupPatchField + ( + "DSigmaCorr" + ); - //- if there are thermal effects - if (patch.boundaryMesh().mesh().objectRegistry:: - foundObject("thermalProperties")) - { - const thermalModel& thermo = - patch.boundaryMesh().mesh().objectRegistry:: - lookupObject("thermalProperties"); + traction += (n & DSigmaCorr); + } - const fvPatchField& threeKalpha = - patch.lookupPatchField - ("((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& DT = - patch.lookupPatchField("DT"); + if + ( + incremental + && nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN + ) + { + // Incremental total Lagrangian - traction -= (n*threeKalpha*DT); - } - else - { - const fvPatchField& T = - patch.lookupPatchField("T"); + const fvPatchTensorField& gradU = + patch.lookupPatchField + ( + "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& gradU = - patch.lookupPatchField("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& sigma = - patch.lookupPatchField("sigma"); - traction += (n & sigma); - } - - //- visco-elastic - if (rheology.viscoActive()) - { - const fvPatchField& DSigmaCorr = - patch.lookupPatchField("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& gradU = - patch.lookupPatchField("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 tractionBoundaryGradient::operator() - ( - const vectorField& traction, - const scalarField& pressure, - const word fieldName, - const fvPatch& patch, - Switch orthotropic, - word nonLinear - ) const - { - // create tmp + +Foam::tmp 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 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("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("rheologyProperties"); - // required fields - vectorField n = patch.nf(); - const diagTensorField Kinv = inv(K); - const fvPatchField& gradField = - patch.lookupPatchField("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 + ( + "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& sigma = - patch.lookupPatchField("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 + ( + "sigma" + ); - sigmaExp = (n*(n&(C && symm(gradField)))) - (K & gradField); - } - //- updated Lagrangian large strain - else if (nonLinear == "updatedLagrangian") - { - const fvPatchField& sigma = - patch.lookupPatchField("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 + ( + "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& mu = - patch.lookupPatchField("mu"); - const fvPatchField& lambda = - patch.lookupPatchField("lambda"); - - // only for nonlinear elastic properties - // if (rheology.type() == plasticityModel::typeName) - // { - // const plasticityModel& plasticity = - // refCast(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& gradField = - patch.lookupPatchField("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& sigma = - patch.lookupPatchField("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& sigma = - patch.lookupPatchField("sigma"); - - tensorField F = I + gradField; - if (fieldName == "DU" && nonLinear == "totalLagrangian") // for incr TL - { - const fvPatchField& gradU = - patch.lookupPatchField("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("rheologyProperties"); - if (rheology.viscoActive()) - { - //Info << "visco active" << endl; - const fvPatchField& DSigmaCorr = - patch.lookupPatchField("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("thermalProperties")) - { - const thermalModel& thermo = - patch.boundaryMesh().mesh().objectRegistry:: - lookupObject("thermalProperties"); - - const fvPatchField& threeKalpha = - patch.lookupPatchField - ("((threeK*rho)*alpha)"); - - if (fieldName == "DU") - { - const fvPatchField& DT = - patch.lookupPatchField("DT"); - - gradient += (n*threeKalpha*DT); - } + sigmaExp = n*(n &(C && symm(gradField))) - (K & gradField); + } else - { - const fvPatchField& T = - patch.lookupPatchField("T"); + { + FatalErrorIn + ( + "tmp 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("mu"); + + const fvPatchScalarField& lambda = + patch.lookupPatchField("lambda"); + + // only for nonlinear elastic properties + // if (rheology.type() == plasticityModel::typeName) + // { + // const plasticityModel& plasticity = + // refCast(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 + ( + "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 + ( + "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 + ( + "sigma" + ); + + tensorField F = I + gradField; + + if + ( + incremental + && nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN + ) + { + // Incremental total Lagrangian + + const fvPatchTensorField& gradU = + patch.lookupPatchField + ( + "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 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("rheologyProperties"); + + if (rheology.viscoActive()) + { + const fvPatchSymmTensorField& DSigmaCorr = + patch.lookupPatchField + ( + "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("thermalProperties") + ) + { + const thermalModel& thermo = + patch.boundaryMesh().mesh().objectRegistry:: + lookupObject("thermalProperties"); + + const fvPatchScalarField& threeKalpha = + patch.lookupPatchField + ( + "((threeK*rho)*alpha)" + ); + + if (!incremental) + { + const fvPatchScalarField& DT = + patch.lookupPatchField("DT"); + + gradient += n*threeKalpha*DT; + } + else + { + const fvPatchScalarField& T = + patch.lookupPatchField("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& gradU = - patch.lookupPatchField("grad(U)"); + if + ( + incremental + && nonLinear == nonLinearGeometry::TOTAL_LAGRANGIAN + ) + { + // gradU is const in a time step + const fvPatchTensorField& gradU = + patch.lookupPatchField("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 // ************************************************************************* // diff --git a/src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.H b/src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.H index d638e8bd1..a004ebb73 100644 --- a/src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.H +++ b/src/solidModels/constitutiveModel/tractionBoundaryGradient/tractionBoundaryGradient.H @@ -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 traction + static tmp 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 operator() + //- Return surface-normal gradient given traction and pressure + static tmp 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 + ); }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/solidModels/contactModels/frictionContactModels/dirichletNeumannFriction/dirichletNeumannFriction.C b/src/solidModels/contactModels/frictionContactModels/dirichletNeumannFriction/dirichletNeumannFriction.C index 38d5c4a55..c750c72e3 100644 --- a/src/solidModels/contactModels/frictionContactModels/dirichletNeumannFriction/dirichletNeumannFriction.C +++ b/src/solidModels/contactModels/frictionContactModels/dirichletNeumannFriction/dirichletNeumannFriction.C @@ -269,18 +269,24 @@ Foam::dirichletNeumannFriction::dirichletNeumannFriction const fvPatch& slavePatch = mesh.boundary()[slavePatchIndex]; const fvPatchField& gradField = slavePatch.lookupPatchField - ("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 diff --git a/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C b/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C index dd60e5fb8..8ebe07e62 100644 --- a/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C +++ b/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C @@ -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& masterFaceZonePatch, - const PrimitivePatch& 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& masterFaceZonePatch, + const PrimitivePatch& 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 > 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(slaveDisp_); - slaveDisp_ = - localSlaveInterpolator.pointToFaceInterpolate(slaveDispPoints); + { + slaveDispPoints = + localSlaveInterpolator.faceToPointInterpolate(slaveDisp_); - // make sure no tangential component - slaveDisp_ = slaveFaceNormals*(slaveFaceNormals & slaveDisp_); - } + slaveDisp_ = localSlaveInterpolator.pointToFaceInterpolate + (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& gradField = - slavePatch.lookupPatchField - ("grad("+fieldName+")"); - slavePressure_ = tractionBoundaryGradient().traction + const fvPatch& slavePatch = mesh.boundary()[slavePatchIndex]; + const fvPatchField& gradField = + slavePatch.lookupPatchField + ( + "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()); - - 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()); - // 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()); - // 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; } diff --git a/src/solidModels/fvPatchFields/fixedDisplacementOrSolidTraction/fixedDisplacementOrSolidTractionFvPatchVectorField.C b/src/solidModels/fvPatchFields/fixedDisplacementOrSolidTraction/fixedDisplacementOrSolidTractionFvPatchVectorField.C index cbafb9380..43854a1b2 100644 --- a/src/solidModels/fvPatchFields/fixedDisplacementOrSolidTraction/fixedDisplacementOrSolidTractionFvPatchVectorField.C +++ b/src/solidModels/fvPatchFields/fixedDisplacementOrSolidTraction/fixedDisplacementOrSolidTractionFvPatchVectorField.C @@ -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 { diff --git a/src/solidModels/fvPatchFields/fixedDisplacementZeroShear/fixedDisplacementZeroShearFvPatchVectorField.C b/src/solidModels/fvPatchFields/fixedDisplacementZeroShear/fixedDisplacementZeroShearFvPatchVectorField.C index 323443ee1..d48138e28 100644 --- a/src/solidModels/fvPatchFields/fixedDisplacementZeroShear/fixedDisplacementZeroShearFvPatchVectorField.C +++ b/src/solidModels/fvPatchFields/fixedDisplacementZeroShear/fixedDisplacementZeroShearFvPatchVectorField.C @@ -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::operator=(vectorField("value", dict, p.size())); + Field::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 normalValue = transform(valueFraction(), refValue()); + this->refValue() = *this; - Field gradValue = - this->patchInternalField() + refGrad()/this->patch().deltaCoeffs(); + Field normalValue = transform(valueFraction(), refValue()); - Field transformGradValue = - transform(I - valueFraction(), gradValue); + Field gradValue = + this->patchInternalField() + refGrad()/this->patch().deltaCoeffs(); - Field::operator=(normalValue + transformGradValue); + Field transformGradValue = + transform(I - valueFraction(), gradValue); + + Field::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(); } diff --git a/src/solidModels/fvPatchFields/fixedDisplacementZeroShearOrSolidTraction/fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField.C b/src/solidModels/fvPatchFields/fixedDisplacementZeroShearOrSolidTraction/fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField.C index e13489cc9..1d9e096f0 100644 --- a/src/solidModels/fvPatchFields/fixedDisplacementZeroShearOrSolidTraction/fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField.C +++ b/src/solidModels/fvPatchFields/fixedDisplacementZeroShearOrSolidTraction/fixedDisplacementZeroShearOrSolidTractionFvPatchVectorField.C @@ -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") diff --git a/src/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C b/src/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C index fb05c7eae..e9b2d2c9b 100644 --- a/src/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C +++ b/src/solidModels/fvPatchFields/solidContact/solidContactFvPatchVectorField.C @@ -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 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 solidContactFvPatchVectorField::Qc() const { - tmp tQc(new scalarField(patch().size(), 0.0)); - scalarField& Qc = tQc(); + tmp 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& gradField = - patch().lookupPatchField("grad("+fieldName_+")"); - // current Cauchy traction for nonlinear - const vectorField curTrac = - tractionBoundaryGradient().traction - ( - gradField, - fieldName_, - patch(), - orthotropic_, - word(nonLinearGeometry::nonLinearNames_[nonLinear_]) - ); + // average force + const fvPatchField& gradField = + patch().lookupPatchField + ( + "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& gradU = - patch().lookupPatchField("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& oldSigma = - patch().lookupPatchField("sigma_0"); - - vectorField oldSf = patch().Sf(); - if (nonLinear_ != nonLinearGeometry::OFF) - { - // current areas - //tensorField F = I + gradField; - if (nonLinear_ == nonLinearGeometry::TOTAL_LAGRANGIAN) - { - const fvPatchField& gradU = + // current areas + const fvPatchField& gradU = patch().lookupPatchField("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& oldSigma = + patch().lookupPatchField("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& gradU = + patch().lookupPatchField("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; } diff --git a/src/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.C b/src/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.C index 1746fda7f..591902529 100644 --- a/src/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.C +++ b/src/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.C @@ -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::operator= + vectorField::operator= ( this->patchInternalField() + (k&gradField.patchInternalField()) + gradient()/this->patch().deltaCoeffs() ); - fvPatchField::evaluate(); + fvPatchVectorField::evaluate(); } // Write diff --git a/src/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.H b/src/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.H index f939ac5f7..ccb8a64c7 100644 --- a/src/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.H +++ b/src/solidModels/fvPatchFields/solidTraction/solidTractionFvPatchVectorField.H @@ -53,7 +53,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class solidTractionFvPatchVectorField Declaration + Class solidTractionFvPatchVectorField Declaration \*---------------------------------------------------------------------------*/ class solidTractionFvPatchVectorField diff --git a/src/solidModels/fvPatchFields/solidTractionFree/solidTractionFreeFvPatchVectorField.C b/src/solidModels/fvPatchFields/solidTractionFree/solidTractionFreeFvPatchVectorField.C index 4381bd73f..889a5dc05 100644 --- a/src/solidModels/fvPatchFields/solidTractionFree/solidTractionFreeFvPatchVectorField.C +++ b/src/solidModels/fvPatchFields/solidTractionFree/solidTractionFreeFvPatchVectorField.C @@ -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::operator= ( this->patchInternalField() - + (k&gradField.patchInternalField()) + + (k & gradField.patchInternalField()) + gradient()/this->patch().deltaCoeffs() ); diff --git a/src/solidModels/fvPatchFields/solidTractionFree/solidTractionFreeFvPatchVectorField.H b/src/solidModels/fvPatchFields/solidTractionFree/solidTractionFreeFvPatchVectorField.H index 1457588d7..740038b58 100644 --- a/src/solidModels/fvPatchFields/solidTractionFree/solidTractionFreeFvPatchVectorField.H +++ b/src/solidModels/fvPatchFields/solidTractionFree/solidTractionFreeFvPatchVectorField.H @@ -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 diff --git a/src/solidModels/fvPatchFields/timeVaryingFixedDisplacementZeroShear/timeVaryingFixedDisplacementZeroShearFvPatchVectorField.C b/src/solidModels/fvPatchFields/timeVaryingFixedDisplacementZeroShear/timeVaryingFixedDisplacementZeroShearFvPatchVectorField.C index a02015dfb..5b24a90c2 100644 --- a/src/solidModels/fvPatchFields/timeVaryingFixedDisplacementZeroShear/timeVaryingFixedDisplacementZeroShearFvPatchVectorField.C +++ b/src/solidModels/fvPatchFields/timeVaryingFixedDisplacementZeroShear/timeVaryingFixedDisplacementZeroShearFvPatchVectorField.C @@ -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(); diff --git a/src/solidModels/nonLinearGeometry/nonLinearGeometry.H b/src/solidModels/nonLinearGeometry/nonLinearGeometry.H index 7a100f693..fce826279 100644 --- a/src/solidModels/nonLinearGeometry/nonLinearGeometry.H +++ b/src/solidModels/nonLinearGeometry/nonLinearGeometry.H @@ -53,7 +53,7 @@ class nonLinearGeometry { public: - //- non-linear solver options + //- Non-linear solver options enum nonLinearType { OFF, diff --git a/tutorials/solidMechanics/elasticIncrSolidFoam/incrPlateHole/0/DU b/tutorials/solidMechanics/elasticIncrSolidFoam/incrPlateHole/0/DU index 64367517e..a19c42bed 100644 --- a/tutorials/solidMechanics/elasticIncrSolidFoam/incrPlateHole/0/DU +++ b/tutorials/solidMechanics/elasticIncrSolidFoam/incrPlateHole/0/DU @@ -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 diff --git a/tutorials/solidMechanics/elasticIncrSolidFoam/incrPlateHole/constant/polyMesh/boundary b/tutorials/solidMechanics/elasticIncrSolidFoam/incrPlateHole/constant/polyMesh/boundary index 8b4279666..1f798d9a7 100644 --- a/tutorials/solidMechanics/elasticIncrSolidFoam/incrPlateHole/constant/polyMesh/boundary +++ b/tutorials/solidMechanics/elasticIncrSolidFoam/incrPlateHole/constant/polyMesh/boundary @@ -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 | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/solidMechanics/elasticSolidFoam/plateHole/0/U b/tutorials/solidMechanics/elasticSolidFoam/plateHole/0/U index 6dcdc00d9..5f1a3c347 100644 --- a/tutorials/solidMechanics/elasticSolidFoam/plateHole/0/U +++ b/tutorials/solidMechanics/elasticSolidFoam/plateHole/0/U @@ -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; diff --git a/tutorials/solidMechanics/elasticSolidFoam/plateHole/constant/polyMesh/boundary b/tutorials/solidMechanics/elasticSolidFoam/plateHole/constant/polyMesh/boundary index 8b4279666..1f798d9a7 100644 --- a/tutorials/solidMechanics/elasticSolidFoam/plateHole/constant/polyMesh/boundary +++ b/tutorials/solidMechanics/elasticSolidFoam/plateHole/constant/polyMesh/boundary @@ -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 | | \*---------------------------------------------------------------------------*/