added different bond laws to linkslip

src/sm/Elements/bondlink3d.C

@@ -85,9 +85,9 @@ BondLink3d :: computeBmatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer, int
  }
 
  //Assemble Bmatrix based on three rigid arm components
- //first node: beam, second node: lattice
+ //first node: beam, second node: solid
  //rigid.at(1) (tangential), rigid.at(2) (lateral), rigid.at(3) (lateral)
- answer.resize(6, 9);
+ answer.resize(3, 9);
  answer.zero();
 
  //Normal displacement jump in x-direction
@@ -115,8 +115,6 @@ BondLink3d :: computeBmatrixAt(GaussPoint *aGaussPoint, FloatMatrix &answer, int
  //Second node
  answer.at(3, 9) = -1.;
 
- answer.at(4,4) = 1.;
-
  return;
 }
 

src/sm/Materials/InterfaceMaterials/bondceb.C

@@ -148,14 +148,15 @@ BondCEBMaterial :: initializeFrom(InputRecord *ir)
  IR_GIVE_OPTIONAL_FIELD(ir, alpha, _IFT_BondCEBMaterial_al);
 
  // dependent parameter
- s0 = pow(pow(s1, -alpha)*taumax/ks, 1./(1.-alpha));
- if ( s0 > s1 ) {
- s0 = s1;
- ks = taumax/s1;
- OOFEM_WARNING("Parameter ks adjusted");
- }
-
- return StructuralInterfaceMaterial :: initializeFrom(ir);
+ //@todo: why not simply if ks <taumax/s1 choose ks = taumax/s1
+ //s0 = pow(pow(s1, -alpha)*taumax/ks, 1./(1.-alpha));
+ if ( s0 > s1 ) {
+ s0 = s1;
+ ks = taumax/s1;
+ OOFEM_WARNING("Parameter ks adjusted");
+ }
+
+ return StructuralInterfaceMaterial :: initializeFrom(ir);
 }
 
 

src/sm/Materials/linkslip.C

