Commit dce1d447 authored by Kévin Moreau's avatar Kévin Moreau

update infinite dam notch test case

parent ce676806
MATERIAL_CLASS = elastic_damage
NAME = PMMA
PLANE_STATE = plane_strain
YOUNG_MODULUS = 3500.e6
YOUNG_MODULUS = 3.5e9
POISSON_RATIO = 0.3
Y_CRIT = 700000.
HARDENING = rational
HARDENING_COEFF = 1.
HARDENING_POWER = 0.25
Y_CRIT = 7.e5
HARDENING = cohesive_linear
CZM_WF = 1.e-5
CZM_SIGMAC = 7.e7
LC = 5.e-5
......@@ -5,26 +5,27 @@ disp_space_order = 1
disp_bc_integ_order = 4
mean_field_space_type = PolyLagrange
damage_shape = Poly2Revert
lc = 1.
lc = 5.e-5
delta_phi_max_ratio = 0.99
do_forced_delocalization = 0
do_crack_cut = 1
fit_to_vertices_ratio = -1.e-2
nb_step_max = 3000
nb_step_max = 600
export_manager = {
disp 10
eps_ref 10
strain 10
damage 10
damage_peak 1
sm_nonlocal 10
analytic_disp 1
analytic_strain 1
analytic_stress 1 }
export_sensors_label = { load_factor 1 force_y 1 singularity_exponent 1 }
export_sensors_label = { load_factor 1 force_y 1 singularity_exponent 1 dissipated_energy 1 nonlocal_bb_max_x 1 }
export_sensors_point = { }
do_infinite_mapping = 1
do_post_pro_filtering = 1
infinite_mapping_offset = 0.
notch_angle = 90.
internal_radius = 5.
external_radius = 6.5
internal_radius = 5.e-4
external_radius = 6.5e-4
// can be changed
r=5.;
r=5.e-4;
R=1.3*r;
costheta=0.70710678118;
sintheta=0.70710678118;
Nn=30;
a=1.;
lc=1.;
lc=5.e-5;
lcmin=lc/Nn;
lcmax=10.*lcmin;
distmin=1.5*lc;
......
This diff is collapsed.
// can be changed
r=0.5e-3;
R=1.3*r;
costheta=0.70710678118;
sintheta=0.70710678118;
Nn=10;
a=1.;
lc=5.e-5;
lcmin=lc/Nn;
lcmax=10.*lcmin;
distmin=1.5*lc;
distmax=3.*lc;
l=r/3.;
s=r/5.;
Point(1) = {0, 0, 0, a};
Point(2) = {l*costheta, 0, 0, a};
Point(3) = {r, 0, 0, a};
Point(4) = {R, 0, 0, a};
Point(5) = {0, l*sintheta, 0, a};
Point(6) = {0, r, 0, a};
Point(7) = {0, R, 0, a};
Point(8) = {0, -l*sintheta, 0, a};
Point(9) = {0, -r, 0, a};
Point(10) = {0, -R, 0, a};
Point(11) = {r-s, 0, 0, a};
Point(12) = {l-s, 0, 0, a};
Line(1) = {2, 5};
Line(2) = {5, 6};
Line(3) = {6, 7};
Line(4) = {2, 8};
Line(5) = {8, 9};
Line(6) = {9, 10};
Circle(7) = {9, 1, 3};
Circle(8) = {3, 1, 6};
Circle(9) = {4, 1, 7};
Circle(10) = {10, 1, 4};
Line(11) = {12, 11};
Line Loop(12) = {2, -8, -7, -5, -4, 1};
Plane Surface(13) = {12};
Line Loop(14) = {8, 3, -9, -10, -6, 7};
Plane Surface(15) = {14};
Physical Line(11) = {9, 10};
Physical Line(12) = {7, 8};
Physical Line(13) = {2, 3, 5, 6};
Physical Surface(101) = {13};
Physical Surface(102) = {15};
Field[1]=Attractor;
Field[1].NNodesByEdge=Nn;
Field[1].EdgesList={11};
Field[2]=Threshold;
Field[2].IField=1;
Field[2].LcMin=lcmin;
Field[2].LcMax=lcmax;
Field[2].DistMin=distmin;
Field[2].DistMax=distmax;
Background Field=2;
// can be changed
r=0.5e-3;
R=1.3*r;
rho=1.e-4;
costheta=0.70710678118;
sintheta=0.70710678118;
Nn=10;
a=1.;
lc=5.e-5;
lcmin=lc/Nn;
lcmax=10.*lcmin;
distmin=1.5*lc;
distmax=3.*lc;
l=r/3.;
s=r/5.;
Point(1) = {0, 0, 0, a};
Point(2) = {l*costheta, 0, 0, a};
Point(3) = {r, 0, 0, a};
Point(4) = {R, 0, 0, a};
Point(5) = {0, l*sintheta, 0, a};
Point(6) = {0, r, 0, a};
Point(7) = {0, R, 0, a};
Point(8) = {0, -l*sintheta, 0, a};
Point(9) = {0, -r, 0, a};
Point(10) = {0, -R, 0, a};
Point(11) = {r-s, 0, 0, a};
Point(12) = {l-s, 0, 0, a};
Point(13) = {l*costheta-rho/sintheta, 0, 0, a};
Point(14) = {l*costheta-rho/sintheta+rho*sintheta, rho*costheta, 0, a};
Point(15) = {l*costheta-rho/sintheta+rho*sintheta, -rho*costheta, 0, a};
Line(1) = {14, 5};
Line(2) = {5, 6};
Line(3) = {6, 7};
Line(4) = {15, 8};
Line(5) = {8, 9};
Line(6) = {9, 10};
Circle(7) = {9, 1, 3};
Circle(8) = {3, 1, 6};
Circle(9) = {4, 1, 7};
Circle(10) = {10, 1, 4};
Circle(11) = {15, 13, 14};
Line(12) = {12, 11};
Line Loop(13) = {2, -8, -7, -5, -4, 1, 11};
Plane Surface(14) = {13};
Line Loop(15) = {8, 3, -9, -10, -6, 7};
Plane Surface(16) = {15};
Physical Line(11) = {9, 10};
Physical Line(12) = {7, 8};
Physical Line(13) = {2, 3, 5, 6};
Physical Surface(101) = {14};
Physical Surface(102) = {16};
Field[1]=Attractor;
Field[1].NNodesByEdge=Nn;
Field[1].EdgesList={12};
Field[2]=Threshold;
Field[2].IField=1;
Field[2].LcMin=lcmin;
Field[2].LcMax=lcmax;
Field[2].DistMin=distmin;
Field[2].DistMax=distmax;
Background Field=2;
// can be changed
r=5.;
r=0.25e-3;
R=1.3*r;
rho=r/20.;
rho=1.e-4;
costheta=0.70710678118;
sintheta=0.70710678118;
Nn=30;
Nn=10;
a=1.;
lc=1.;
lc=5.e-5;
lcmin=lc/Nn;
lcmax=10.*lcmin;
distmin=1.5*lc;
distmax=3.*lc;
distmax=2.*lc;
Point(1) = {0, 0, 0, a};
Point(2) = {-r*costheta, r*sintheta, 0, a};
......@@ -22,7 +22,9 @@ Point(7) = {-R*costheta, -R*sintheta, 0, a};
Point(8) = {r-distmax, 0, 0, a};
Point(9) = {-rho*costheta, rho*sintheta, 0, a};
Point(10) = {-rho*costheta, -rho*sintheta, 0, a};
Point(11) = {rho, 0, 0, a};
Point(11) = {-rho/sintheta+rho, 0, 0, a};
Point(12) = {-rho/sintheta, 0, 0, a};
Point(13) = {-rho/sintheta+rho, 0, 0, a};
Line(1) = {2, 9};
Line(2) = {2, 5};
Line(3) = {4, 10};
......@@ -31,10 +33,10 @@ Circle(5) = {2, 1, 3};
Circle(6) = {3, 1, 4};
Circle(7) = {5, 1, 6};
Circle(8) = {6, 1, 7};
Circle(9) = {10, 1, 11};
Circle(10) = {11, 1, 9};
Circle(9) = {9, 12, 13};
Circle(10) = {10, 12, 13};
Line(11) = {11, 8};
Line Loop(10) = {6, 3, 9, 10, -1, 5};
Line Loop(10) = {1, 9, -10, -3, -6, -5};
Plane Surface(11) = {10};
Line Loop(12) = {7, 8, -4, -6, -5, 2};
Plane Surface(13) = {12};
......
// can be changed
r=0.25e-3;
R=1.3*r;
rho=1.e-4;
s=1.e-4;
costheta=0.70710678118;
sintheta=0.70710678118;
Nn=10;
a=1.;
lc=5.e-5;
lcmin=lc/Nn;
lcmax=10.*lcmin;
distmin=1.5*lc;
distmax=2.*lc;
r0=0.;
Point(1) = {0, 0, 0, a};
Point(2) = {r0, 0, 0, a};
Point(3) = {r0-rho+rho/costheta, 0, 0, a};
Point(4) = {r0-rho, 0, 0, a};
Point(5) = {-r*costheta+r0-rho+rho/costheta, r*sintheta, 0, a};
Point(6) = {-R*costheta+r0-rho+rho/costheta, R*sintheta, 0, a};
Point(7) = {-r*costheta+r0-rho+rho/costheta, -r*sintheta, 0, a};
Point(8) = {-R*costheta+r0-rho+rho/costheta, -R*sintheta, 0, a};
Point(9) = {r+r0-rho+rho/costheta, 0, 0, a};
Point(10) = {R+r0-rho+rho/costheta, 0, 0, a};
Point(11) = {r0-rho+rho*costheta, rho*sintheta, 0, a};
Point(12) = {r0-rho+rho*costheta, -rho*sintheta, 0, a};
Point(13) = {r+r0-rho+rho/costheta-s, 0, 0, a};
Circle(1) = {11, 4, 2};
Circle(2) = {2, 4, 12};
Line(3) = {11, 5};
Line(4) = {5, 6};
Line(5) = {12, 7};
Line(6) = {7, 8};
Circle(7) = {7, 3, 9};
Circle(8) = {9, 3, 5};
Circle(9) = {8, 3, 10};
Circle(10) = {10, 3, 6};
Line Loop(11) = {3, -8, -7, -5, -2, -1};
Plane Surface(12) = {11};
Line Loop(13) = {7, 8, 4, -10, -9, -6};
Plane Surface(14) = {13};
Physical Line(11) = {7, 8, 4, 6};
Physical Line(12) = {9, 10};
Physical Surface(101) = {12};
Physical Surface(102) = {14};
Line(13) = {1, 13};
Field[1]=Attractor;
Field[1].NNodesByEdge=Nn;
Field[1].EdgesList={13};
Field[2]=Threshold;
Field[2].IField=1;
Field[2].LcMin=lcmin;
Field[2].LcMax=lcmax;
Field[2].DistMin=distmin;
Field[2].DistMax=distmax;
Background Field=2;
// can be changed
r=0.5e-3;
R=1.3*r;
rho=1.e-4;
costheta=0.70710678118;
sintheta=0.70710678118;
Nn=10;
a=1.;
lc=5.e-5;
lcmin=lc/Nn;
lcmax=10.*lcmin;
distmin=1.5*lc;
distmax=3.*lc;
Point(1) = {0, 0, 0, a};
Point(2) = {-r*costheta, r*sintheta, 0, a};
Point(3) = {r, 0, 0, a};
Point(4) = {-r*costheta, -r*sintheta, 0, a};
Point(5) = {-R*costheta, R*sintheta, 0, a};
Point(6) = {R, 0, 0, a};
Point(7) = {-R*costheta, -R*sintheta, 0, a};
Point(8) = {r-distmax, 0, 0, a};
Point(9) = {-rho/sintheta, 0, 0, a};
Point(10) = {-rho/sintheta+rho*sintheta, rho*costheta, 0, a};
Point(11) = {-rho/sintheta+rho*sintheta, -rho*costheta, 0, a};
Line(1) = {10, 2};
Line(2) = {2, 5};
Line(3) = {4, 11};
Line(4) = {4, 7};
Circle(5) = {2, 1, 3};
Circle(6) = {3, 1, 4};
Circle(7) = {5, 1, 6};
Circle(8) = {6, 1, 7};
Circle(9) = {10, 9, 11};
Line Loop(10) = {1, 5, 6, 3, -9};
Plane Surface(11) = {10};
Line Loop(12) = {7, 8, -4, -6, -5, 2};
Plane Surface(13) = {12};
Line(10) = {9, 8};
Physical Point(1) = {1};
Physical Line(11) = {7, 8};
Physical Line(12) = {5, 6};
Physical Surface(101) = {11};
Physical Surface(102) = {13};
Field[1]=Attractor;
Field[1].NNodesByEdge=Nn;
Field[1].EdgesList={10};
Field[2]=Threshold;
Field[2].IField=1;
Field[2].LcMin=lcmin;
Field[2].LcMax=lcmax;
Field[2].DistMin=distmin;
Field[2].DistMax=distmax;
Background Field=2;
......@@ -10,10 +10,18 @@
// Leguillon
class EvalAsymptoticStress : public xfem::xEval<xfem::xTensor2> {
public:
EvalAsymptoticStress(const double omega, const double alpha, const double coeff) :
omega(omega), alpha(alpha), coeff(coeff), chi(sin(omega*(1.+alpha)/2.)/sin(omega*(1.-alpha)/2.)) {}
EvalAsymptoticStress(const double omega, const double alpha, const double chi, const double coeff, const double nu) :
omega(omega), alpha(alpha), chi(chi), coeff(coeff), nu(nu) {}
void operator()(const xfem::xGeomElem* appro, const xfem::xGeomElem* integ, xfem::xTensor2& res) const {
const Trellis_Util::mPoint& pt=integ->getXYZ();
res(1,1)=1.;
res(2,2)=nu;
return;
Trellis_Util::mPoint pt=integ->getXYZ();
const double omega_leguillon=2.*M_PI-omega;
const double rho=1.e-4;
const double q=(2.*M_PI-omega_leguillon)/M_PI;
const double r0=(1.-1./q)*rho;
pt(0)+=rho/sin(omega_leguillon/2.)-rho+r0;
const double r=sqrt(pt(0)*pt(0)+pt(1)*pt(1));
const double theta=atan2(pt(1), pt(0));
xfem::xTensor2 sig_rtheta;
......@@ -21,12 +29,14 @@ public:
sig_rtheta(1,1)=-cos((1.+alpha)*theta)+(1.+alpha)/(1.-alpha)*chi*cos((1.-alpha)*theta);
sig_rtheta(0,1)=-sin((1.+alpha)*theta)+chi*sin((1.-alpha)*theta);
sig_rtheta(1,0)=sig_rtheta(0,1);
sig_rtheta*=coeff*std::pow(std::max(r, 1.e-3), alpha-1.);
sig_rtheta*=coeff*std::pow(r, alpha-1.);
sig_rtheta(2,2)=-nu*(sig_rtheta(0,0)+sig_rtheta(1,1));
xfem::xTensor2 P;
P(0,0)=cos(theta);
P(0,1)=sin(theta);
P(1,0)=-sin(theta);
P(1,1)=cos(theta);
P(2,2)=1.;
xfem::xTensor2 P_t(P);
P_t.transpose();
res=P_t*sig_rtheta*P;
......@@ -34,19 +44,24 @@ public:
private:
const double omega;
const double alpha;
const double coeff;
const double chi;
const double coeff;
const double nu;
};
class EvalAsymptoticDisp : public xfem::xEval<xfem::xVector> {
public:
EvalAsymptoticDisp(const double omega, const double alpha, const double coeff, const double lambda, const double mu) :
omega(omega), alpha(alpha), coeff(coeff),
chi(sin(omega*(1.+alpha)/2.)/sin(omega*(1.-alpha)/2.)),
EvalAsymptoticDisp(const double omega, const double alpha, const double chi, const double coeff, const double lambda, const double mu) :
omega(omega), alpha(alpha), chi(chi), coeff(coeff),
A((lambda+3.*mu-alpha*(lambda+mu))/((lambda+mu)*(1.-alpha))),
B((lambda+3.*mu+alpha*(lambda+mu))/((lambda+mu)*(1.-alpha))) {}
void operator()(const xfem::xGeomElem* appro, const xfem::xGeomElem* integ, xfem::xVector& res) const {
const Trellis_Util::mPoint& pt=integ->getXYZ();
Trellis_Util::mPoint pt=integ->getXYZ();
const double omega_leguillon=2.*M_PI-omega;
const double rho=1.e-4;
const double q=(2.*M_PI-omega_leguillon)/M_PI;
const double r0=(1.-1./q)*rho;
pt(0)+=rho/sin(omega_leguillon/2.)-rho+r0;
const double r=sqrt(pt(0)*pt(0)+pt(1)*pt(1));
const double theta=atan2(pt(1), pt(0));
xfem::xVector u_rtheta;
......@@ -63,8 +78,8 @@ public:
private:
const double omega;
const double alpha;
const double coeff;
const double chi;
const double coeff;
const double A;
const double B;
};
......@@ -98,6 +113,7 @@ int main(int argc, char* argv[]) {
xParseData parse_data;
parse_data.registerInt("do_infinite_mapping");
parse_data.registerInt("do_post_pro_filtering");
parse_data.registerInt("do_superimposed");
parse_data.registerDouble("infinite_mapping_offset");
parse_data.registerDouble("notch_angle");
parse_data.registerDouble("internal_radius");
......@@ -125,8 +141,8 @@ int main(int argc, char* argv[]) {
const xfem::xTensors& props = *(xfem::xMaterialManagerSingleton::instance().getMaterial(&geo_integ)->getProperties());
const double E=props.scalar("YOUNG_MODULUS");
const double nu=props.scalar("POISSON_RATIO");
const double lambda=E*nu/(1.+nu)/(1.-2.*nu);
const double mu=E/2./(1.+nu);
double lambda=E*nu/(1.+nu)/(1.-2.*nu);
double mu=E/2./(1.+nu);
std::cout<<"omega "<<omega<<" "<<omega*180/M_PI<<std::endl;
double alpha=0.5;
......@@ -143,22 +159,27 @@ int main(int argc, char* argv[]) {
post_pro.exportOnTime("singularity_exponent", 0., alpha);
const double coeff1=1./(((1.+alpha)*sin(omega*(1.+alpha)/2.))/((1.-alpha)*sin(omega*(1.-alpha)/2.))-1.)/std::pow(2*M_PI, 1.-alpha)*1.e9;
const double coeff2=coeff1/(2.*mu*alpha);
double chi=sin(omega*(1.+alpha)/2.)/sin(omega*(1.-alpha)/2.);
double coeff1=1./((1.+alpha)/(1.-alpha)*chi-1.)/std::pow(2.*M_PI, 1.-alpha);
if(parse_data.getDouble("notch_angle")==180) {
chi=0.;
coeff1=-1.;
}
double coeff2=coeff1/(2.*mu*alpha);
EvalAsymptoticDisp eval_disp(omega, alpha, coeff2, lambda, mu);
EvalAsymptoticStress eval_stress(omega, alpha, coeff1);
EvalAsymptoticDisp eval_disp(omega, alpha, chi, coeff2, lambda, mu);
EvalAsymptoticStress eval_stress(omega, alpha, chi, coeff1, nu);
EvalAsymptoticStrain eval_strain(eval_stress, lambda, mu);
auto integ_rule=tls_geom.getIntegRuleBasic(0);
auto begin=tls_geom.begin();
auto end=tls_geom.end();
post_pro.exportOnSpace("analytic_stress", 0, eval_stress, integ_rule, begin, end);
post_pro.exportOnSpace("analytic_strain", 0, eval_strain, integ_rule, begin, end);
post_pro.exportOnSpace("analytic_strain", 0, *eval_strain, integ_rule, begin, end);
post_pro.exportOnSpace("analytic_disp", 0, eval_disp, integ_rule, begin, end);
TLSSolver tls_solver(tls_geom, parse_data, pre_pro, post_pro);
QSFormulationSuperimposed formulation(tls_geom, tls_solver, data, parse_data, pre_pro, post_pro, eval_strain, eval_stress);
QSFormulationSuperimposed formulation(tls_geom, tls_solver, data, parse_data, pre_pro, post_pro, *eval_strain, eval_stress);
Algorithm algorithm(formulation, parse_data);
algorithm.run(pre_pro.getStep());
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment