Closed luchete80 closed 8 months ago
Check if is related to #341
At 7f58e66 was working ok
Here also working ok 069208ae175d83fa4a1e69172f456269356e0994
Eorking here: a9b6d115b065b23a03680663f383be88bb260924
Also working ehere : ba5aff42586c630a1e5a53f0c85f2595b1905e4d
Here is breaking: 1833b2a89119384e46aa89e8db2223c5b4f46515
$ git diff ba5aff4
diff --git a/DEVLOG.md b/DEVLOG.md
index e70def6..a538fab 100644
--- a/DEVLOG.md
+++ b/DEVLOG.md
@@ -463,4 +463,7 @@
- Fixed crash in Johnson Cook damage calc
- Prevented sigma asterisk for damage be NAN at sigma_eq = 0
- Added damage count
- - Added damage output
\ No newline at end of file
+ - Added damage output
+ - Corrected Johnson Cook damage expression
+ - Calculate fracture strain and stresses per particle and not per pair
+ - Corrected Acceleration expression
\ No newline at end of file
diff --git a/Source/Damage.cpp b/Source/Damage.cpp
index 571456d..ba85671 100644
--- a/Source/Damage.cpp
+++ b/Source/Damage.cpp
@@ -20,12 +20,23 @@ inline void Domain::CalcDamage(){
//USED IN JOHNSON COOK DAMAGE
#pragma omp parallel for schedule(static) num_threads(Nproc)
#ifdef __GNUC__
- for (size_t i=0; i< solid_part_count ; i++) //Like in Domain::Move
+ for (size_t i=0; i< solid_part_count ; i++) { //Like in Domain::Move
#else
- for (int i=0; i<Particles.Size(); i++)//Like in Domain::Move
+ for (int i=0; i<solid_part_count; i++) {//Like in Domain::Move
#endif
- Particles[i]->CalculateEquivalentStress();
-
+ if (Particles[i]->delta_pl_strain>0.0) {
+ Particles[i]->CalculateEquivalentStress(); //Set Particles[i]->Sigma_eq
+ if (Particles[i]->Sigma_eq>1.0e-3) {
+ sig_as = 1.0/3.0* (Particles[i]->Sigma(0,0)+Particles[i]->Sigma(1,1)+Particles[i]->Sigma(2,2))/Particles[i]->Sigma_eq; //Stress triaxiality sig_m / sig_eff
+ // Particles[i]->eps_f = Particles[i]->mat->damage->CalcFractureStrain(Particles[i]->eff_strain_rate, sig_as, PP[i]->T);
+ Particles[i]->dam_D += Particles[i]->delta_pl_strain/Particles[i]->mat->damage->CalcFractureStrain(Particles[i]->eff_strain_rate, sig_as, Particles[i]->T);
+ // if (Particles[i]->dam_D > 0.0 && Particles[i]->pl_strain ==0.0)
+ // cout <<"dam " <<Particles[i]->dam_D<<", Pl Strain: "<<Particles[i]->pl_strain<<endl;
+ if (Particles[i]->dam_D > 1.0) Particles[i]->dam_D = 1.0;
+ }
+ }//if incremental pl_strain
+ }
+
#pragma omp parallel for schedule (static) private (PP, sig_as, i) num_threads(Nproc)
#ifdef __GNUC__
for (size_t k=0; k<Nproc;k++)
@@ -97,21 +108,21 @@ inline void Domain::CalcDamage(){
//FIRST WE IMPLEMENT JOHNSON COOK FAILURE AS ISLAM 2017
if (dam_D[k][p] <1.0){ //OR dam_df0[][] == 0.0
- for (i=0;i<2 ;i++){
- //IN CASE OF NO DAMAGE INCLUDED; IT IS NOT CALCULATED
- // cout << "sig_00"<<PP[i]->Sigma(0,0)<<endl;
- if (PP[i]->Sigma_eq>1.0e-3)
- sig_as = 1.0/3.0* (PP[i]->Sigma(0,0)+PP[i]->Sigma(1,1)+PP[i]->Sigma(2,2))/PP[i]->Sigma_eq; //Stress triaxiality sig_m / sig_eff
- else
- sig_as = 0.0;
- // cout << "Fr strain "<<PP[i]->mat->damage->CalcFractureStrain(PP[i]->eff_strain_rate, sig_as, PP[i]->T)<<endl;
- // cout << "str_rate, sig_as, ppi "<< PP[i]->eff_strain_rate<<", " <<sig_as<< ", "<<PP[i]->T<<endl;
- omp_set_lock(&PP[i]->my_lock);
- PP[i]->dam_D += PP[i]->delta_pl_strain/PP[i]->mat->damage->CalcFractureStrain(PP[i]->eff_strain_rate, sig_as, PP[i]->T);
- omp_unset_lock(&PP[i]->my_lock);
- //cout << "PP[i]->dam_D: "<<PP[i]->dam_D<<endl;
- //dam_D[k][p] += ;
- }
+ // for (i=0;i<2 ;i++){
+ // //IN CASE OF NO DAMAGE INCLUDED; IT IS NOT CALCULATED
+ // // cout << "sig_00"<<PP[i]->Sigma(0,0)<<endl;
+ // if (PP[i]->Sigma_eq>1.0e-3) {
+ // sig_as = 1.0/3.0* (PP[i]->Sigma(0,0)+PP[i]->Sigma(1,1)+PP[i]->Sigma(2,2))/PP[i]->Sigma_eq; //Stress triaxiality sig_m / sig_eff
+ // // cout << "Fr strain "<<PP[i]->mat->damage->CalcFractureStrain(PP[i]->eff_strain_rate, sig_as, PP[i]->T)<<endl;
+ // // cout << "str_rate, sig_as, ppi "<< PP[i]->eff_strain_rate<<", " <<sig_as<< ", "<<PP[i]->T<<endl;
+ // // omp_set_lock(&PP[i]->my_lock);
+ // // PP[i]->dam_D += PP[i]->delta_pl_strain/PP[i]->eps_f;
+ // // if (PP[i]->dam_D > 1.0) PP[i]->dam_D = 1.0;
+ // // omp_unset_lock(&PP[i]->my_lock);
+ // }
+ // //cout << "PP[i]->dam_D: "<<PP[i]->dam_D<<endl;
+ // //dam_D[k][p] += ;
+ // }
dam_D[k][p] = (PP[0]->dam_D + PP[1]->dam_D) /2.0;
} else {
dam_D[k][p] = 1.0;
@@ -122,9 +133,9 @@ inline void Domain::CalcDamage(){
}//for pair p
} //proc k
- for (size_t i=0; i< solid_part_count ; i++) //Like in Domain::Move
- if (Particles[i]->dam_D==1.0)
- cout << "DAMAGE 1!"<<endl;
+ // for (size_t i=0; i< solid_part_count ; i++) //Like in Domain::Move
+ // if (Particles[i]->dam_D==1.0)
+ // cout << "DAMAGE 1!"<<endl;
}
diff --git a/Source/InteractionAlt.cpp b/Source/InteractionAlt.cpp
index 2b49a33..13e4493 100644
--- a/Source/InteractionAlt.cpp
+++ b/Source/InteractionAlt.cpp
@@ -13,7 +13,7 @@ inline void Domain::CalcAccel() {
Particle *P1, *P2;
double dam_f = 1.0; //if not damage
- #pragma omp parallel for schedule (static) private (P1,P2) num_threads(Nproc)
+ #pragma omp parallel for schedule (static) private (P1,P2,dam_f) num_threads(Nproc)^M
#ifdef __GNUC__
for (size_t k=0; k<Nproc;k++)
#else
@@ -190,10 +190,10 @@ inline void Domain::CalcAccel() {
// Locking the particle 1 for updating the properties
omp_set_lock(&P1->my_lock);
if (!gradKernelCorr){
- P1->a += mj * temp;
+ P1->a += dam_f*mj * temp;^M
//P1->dDensity += mj * (di/dj) * temp1;
} else{
- P1->a += mj * temp_c[0];
+ P1->a += dam_f *mj * temp_c[0];^M
//P1->dDensity += mj * (di/dj) * temp1_c[0];
}
@@ -202,9 +202,9 @@ inline void Domain::CalcAccel() {
// Locking the particle 2 for updating the properties
omp_set_lock(&P2->my_lock);
if (!gradKernelCorr){
- P2->a -= mi * temp;
+ P2->a -= dam_f * mi * temp; ^M
}else {
- P2->a -= mi * temp_c[1];
+ P2->a -= dam_f * mi * temp_c[1];^M
}
omp_unset_lock(&P2->my_lock);
//#endif
@@ -462,7 +462,7 @@ inline void Domain::RateTensorsReduction(){
inline void Domain::CalcDensInc() {
double dam_f = 1.0; //if not damage
Particle *P1, *P2;
- #pragma omp parallel for schedule (static) private (P1,P2) num_threads(Nproc)
+ #pragma omp parallel for schedule (static) private (P1,P2,dam_f) num_threads(Nproc)^M
#ifdef __GNUC__
for (size_t k=0; k<Nproc;k++)
#else
diff --git a/Source/Material.h b/Source/Material.h
index 9676a1f..7ff8ba9 100644
--- a/Source/Material.h
+++ b/Source/Material.h
@@ -56,7 +56,7 @@ public DamageModel {
}
std::string getDamageCriterion() {return string("JohnsonCook");}
double CalcFractureStrain(const double &str_rate, const double &sig_as,const double &T) { //sig_as is stress triaxiality
- return ( (m_D1 + pow(m_D2, m_D3*sig_as)) * (pow (1.0 + (str_rate/m_eps0), m_D4)) ); //Islam 2017 eq. 35,
+ return ( (m_D1 + m_D2 *exp(m_D3*sig_as)) * pow (1.0 + (str_rate/m_eps0), m_D4) * (1.0 + m_D5 * T)); //Islam 2017 eq. 35, ^M
}
virtual ~JohnsonCookDamage(){}
};
diff --git a/Source/Output.cpp b/Source/Output.cpp
index 0cce6a1..0f3921a 100644
--- a/Source/Output.cpp
+++ b/Source/Output.cpp
@@ -275,6 +275,12 @@ inline void Domain::WriteXDMF (char const * FileKey)
dsname.Printf("ps_en");
H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,ps_en);
+^M
+ if (model_damage){^M
+ dsname.Printf("Damage");^M
+ H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,Damage); ^M
+ ^M
+ }^M
dims[0] = 3*Particles.Size();
dsname.Printf("Displacement");
@@ -282,11 +288,7 @@ inline void Domain::WriteXDMF (char const * FileKey)
dsname.Printf("Contact Force");
H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,ContForce);
- if (model_damage){
- // dsname.Printf("Damage");
- // H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,Damage);
-
- }
+^M
// dsname.Printf("Tg Dir");
// H5LTmake_dataset_float(file_id,dsname.CStr(),1,dims,TgDir);
@@ -478,11 +480,11 @@ inline void Domain::WriteXDMF (char const * FileKey)
oss << " </DataItem>\n";
oss << " </Attribute>\n";
if (model_damage){
- // oss << " <Attribute Name=\"damage\" AttributeType=\"Scalar\" Center=\"Node\">\n";
- // oss << " <DataItem Dimensions=\"" << Particles.Size() << "\" NumberType=\"Float\" Precision=\"10\" Format=\"HDF\">\n";
- // oss << " " << fn.CStr() <<":/damage \n";
- // oss << " </DataItem>\n";
- // oss << " </Attribute>\n";
+ oss << " <Attribute Name=\"Damage\" AttributeType=\"Scalar\" Center=\"Node\">\n";^M
+ oss << " <DataItem Dimensions=\"" << Particles.Size() << "\" NumberType=\"Float\" Precision=\"10\" Format=\"HDF\">\n";^M
+ oss << " " << fn.CStr() <<":/Damage \n";^M
+ oss << " </DataItem>\n";^M
+ oss << " </Attribute>\n"; ^M
}
// oss << " <Attribute Name=\"Tg Dir\" AttributeType=\"Vector\" Center=\"Node\">\n";
// oss << " <DataItem Dimensions=\"" << Particles.Size() << " 3\" NumberType=\"Float\" Precision=\"10\" Format=\"HDF\">\n";
diff --git a/Source/SolverLeapfrog.cpp b/Source/SolverLeapfrog.cpp
index 0884f1b..eb939cf 100644
--- a/Source/SolverLeapfrog.cpp
+++ b/Source/SolverLeapfrog.cpp
@@ -395,7 +395,7 @@ inline void Domain::SolveDiffUpdateLeapFrog (double tf, double dt, double dtOut,
for (int p=0;p<Particles.Size();p++){
if (Particles[p]->dam_D>max_dam)
max_dam = Particles[p]->dam_D;
- if (Particles[p]->dam_D>1.0)
+ if (Particles[p]->dam_D==1.0) ^M
dam_count++;
}
Total CPU time: 6.56079 Calculation Times Accel: 0.26%, Density: 0.17%, Stress: 0.27%, Energy: 0.085%, Contact: 0%, Nb: 0.0072%, Update: 0.00046%, Output No. 3 at 0.0001 has been generated Current Time Step = 2.3e-06 Max plastic strain: 0in particle10764 Max Displacements (No Cont Surf): 3 [ 1.1549e-061.1549e-069.36213e-06 ] Int Energy: 0.51654, Kin Energy: -0.0618576
This could be due to material modifications due to damage