Line data Source code
1 : #ifndef MODEL_INTERFACE_CRACK_H
2 : #define MODEL_INTERFACE_CRACK_H
3 :
4 : #include <AMReX.H>
5 : #include <AMReX_AmrCore.H>
6 :
7 : #include <iostream>
8 : #include <fstream>
9 :
10 : namespace Model
11 : {
12 : namespace Interface
13 : {
14 : namespace Crack
15 : {
16 : class Crack
17 : {
18 : public:
19 : enum GType {GSQUARE, GMULTIWELL, GPHI4C3, GCUBICM, GSQUAREP, GSQUAREPM, GWULINEAR};
20 : enum G2Type {G2WULINEAR};
21 : enum WType {WSQUARE, WMULTIWELL, WMULTIWELL2, WPHI4C3, WU};
22 :
23 0 : Crack() {};
24 :
25 : // Geometric function involved in representing sharp crack as a band.
26 : // Often this function is given in terms of a variable d = 1-c.
27 : // Here we are implementing w as a function of c.
28 : AMREX_FORCE_INLINE
29 0 : virtual Set::Scalar w_phi(Set::Scalar c, Set::Scalar /*p=0.*/)
30 : {
31 0 : switch(w_type)
32 : {
33 0 : case WSQUARE: return (1.-c)*(1.-c);
34 0 : case WMULTIWELL: return (1.-c)*(1.-c)*c*c;
35 0 : case WMULTIWELL2: return (1.+c)*(1.+c)*(1.-c)*(1.-c);
36 0 : case WPHI4C3: return 1. - 4.*c*c*c + 3.*c*c*c*c;
37 0 : case WU: return 2.0*(1.-c) - (1.-c)*(1.-c);
38 :
39 0 : default: return (1.-c)*(1.-c);
40 : }
41 : }
42 : AMREX_FORCE_INLINE
43 0 : virtual Set::Scalar g_phi(Set::Scalar c, Set::Scalar p=0.)
44 : {
45 0 : switch(g_type)
46 : {
47 0 : case GSQUARE: return c*c;
48 0 : case GMULTIWELL: return (2.-c)*(2.-c)*c*c;
49 0 : case GPHI4C3: return 4.*c*c*c - 3.*c*c*c*c;
50 0 : case GSQUAREP: return std::pow(c,2.*p);
51 0 : case GSQUAREPM: return std::pow(c,2.*(std::pow(p,m_d_exponent)));
52 0 : case GCUBICM: return m_d_exponent*(c*c*c - c*c) + 3.*c*c - 2.*c*c*c;
53 0 : case GWULINEAR: {
54 : // this is for a linear softening curve for the CZM model
55 : // Obtain from Wu 2018 CMAME
56 : // w(c) = c^p / (c^p + Q(c))
57 : // where, Q(c) = a1 (1-c) + a1 a2 (1-c)^2 + a1 a2 a3 (1-c)^3
58 : // for linear law with p=2, a1 = 4 * l_w / (pi * b)
59 : // a2 = -1/2, a3 = 0
60 0 : Set::Scalar d = 1.0-c;
61 0 : return (c*c) / (c*c + ((2.0 * m_pf_czm_a0 * d) - (m_pf_czm_a0 * d * d)));
62 : // if (c < 1.e-5) c = 1.e-5;
63 : // double phi = m_pf_czm_a0 * m_pf_czm_order* std::sqrt(w_phi(c,0.));
64 : // phi *= std::sqrt(1 - std::pow(c, 2*m_pf_czm_order));
65 : // phi /= std::pow(c, m_pf_czm_order+1);
66 : // return 1.0 / (1.0 + phi);
67 : }
68 0 : default: return c*c;
69 : }
70 : }
71 :
72 : // The derivative of the geometric function is often written in terms of w'(d).
73 : // Here, we are taking w'(c). It's better to just take the derivative of the function
74 : // we implemented in w_phi(c) with respect to c.
75 : AMREX_FORCE_INLINE
76 0 : virtual Set::Scalar Dw_phi(Set::Scalar c, Set::Scalar /*p=0.*/)
77 : {
78 0 : switch(w_type)
79 : {
80 0 : case WSQUARE: return -2.*(1.-c);
81 0 : case WMULTIWELL: return 4.*c*c*c - 6.*c*c + 2.*c;
82 0 : case WMULTIWELL2: return 4.*c*c*c - 4.*c;
83 0 : case WPHI4C3: return 12.*(c-1.)*c*c;
84 0 : case WU: return -2.*c;
85 :
86 0 : default: return -2.*(1.-c);
87 : }
88 : }
89 : AMREX_FORCE_INLINE
90 0 : virtual Set::Scalar Dg_phi(Set::Scalar c, Set::Scalar p=0.)
91 : {
92 0 : switch(g_type)
93 : {
94 0 : case GSQUARE: return 2.*c;
95 0 : case GMULTIWELL: return 4.*c*c*c - 12.*c*c + 8.*c;
96 0 : case GPHI4C3: return 12.*(1.-c)*c*c;
97 0 : case GSQUAREP: return 2.*p*std::pow(c,2*p -1.);
98 0 : case GSQUAREPM: return 2.*std::pow(p,m_d_exponent)*(std::pow(c, 2*std::pow(p,m_d_exponent)-1));
99 0 : case GCUBICM: return m_d_exponent*(3.*c*c - 2.*c) + 6.*c - 6.*c*c;
100 0 : case GWULINEAR: {
101 : // Derivative of c^2 / (c^2 + Q(c))
102 : // 2c / (c^2 + Q(c)) - c^2 (2c + Q'(c)) / (c^2+Q(c))^2
103 : // ( 2c (c^2 + Q(c)) - c^2 (2c + Q'(c)) ) / (c^2 + Q(c))^2
104 : // ( 2c Q(c) - c^2 Q'(c) ) / (c^2 + Q(c))^2
105 : // Q(d) = a1 d + a1 a2 d^2 + a1 a2 a3 d^3
106 : // In this case, a = 2*a0, a2 = -0.5, a3 = 0
107 : // Q(d) = 2 a0 d - a0 d^2
108 : // Q'(d) = 2 a0 - 2 a0 d
109 : // Q'(c) - -Q'(d) = -2 a0 + 2 a0 d
110 0 : Set::Scalar d = 1.0 -c;
111 0 : Set::Scalar dQ_c = (-2.0*m_pf_czm_a0) + (2.0*m_pf_czm_a0*d);
112 0 : Set::Scalar Q_c = (2.0*m_pf_czm_a0*d) - (m_pf_czm_a0*d*d);
113 0 : return (2.0*c*Q_c - c*c*dQ_c) / ((c*c + Q_c)*(c*c + Q_c));
114 : }
115 :
116 0 : default: return 2*c;
117 : }
118 : }
119 :
120 : AMREX_FORCE_INLINE
121 0 : virtual Set::Scalar Dg2_phi(Set::Scalar c)
122 : {
123 0 : switch (g2_type)
124 : {
125 0 : case G2WULINEAR: {
126 0 : if (c < 1.e-5) c = 1.e-5;
127 : // Set::Scalar mu = m_pf_czm_a0 * w_phi(c,0.) / (std::pow(c, 2*m_pf_czm_order));
128 0 : Set::Scalar Dmu = m_pf_czm_a0 * ( (2.*m_pf_czm_order*w_phi(c,0.)) + (c * -1. * Dw_phi(c,0.)) ) / (std::pow(c, 1 + 2*m_pf_czm_order));
129 0 : return Dmu * g_phi(c,0.) * g_phi(c,0.);
130 : }
131 :
132 0 : default: Util::Abort(INFO, "Dissipation function type must be specified");
133 : }
134 0 : return NAN;
135 : }
136 :
137 : virtual Set::Scalar Gc(Set::Scalar theta) = 0;
138 : virtual Set::Scalar DGc(Set::Scalar theta) = 0;
139 : virtual Set::Scalar DDGc(Set::Scalar theta) = 0;
140 : virtual Set::Scalar Zeta(Set::Scalar theta) = 0;
141 : virtual Set::Scalar Mobility(Set::Scalar theta) = 0;
142 : virtual Set::Scalar DrivingForceThreshold(Set::Scalar theta) = 0;
143 :
144 : void ExportToFile(std::string filename, amrex::Real dTheta)
145 : {
146 : std::ofstream outFile;
147 : outFile.open(filename);
148 :
149 : for(amrex::Real theta=0; theta<2*pi; theta=theta+dTheta)
150 : {
151 : outFile << theta << " " << Gc(theta) << std::endl;
152 : }
153 : outFile.close();
154 :
155 : }
156 0 : void SetGType(const GType a_type)
157 : {
158 0 : g_type = a_type;
159 0 : }
160 :
161 0 : void SetWType(const WType a_type)
162 : {
163 0 : w_type = a_type;
164 0 : }
165 :
166 : void SetG2Type (const G2Type /*a_type*/)
167 : {
168 :
169 : }
170 :
171 : void SetDuctileExponent(const Set::Scalar m)
172 : {
173 : m_d_exponent = m;
174 : }
175 :
176 0 : void SetPFCZMConstants(const Set::Scalar a_pf_czm_order, const Set::Scalar a_pf_czm_a0)
177 : {
178 0 : m_pf_czm_order = a_pf_czm_order;
179 0 : m_pf_czm_a0 = a_pf_czm_a0;
180 0 : }
181 :
182 : protected:
183 : static constexpr Set::Scalar pi = 3.14159265359;
184 : GType g_type = GType::GSQUARE;
185 : G2Type g2_type = G2Type::G2WULINEAR;
186 : WType w_type = WType::WSQUARE;
187 : Set::Scalar m_d_exponent = 1.;
188 : Set::Scalar m_pf_czm_order = 1.;
189 : Set::Scalar m_pf_czm_a0 = 1.0;
190 : };
191 : }
192 : }
193 : }
194 :
195 : #endif
|
