Line data Source code
1 : #ifndef MODEL_INTERFACE_CRACK_PFCZM_H
2 : #define MODEL_INTERFACE_CRACK_PFCZM_H
3 :
4 : #include <iostream>
5 : #include <fstream>
6 : #include <cmath>
7 :
8 : #include "AMReX.H"
9 : #include "Crack.H"
10 : #include "Set/Set.H"
11 : #include "Util/Util.H"
12 :
13 : namespace Model
14 : {
15 : namespace Interface
16 : {
17 : namespace Crack
18 : {
19 : class PFCZM : public Crack
20 : {
21 : public:
22 : enum FSType {WANG2023, MC};
23 :
24 0 : PFCZM() = default;
25 0 : virtual ~PFCZM() = default;
26 :
27 : // In Wu 2024 formulation, Gc is scaled by a constant C_alpha.
28 : // $C_\alpha = 4 \int_0^1 \sqrt{\alpha(\beta)} d\beta$
29 : // where $\alpha$ is the geometric function defined in the code as w_phi
30 : // For the optimal function chosen in Wu 2024, $C_\alpha = \pi$
31 0 : Set::Scalar Gc (Set::Scalar /*ratio*/) { return _Gc / _c_alpha; }
32 0 : Set::Scalar c_alpha () { return _c_alpha; }
33 0 : Set::Scalar l_w() { return _l_w; }
34 :
35 : // this hould check if we need to do mixed mode loading.
36 0 : bool mixed_mode() { return _mixed_mode; }
37 0 : FSType failure_surface() { return _failure_surface; }
38 :
39 : // added functions for mixed mode - Wang 2023 surface
40 : Set::Scalar GcI_bar () { return _GcI_bar / _c_alpha; }
41 : Set::Scalar GcII_bar () { return _GcII_bar / _c_alpha; }
42 0 : Set::Scalar chi() { return _chi; }
43 0 : Set::Scalar beta_bar() { return _beta_bar; }
44 0 : Set::Scalar sig_t() { return _sig_t; }
45 0 : Set::Scalar tau_s() { return _tau_s; }
46 :
47 : // added functions for mixed mode - mohr columb surface
48 0 : Set::Scalar cohesion() { return _cohesion; }
49 0 : Set::Scalar friction() { return _friction; }
50 :
51 : // For now, we are storing things like youngs modulus stuff here.
52 : // ideally, they should be moved to Material -> Solid -> Affine.
53 0 : Set::Scalar DGc (Set::Scalar /*theta*/) { return 0.0; }
54 0 : Set::Scalar DDGc (Set::Scalar /*theta*/) { return 0.0; }
55 0 : Set::Scalar Zeta(Set::Scalar /*theta*/) { return zeta; }
56 0 : Set::Scalar Mobility (Set::Scalar /*theta*/) {return mobility;}
57 0 : Set::Scalar DrivingForceThreshold (Set::Scalar /*theta*/) {return threshold;}
58 :
59 : private:
60 : Set::Scalar _Gc = NAN;
61 : Set::Scalar _c_alpha = pi; // scaling factor.
62 : Set::Scalar _l_w = 1.e-5; //irwing length
63 : Set::Scalar zeta = NAN;
64 : Set::Scalar mobility = NAN;
65 : Set::Scalar threshold = NAN;
66 :
67 : Set::Scalar czm_a0 = 1.0;
68 : Set::Scalar czm_order = NAN;
69 :
70 : bool _mixed_mode = false;
71 : FSType _failure_surface = FSType::WANG2023;
72 :
73 : Set::Scalar _GcI_bar = 1.e3;
74 : Set::Scalar _GcII_bar = 1.e3;
75 :
76 : // Parameters specific to WANG 2023 failure surface
77 : Set::Scalar _chi = NAN;
78 : Set::Scalar k1 = 1.0;
79 : Set::Scalar _beta_bar = -1.0;
80 : Set::Scalar _sig_t = 1.0;
81 : Set::Scalar _tau_s = 1.0;
82 :
83 : // Parameters specific to MC failure surface
84 : Set::Scalar _cohesion = NAN;
85 : Set::Scalar _friction = NAN;
86 :
87 : public:
88 0 : static void Parse(PFCZM & value, IO::ParmParse & pp)
89 : {
90 0 : pp.query_default("czm_order", value.czm_order, 1.0); //Order 1, 1.5, or 2 (Wu 2024, JMPS)
91 0 : pp.query_default("G_c",value._Gc,1.0e3); // Fracture energy in N/m
92 0 : pp.query_default("zeta",value.zeta,1.e-5); // Lengthscale regularization in m
93 :
94 0 : Set::Scalar fracture_strength = NAN;
95 :
96 0 : Set::Scalar youngs = NAN;
97 0 : Set::Scalar nu = NAN;
98 0 : Set::Scalar tensile_strength = 2.e8;
99 0 : Set::Scalar shear_strength = 1.41e8;
100 0 : Set::Scalar compressive_strength = NAN;
101 :
102 : // Toggle to determine mixed mode
103 0 : pp.query_default("mixed_mode", value._mixed_mode, false);
104 :
105 0 : if (value._mixed_mode)
106 : {
107 : // Young's modulus
108 0 : pp.query_default("E", youngs, 2.e11);
109 :
110 : // Poisson's ratio $\nu$
111 0 : pp.query_default("nu", nu, 0.3);
112 :
113 : // shear modulus
114 0 : Set::Scalar mu = 0.5 * youngs / (1. + nu);
115 :
116 0 : std::string fail_surf_type = "";
117 : // type of failure surface
118 0 : pp.query_validate("failure_surface", fail_surf_type, {"wang2023", "mohr", "mc", "columb"}); // degradation function. For now, we only have linear softening
119 0 : std::map<std::string,FSType> fs_map;
120 0 : fs_map["wang2023"] = FSType::WANG2023;
121 0 : fs_map["mohr"] = FSType::MC;
122 0 : fs_map["mc"] = FSType::MC;
123 0 : fs_map["columb"] = FSType::MC;
124 0 : value._failure_surface = fs_map[fail_surf_type];
125 :
126 0 : if (value._failure_surface == FSType::WANG2023)
127 : {
128 : // Fracture strength
129 0 : pp.query_default("fracture_strength", fracture_strength, 1.e6);
130 :
131 : // $\chi$ value
132 0 : pp.query_default("chi", value._chi, 1.0);
133 :
134 : // Compressive strength
135 0 : pp.query_default("compressive_strength", compressive_strength,4.e8);
136 0 : value.k1 = ( value._chi < 1/std::sqrt(2) ) ? ( 2 * value._chi * std::sqrt(1.0 - value._chi*value._chi) ) : 1.0;
137 :
138 0 : tensile_strength = fracture_strength / value.k1;
139 0 : shear_strength = value._chi * tensile_strength;
140 0 : value._beta_bar = (1.0 / (value._chi * value._chi)) - (compressive_strength * compressive_strength / (4.0 * value._chi * value._chi * shear_strength * shear_strength));
141 :
142 0 : value._GcI_bar = value._Gc/(value.k1*value.k1);
143 0 : value._GcII_bar = value._GcI_bar * youngs * value._chi * value._chi / mu ;
144 0 : fracture_strength = tensile_strength * value.k1;
145 :
146 0 : Util::Message(INFO, "Chi = ", value._chi);
147 0 : Util::Message(INFO, "mu = ", mu);
148 0 : Util::Message(INFO, "k1 = ", value.k1);
149 0 : Util::Message(INFO, "beta_bar = ", value._beta_bar);
150 0 : Util::Message(INFO, "GcI_bar = ", value._GcI_bar);
151 0 : Util::Message(INFO, "GcII_bar = ", value._GcII_bar);
152 0 : Util::Message(INFO, "f_t = ", fracture_strength);
153 0 : Util::Message(INFO, "f_c = ", compressive_strength);
154 0 : Util::Message(INFO, "sig_t = ", tensile_strength);
155 0 : Util::Message(INFO, "tau_s = ", shear_strength);
156 :
157 0 : value._sig_t = tensile_strength;
158 0 : value._tau_s = shear_strength;
159 : }
160 0 : else if (value._failure_surface == FSType::MC)
161 : {
162 : // cohesion coefficient for MC FS type
163 0 : pp.query_default("cohesion", value._cohesion, 1.0);
164 : // friction coefficient for MC FS type
165 0 : pp.query_default("friction", value._friction, 0.5);
166 :
167 0 : Set::Scalar phi = std::atan(value._friction);
168 0 : fracture_strength = 2.0 * value._cohesion * std::cos(phi) / (1.0 + std::sin(phi));
169 : }
170 :
171 0 : }
172 : else
173 : {
174 : // Young's modulus
175 0 : pp.query_default("E", youngs, 2.e11);
176 : // fracture strength beyond which softening starts [Pa]
177 0 : pp.query_default("fracture_strength", fracture_strength,1.e6);
178 : }
179 :
180 0 : Set::Scalar irwing_length = youngs * value._Gc / (fracture_strength * fracture_strength);
181 0 : value.czm_a0 = 4.0 * irwing_length / (pi * value.zeta);
182 0 : Util::Message(INFO, "irwing_length = ", irwing_length, ", a0 = ", value.czm_a0);
183 0 : value._l_w = irwing_length;
184 :
185 0 : if (value.zeta > 0.33 * irwing_length) Util::Warning(INFO, "Zeta value is greater than irwing length. Consider reducing it.");
186 :
187 0 : pp.query_default("mobility",value.mobility,1.0e-2); // Mobility (speed)
188 0 : pp.query_default("threshold", value.threshold,0.0); // Threshold
189 :
190 0 : std::string gtype = "";
191 0 : std::string wtype = "";
192 0 : pp.query_validate("gtype", gtype, {"wu_linear"}); // degradation function. For now, we only have linear softening
193 0 : pp.query_validate("wtype", wtype, {"wu"}); // geometric function. no choice except for the optimal one for now.
194 :
195 0 : std::map<std::string,Model::Interface::Crack::Crack::GType> g_map;
196 0 : g_map["wu_linear"] = Model::Interface::Crack::Crack::GType::GWULINEAR;
197 :
198 0 : std::map<std::string,Model::Interface::Crack::Crack::WType> w_map;
199 0 : w_map["wu"] = Model::Interface::Crack::Crack::WType::WU;
200 :
201 0 : value.SetGType(g_map[gtype]);
202 0 : value.SetWType(w_map[wtype]);
203 :
204 0 : if (w_map[wtype] == Model::Interface::Crack::Crack::WType::WU) value._c_alpha = pi;
205 :
206 0 : value.SetPFCZMConstants(value.czm_order, value.czm_a0);
207 0 : }
208 : };
209 : }
210 : }
211 : }
212 : #endif
|