@@ -57,11 +57,8 @@ namespace oofem {
  REGISTER_Material(LinkSlip);
 
  /// constructor which creates a dummy material without a status and without random extension interface
- LinkSlip :: LinkSlip(int n, Domain *d, double e0, double a1, double a2) : LinearElasticMaterial(n, d)
+ LinkSlip :: LinkSlip(int n, Domain *d) : LinearElasticMaterial(n, d)
  {
- eNormalMean = e0;
- alphaOne = a1;
- alphaTwo = a2;
  }
 
 
@@ -88,23 +85,40 @@ namespace oofem {
  {
  IRResultType result; // Required by IR_GIVE_FIELD macro
 
-
  LinearElasticMaterial :: initializeFrom(ir);
 
  //axial stiffness
- IR_GIVE_FIELD(ir, eNormalMean, _IFT_LinkSlip_eNormal); // Macro
+ IR_GIVE_FIELD(ir, kNormal, _IFT_LinkSlip_kn); // Macro
 
  //Ratio of lateral to axial stiffness
  alphaOne = 1000.;
- IR_GIVE_OPTIONAL_FIELD(ir, alphaOne, _IFT_LinkSlip_alphaOne); // Macro
+ IR_GIVE_OPTIONAL_FIELD(ir, alphaOne, _IFT_LinkSlip_a1); // Macro
 
- //Ratio of lateral to axial stiffness
- alphaTwo = 0.;
- IR_GIVE_OPTIONAL_FIELD(ir, alphaTwo, _IFT_LinkSlip_alphaTwo); // Macro
+ //Two models available
+ //0 default. Linear elastic-perfect plastic
+ //1 bond slip according first part of CEB model and then constant
+ //2 bond slip according to CEB model
+ IR_GIVE_OPTIONAL_FIELD(ir, type, _IFT_LinkSlip_type); // Macro
+
+ IR_GIVE_FIELD(ir, tauMax, _IFT_LinkSlip_t0); // Macro
 
+ if(type == 1 || type == 2){
+ IR_GIVE_FIELD(ir, s1, _IFT_LinkSlip_s1);
+ IR_GIVE_FIELD(ir, alpha, _IFT_LinkSlip_alpha);
+ }
 
- //Parameter which limits the stress in slip direction.
- IR_GIVE_FIELD(ir, tauZero, _IFT_LinkSlip_tauZero); // Macro
+ if(type == 2){
+ IR_GIVE_FIELD(ir, s2, _IFT_LinkSlip_s2);
+ IR_GIVE_FIELD(ir, s3, _IFT_LinkSlip_s3);
+ IR_GIVE_FIELD(ir, tauFinal, _IFT_LinkSlip_tf);
+ }
+
+ if(type == 1 || type == 2){
+ if ((type == 1 || type == 2) && this->kNormal < this->tauMax/this->s1) {
+ this->kNormal = this->tauMax/this->s1;
+ OOFEM_WARNING("Parameter kN adjusted");
+ }
+ }
 
  return IRRT_OK;
  }
@@ -118,6 +132,44 @@ namespace oofem {
  return answer;
  }
 
+ double
+ LinkSlip :: evaluateBondStress(const double kappa) const
+ {
+ if(this->type == 0){//elastic-perfectly plastic
+ return this->tauMax;
+ }
+ else if(this->type == 1){//modified bond model 
+ if ( kappa <= 0. ){
+ return 0.;
+ }
+ if ( kappa <= s1 ){
+ return tauMax * pow(kappa/s1, alpha);
+ } 
+ return tauMax;
+ }
+ else if(this->type == 2){
+ if ( kappa <= 0. ){
+ return 0.;
+ }
+ if ( kappa <= s1 ){
+ return tauMax * pow(kappa/s1, alpha);
+ }
+ if ( kappa <= s2 ){
+ return tauMax;
+ }
+ if ( kappa <= s3 ){
+ return tauMax - (tauMax-tauFinal) * (kappa-s2) / (s3-s2);
+ }
+ return tauFinal;
+ }
+ else{//unknown type
+ OOFEM_ERROR("Unknown bond model type. Type should be 0, 1 or 2.");
+ }
+
+ return 0.; //Should not be here. 
+ }
+
+
 
  void
  LinkSlip :: giveRealStressVector_3d(FloatArray &answer,
@@ -132,32 +184,34 @@ namespace oofem {
  FloatArray strainVector;
 
  double tempPlasticStrain = status->givePlasticStrain();
-
- // this->initGpForNewStep(gp);
+ double tempKappa = status->giveKappa();
 
  FloatMatrix stiffnessMatrix;
  this->giveStiffnessMatrix(stiffnessMatrix, ElasticStiffness, gp, atTime);
 
- answer.resize(6);
+ answer.resize(3);
  answer.zero();
 
  /*First component is the slip one for which the stress should be limited using plasiticity (frictional slip between fibre and matrix). The other components are kept elastic. */
  answer.at(1) = (totalStrain.at(1)-tempPlasticStrain)*stiffnessMatrix.at(1, 1);
-
- double f = fabs(answer.at(1))-tauZero;
+
+
+ double f = fabs(answer.at(1)) - evaluateBondStress(tempKappa);
 
  if(f>0){//plastic response.
  //Reduced stress by increasing plastic strain.
  tempPlasticStrain = tempPlasticStrain + sgn(answer.at(1))*f/stiffnessMatrix.at(1, 1);
+ tempKappa = tempKappa + f/stiffnessMatrix.at(1,1);
  answer.at(1) = (totalStrain.at(1)-tempPlasticStrain)*stiffnessMatrix.at(1, 1);
  }
 
  //Compute the final stress components
- for ( int i = 2; i <= 6; i++ ) { // only diagonal terms matter
+ for ( int i = 2; i <= 3; i++ ) { // only diagonal terms matter
  answer.at(i) = stiffnessMatrix.at(i, i) * totalStrain.at(i);
  }
 
  //Set temp values in status needed for dissipation
+ status->letTempKappaBe(tempKappa);
  status->letTempPlasticStrainBe(tempPlasticStrain);
  status->letTempStrainVectorBe(totalStrain);
  status->letTempStressVectorBe(answer);
@@ -178,20 +232,13 @@ namespace oofem {
  {
 
  /* Returns elastic moduli in reduced stress-strain space*/
- answer.resize(6, 6);
+ answer.resize(3, 3);
  answer.zero();
 
  answer.at(1, 1) = 1.;
  answer.at(2, 2) = this->alphaOne; // shear
- answer.at(3, 3) = this->alphaOne; // shear
- answer.at(4, 4) = this->alphaTwo; // shear
- answer.at(5, 5) = this->alphaTwo; // shear
- answer.at(6, 6) = this->alphaTwo; // shear
-
-
- answer.times(this->eNormalMean);
-
-
+ answer.at(3, 3) = this->alphaOne; // shear 
+ answer.times(this->kNormal); 
  }
 
  LinkSlipStatus :: LinkSlipStatus(GaussPoint *g) : StructuralMaterialStatus(g)
@@ -216,14 +263,14 @@ namespace oofem {
  LinkSlip :: giveThermalDilatationVector(FloatArray &answer,
  GaussPoint *gp, TimeStep *tStep)
  //
- // returns a FloatArray(6) of initial strain vector
+ // returns a FloatArray(3) of initial strain vector
  // caused by unit temperature in direction of
  // gp (element) local axes
  //
  {
  double alpha = this->give(tAlpha, gp);
 
- answer.resize(6);
+ answer.resize(3);
  answer.zero();
 
  answer.at(1) = alpha;
@@ -234,37 +281,23 @@ namespace oofem {
  {
  MaterialStatus :: printOutputAt(file, tStep);
 
-
- fprintf(file, " slip ");
- for ( auto &var : strainVector ) {
- fprintf( file, " %.4e", var );
+ fprintf(file, " jump ");
+ for ( auto &val : this->strainVector ) {
+ fprintf(file, " %.4e", val );
  }
 
- fprintf(file, "\n stress");
- for ( auto &var : stressVector ) {
- fprintf( file, " %.4e", var );
+ fprintf(file, "\n traction ");
+ for ( auto &val : this->stressVector ) {
+ fprintf(file, " %.4e", val );
  }
  fprintf(file, "\n");
-
 
- fprintf(file, "plasticStrain %.8e\n", this->plasticStrain);
+ fprintf(file, "plasticStrain %.8e, kappa %.8e\n", this->plasticStrain, this->kappa);
  return;
  }
 
- double
- LinkSlip :: give(int aProperty, GaussPoint *gp)
- {
- if ( aProperty == eNormal_ID ) {
- return 1.;
- } else {
- return LinearElasticMaterial :: give(aProperty, gp);
- }
- }
-
-
-
- contextIOResultType
- LinkSlipStatus :: saveContext(DataStream &stream, ContextMode mode, void *obj)
+ void
+ LinkSlipStatus :: saveContext(DataStream &stream, ContextMode mode)
  //
  // saves full information stored in this Status
  // no temp variables stored
@@ -278,14 +311,17 @@ namespace oofem {
  if ( !stream.write(plasticStrain) ) {
  THROW_CIOERR(CIO_IOERR);
  }
-
-
- return CIO_OK;
+
+ // write a raw data
+ if ( !stream.write(kappa) ) {
+ THROW_CIOERR(CIO_IOERR);
+ }
+
  }
 
 
- contextIOResultType
- LinkSlipStatus :: restoreContext(DataStream &stream, ContextMode mode, void *obj)
+ void
+ LinkSlipStatus :: restoreContext(DataStream &stream, ContextMode mode)
  //
  // restores full information stored in stream to this Status
  //
@@ -296,9 +332,12 @@ namespace oofem {
  if ( !stream.read(plasticStrain) ) {
  THROW_CIOERR(CIO_IOERR);
  }
-
-
- return CIO_OK;
+
+ // read raw data
+ if ( !stream.read(kappa) ) {
+ THROW_CIOERR(CIO_IOERR);
+ }
+
  }
 
  void
@@ -311,20 +350,17 @@ namespace oofem {
  {
  StructuralMaterialStatus :: updateYourself(atTime);
  this->plasticStrain = this->tempPlasticStrain;
+ this->kappa = this->tempKappa;
  }
 
-
-
+
  int
  LinkSlip :: giveIPValue(FloatArray &answer,
  GaussPoint *gp,
  InternalStateType type,
  TimeStep *atTime)
  {
-
  // return LinearElasticMaterial :: giveIPValue(answer, gp, type, atTime);
  return Material :: giveIPValue(answer, gp, type, atTime);
-
-
  }
 }