Line data Source code
1 : /*************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* */
4 : /* MASTODON */
5 : /* */
6 : /* (c) 2015 Battelle Energy Alliance, LLC */
7 : /* ALL RIGHTS RESERVED */
8 : /* */
9 : /* Prepared by Battelle Energy Alliance, LLC */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* See COPYRIGHT for full restrictions */
13 : /*************************************************/
14 :
15 : // MASTODON includes
16 : #include "ComputeLRIsolatorElasticity.h"
17 :
18 : // MOOSE includes
19 : #include "MooseMesh.h"
20 : #include "Assembly.h"
21 : #include "NonlinearSystem.h"
22 : #include "MooseVariable.h"
23 : #include "MathUtils.h"
24 :
25 : // libmesh includes
26 : #include "libmesh/quadrature.h"
27 : #include "libmesh/utility.h"
28 :
29 : registerMooseObject("MastodonApp", ComputeLRIsolatorElasticity);
30 :
31 : InputParameters
32 66 : ComputeLRIsolatorElasticity::validParams()
33 : {
34 66 : InputParameters params = Material::validParams();
35 66 : params.addClassDescription(
36 : "Compute the forces and the stiffness matrix for an LR isolator element.");
37 : // Switches
38 132 : params.addParam<bool>("cavitation", false, "Switch for modeling cavitation and post-cavitation.");
39 132 : params.addParam<bool>("buckling_load_variation",
40 132 : false,
41 : "Switch for modeling buckling load variation during the analysis.");
42 132 : params.addParam<bool>(
43 : "horizontal_stiffness_variation",
44 132 : false,
45 : "Switch for modeling variation of horizontal stiffness during the analysis.");
46 132 : params.addParam<bool>("vertical_stiffness_variation",
47 132 : false,
48 : "Switch for modeling variation of vertical stiffness during the analysis.");
49 132 : params.addParam<bool>(
50 : "strength_degradation",
51 132 : false,
52 : "Switch for modeling strength degradation due to lead core heating."); // Strength degradation
53 : // due to heating
54 : // Material properties
55 132 : params.addRequiredParam<Real>("fy", "Yield strength of the bearing.");
56 132 : params.addRequiredParam<Real>("alpha",
57 : "Ratio of post-yield shear stiffness to the initial elastic shear "
58 : "stiffness of the bearing. This is dimensionless");
59 132 : params.addRequiredParam<Real>("G_rubber", "Shear modulus of rubber.");
60 132 : params.addRequiredParam<Real>("K_rubber", "Bulk modulus of rubber.");
61 132 : params.addRequiredParam<Real>("D1", "Diameter of lead core.");
62 132 : params.addRequiredParam<Real>("D2", "Outer diameter of the bearing.");
63 132 : params.addRequiredParam<Real>("ts", "Thickness of a single steel shim.");
64 132 : params.addRequiredParam<Real>("tr", "Thickness of a single rubber layer.");
65 132 : params.addRequiredParam<Real>("n", "Number of rubber layers.");
66 132 : params.addRequiredParam<Real>("tc", "Thickness of the rubber cover of the bearing.");
67 132 : params.addRequiredParam<Real>("gamma", "Gamma parameter of Newmark algorithm.");
68 132 : params.addRequiredParam<Real>("beta", "Beta parameter of Newmark algorithm.");
69 132 : params.addParam<Real>("rho_lead", 11200.0, "Density of lead. Defaults to 11200 kg/m3.");
70 132 : params.addParam<Real>(
71 132 : "c_lead", 130.0, "Specific heat capacity of lead. Defaults to 130.0 N-m/kg oC.");
72 132 : params.addParam<Real>(
73 132 : "k_steel", 50.0, "Thermal conductivity of steel. Defaults to 50.0 W/(m-oC).");
74 132 : params.addParam<Real>(
75 132 : "a_steel", 1.41e-05, "Thermal diffusivity of steel. Defaults to 1.41e-05 m2/s.");
76 132 : params.addParam<Real>("kc", 20.0, "Cavitation parameter.");
77 132 : params.addParam<Real>("phi_m", 0.75, "Damage index.");
78 132 : params.addParam<Real>("ac", 1.0, "Strength degradation parameter.");
79 132 : params.addParam<Real>("cd", 0.0, "Viscous damping parameter.");
80 132 : params.set<MooseEnum>("constant_on") = "ELEMENT"; // _qp = 0
81 66 : return params;
82 0 : }
83 :
84 49 : ComputeLRIsolatorElasticity::ComputeLRIsolatorElasticity(const InputParameters & parameters)
85 : : Material(parameters),
86 49 : _cavitation(getParam<bool>("cavitation")),
87 98 : _buckling_load_variation(getParam<bool>("buckling_load_variation")),
88 98 : _horizontal_stiffness_variation(getParam<bool>("horizontal_stiffness_variation")),
89 98 : _vertical_stiffness_variation(getParam<bool>("vertical_stiffness_variation")),
90 98 : _strength_degradation(getParam<bool>("strength_degradation")),
91 98 : _fy(getParam<Real>("fy")),
92 98 : _alpha(getParam<Real>("alpha")),
93 98 : _Gr(getParam<Real>("G_rubber")),
94 98 : _Kr(getParam<Real>("K_rubber")),
95 98 : _d1(getParam<Real>("D1")),
96 98 : _d2(getParam<Real>("D2")),
97 98 : _ts(getParam<Real>("ts")),
98 98 : _tr(getParam<Real>("tr")),
99 98 : _n(getParam<Real>("n")),
100 98 : _tc(getParam<Real>("tc")),
101 98 : _gamma(getParam<Real>("gamma")),
102 98 : _beta(getParam<Real>("beta")),
103 98 : _rhoL(getParam<Real>("rho_lead")), // qL in opensees
104 98 : _cL(getParam<Real>("c_lead")),
105 98 : _kS(getParam<Real>("k_steel")),
106 98 : _aS(getParam<Real>("a_steel")),
107 98 : _kc(getParam<Real>("kc")),
108 98 : _phi_m(getParam<Real>("phi_m")),
109 98 : _ac(getParam<Real>("ac")),
110 98 : _cd(getParam<Real>("cd")),
111 49 : _sD(0.5),
112 98 : _basic_def(getMaterialPropertyByName<ColumnMajorMatrix>("deformations")),
113 98 : _basic_def_old(getMaterialPropertyByName<ColumnMajorMatrix>("old_deformations")),
114 98 : _basic_vel(getMaterialPropertyByName<ColumnMajorMatrix>("deformation_rates")),
115 98 : _basic_vel_old(getMaterialPropertyByName<ColumnMajorMatrix>("old_deformation_rates")),
116 49 : _Fb(declareProperty<ColumnMajorMatrix>("basic_forces")),
117 98 : _Fb_old(getMaterialPropertyOld<ColumnMajorMatrix>("basic_forces")),
118 49 : _Fl(declareProperty<ColumnMajorMatrix>("local_forces")),
119 49 : _Fg(declareProperty<ColumnMajorMatrix>("global_forces")),
120 49 : _Kb(declareProperty<ColumnMajorMatrix>("basic_stiffness_matrix")),
121 49 : _Kl(declareProperty<ColumnMajorMatrix>("local_stiffness_matrix")),
122 49 : _Kg(declareProperty<ColumnMajorMatrix>("global_stiffness_matrix")),
123 98 : _total_gl(getMaterialPropertyByName<ColumnMajorMatrix>("total_global_to_local_transformation")),
124 98 : _total_lb(getMaterialPropertyByName<ColumnMajorMatrix>("total_local_to_basic_transformation")),
125 98 : _length(getMaterialPropertyByName<Real>("initial_isolator_length")),
126 49 : _pi(libMesh::pi),
127 49 : _TL_trial(0.0),
128 49 : _TLC(0.0),
129 49 : _z(declareProperty<ColumnMajorMatrix>("hysteresis_parameter")),
130 98 : _z_old(getMaterialPropertyOld<ColumnMajorMatrix>("hysteresis_parameter")),
131 49 : _umax(declareProperty<ColumnMajorMatrix>("max_tensile_deformation")),
132 98 : _umax_old(getMaterialPropertyOld<ColumnMajorMatrix>("max_tensile_deformation")),
133 49 : _ucn(declareProperty<ColumnMajorMatrix>("initial_cavitation_deformation")),
134 98 : _ucn_old(getMaterialPropertyOld<ColumnMajorMatrix>("initial_cavitation_deformation")),
135 49 : _Fcrmin(declareProperty<ColumnMajorMatrix>("initial_buckling_load")),
136 147 : _Fcrmin_old(getMaterialPropertyOld<ColumnMajorMatrix>("initial_buckling_load"))
137 :
138 : {
139 : // Bearing material and geometric parameters
140 49 : _A = (_pi / 4.0) * ((_d2 + _tc) * (_d2 + _tc) - _d1 * _d1); // Bonded rubber area of bearing
141 49 : Real S = (_d2 * _d2 - _d1 * _d1) / (4 * _d2 * _tr); // Shape factor for case with lead core
142 49 : _Tr = _n * _tr; // height of rubber in the bearing
143 49 : _h = _Tr + (_n - 1) * _ts; // height of rubber + shims
144 : Real F; // Dimension modification factor
145 49 : if (_d1 == 0) // If no lead core, i.e., elastomeric bearing
146 : F = 1.0;
147 : else
148 : {
149 49 : Real r = _d2 / _d1; // Outer to inner diameter ratio
150 49 : F = (r * r + 1) / ((r - 1) * (r - 1)) +
151 49 : (1 + r) / ((1 - r) * log(r)); // Dimension modification factor
152 : }
153 49 : Real Ec = 1.0 / ((1.0 / (6 * _Gr * S * S * F)) +
154 49 : (4.0 / 3.0) * (1.0 / _Kr)); // Compressive modulus of elastomeric bearing
155 49 : Real I = (_pi / 64.0) * (pow((_d2 + _tc), 4) - pow(_d1, 4)); // Moment of inertia of bearing
156 49 : _rg = sqrt(I / _A); // Radius of gyration
157 :
158 : // Bearing shear behavior parameters
159 49 : if (_alpha < 0.0 || _alpha >= 1.0)
160 1 : mooseError("In ComputeLRIsolatorElasticity block, ",
161 : name(),
162 : ". Parameter, '_alpha' should be >= 0.0 and < 1.0.");
163 48 : _qYield0 = _fy * (1 - _alpha);
164 48 : _qYield = _qYield0; // This yield stress changes with time
165 48 : _ke = _Gr * _A / _Tr; // Stiffness of elastic component (due to rubber)
166 48 : _k0 = (1.0 / _alpha - 1.0) * _ke; // Initial stiffness of hysteretic component (due to lead)
167 :
168 : // Axial parameters: compression
169 48 : Real Erot = Ec / 3.0; // Rotation modulus of bearing
170 48 : Real As = _A * _h / _Tr; // Adjusted shear area of bearing
171 48 : Real Is = I * _h / _Tr; // Adjusted moment of intertia of bearing
172 48 : Real Pe = _pi * _pi * Erot * Is / (_h * _h); // Euler buckling load of bearing
173 48 : _kv0 = _A * Ec / _Tr; // Pre-cavitation tensile stiffness at zero lateral displacement
174 48 : _kv = _kv0; // Pre-cavitation stiffness initialized to that at zero displacement
175 48 : _Fcr = -sqrt(Pe * _Gr * As); // Critical buckling load in compression
176 48 : _Fcrn = _Fcr; // Current value of critical buckling load
177 48 : _ucr = _Fcr / _kv0; // Initial value of critical buckling deformation
178 48 : _ucrn = _ucr; // Current value of critical buckling deformation
179 :
180 : // Axial parameters: tension
181 48 : _Fc = 3.0 * _Gr * _A; // Force that initiates cavitation
182 48 : _Fcn = _Fc; // Initial value of cavitation force (will be updated each time step)
183 48 : _uc = _Fc / _kv0; // Deformation at which cavitation is first initiated
184 48 : _Fmax = _Fc; // Initial value of maximum tensile force (will be updated each time step)
185 :
186 48 : _dzdu.reshape(2, 2);
187 : _dzdu.zero();
188 : _ubC.zero();
189 48 : _tC = _t;
190 48 : }
191 :
192 : void
193 44 : ComputeLRIsolatorElasticity::initQpStatefulProperties()
194 : {
195 44 : _z[_qp].reshape(2, 1);
196 44 : _z[_qp].zero();
197 44 : _umax[_qp].reshape(1, 1);
198 44 : _umax[_qp](0) = _uc;
199 44 : _ucn[_qp].reshape(1, 1);
200 44 : _ucn[_qp](0) = _uc;
201 44 : _Fcrmin[_qp].reshape(1, 1);
202 44 : _Fcrmin[_qp](0) = _Fcr;
203 44 : }
204 :
205 : void
206 32108 : ComputeLRIsolatorElasticity::initializeLRIsolator()
207 : {
208 32108 : Real I = (_pi / 64.0) * (pow((_d2 + _tc), 4) - pow(_d1, 4)); // Moment of inertia of bearing
209 32108 : Real Is = I * _h / _Tr; // Adjusted moment of intertia of bearing
210 32108 : Real Er = 3.0 * _Gr; // Elastic modulus of rubber (assuming nu = 0.5)
211 :
212 : // Initializing stiffness matrices
213 32108 : _Kb[_qp].reshape(6, 6);
214 32108 : _Kb[_qp].identity();
215 32108 : _Kb[_qp](0, 0) = _kv0;
216 32108 : _Kb[_qp](1, 1) = _k0 + _ke;
217 32108 : _Kb[_qp](2, 2) = _k0 + _ke;
218 32108 : _Kb[_qp](3, 3) = _Gr * (2 * Is) / _h; // torsional stiffness
219 32108 : _Kb[_qp](4, 4) = Er * Is / _h; // rotational stiffness
220 32108 : _Kb[_qp](5, 5) = Er * Is / _h; // rotational stiffness
221 :
222 : // Initializing forces
223 32108 : _Fb[_qp].reshape(6, 1); // forces in the basic system
224 32108 : _Fb[_qp].zero();
225 32108 : _Fl[_qp].reshape(12, 1); // forces in local system (6+6 = 12dofs)
226 32108 : _Fl[_qp].zero();
227 32108 : _Fg[_qp] = _Fl[_qp]; // forces in global system
228 :
229 32108 : _Fb[_qp] = _Kb[_qp] * _basic_def[_qp];
230 32108 : }
231 :
232 : void
233 32108 : ComputeLRIsolatorElasticity::computeQpProperties()
234 : {
235 32108 : initializeLRIsolator();
236 :
237 : // Compute axial forces and stiffness terms
238 32108 : computeAxial();
239 :
240 : // Computing shear forces and stiffness terms
241 32108 : computeShear();
242 :
243 : // Compute rotational components
244 : // basic x direction
245 32108 : _Fb[_qp](3) = _Kb[_qp](3, 3) * _basic_def[_qp](3);
246 : // basic y direction
247 32108 : _Fb[_qp](4) = _Kb[_qp](4, 4) * _basic_def[_qp](4);
248 : // basic z direction
249 32108 : _Fb[_qp](5) = _Kb[_qp](5, 5) * _basic_def[_qp](5);
250 :
251 : // Finalize forces and stiffness matrix and convert them into global co-ordinate system
252 32108 : finalize();
253 32108 : }
254 :
255 : void
256 32108 : ComputeLRIsolatorElasticity::computeAxial()
257 : {
258 32108 : Real uh = sqrt(_basic_def_old[_qp](1, 0) * _basic_def_old[_qp](1, 0) +
259 32108 : _basic_def_old[_qp](2, 0) * _basic_def_old[_qp](2, 0));
260 :
261 : // Updating Vertical stiffness based on previous step horizontal displacements
262 32108 : if (_vertical_stiffness_variation)
263 : {
264 32108 : _kv = _kv0 * (1.0 / (1.0 + (3.0 / (_pi * _pi)) * (uh / _rg) * (uh / _rg)));
265 32108 : if (uh > 0.0)
266 23380 : _uc = _Fc / _kv;
267 : }
268 :
269 : // Updating current value of critical Buckling load bassed on previous step horizontal
270 : // displacements
271 32108 : if (_buckling_load_variation)
272 : {
273 32108 : if (uh / _d2 > 1.0)
274 0 : mooseError("Horizontal displacement is greater than the outer diameter of the isolator!\n");
275 32108 : Real delta = 2.0 * acos(uh / _d2); // delta is undefined for uh/_d2 > 1.0
276 32108 : Real _Ar = ((_d2 + _tc) * (_d2 + _tc) - _d1 * _d1) / 4.0 *
277 32108 : (delta - sin(delta)); // _Ar does not include lead core
278 32108 : if (_Ar / _A < 0.2 || uh / _d2 >= 1.0)
279 0 : _Fcrn = 0.2 * _Fcr;
280 : else
281 32108 : _Fcrn = _Fcr * _Ar / _A;
282 32108 : if (_Fcrn > _Fcrmin_old[_qp](0))
283 21788 : _Fcrmin[_qp](0) = _Fcrn;
284 : }
285 32108 : _ucrn = _Fcrn / _kv;
286 :
287 : // Calculating force and stiffness in basic x-direction
288 : // Post buckling phase
289 32108 : if (_basic_def[_qp](0) <= _ucrn)
290 : {
291 574 : if (_buckling_load_variation)
292 : {
293 574 : _Kb[_qp](0, 0) = _kv / 1000.0;
294 574 : _Fb[_qp](0, 0) = _Fcrmin[_qp](0) + _Kb[_qp](0, 0) * (_basic_def[_qp](0) - _ucrn);
295 : }
296 : else
297 : {
298 0 : _Kb[_qp](0, 0) = _kv;
299 0 : _Fb[_qp](0, 0) = _Kb[_qp](0, 0) * _basic_def[_qp](0, 0);
300 : }
301 : }
302 :
303 32108 : if (_basic_def[_qp](0, 0) > _ucrn)
304 : {
305 31534 : if (_cavitation) // simulating cavitation effects
306 : {
307 : // Linear Elastic Phase
308 31534 : if (_basic_def[_qp](0, 0) <= _ucn_old[_qp](0))
309 : {
310 29468 : _Kb[_qp](0, 0) = _kv;
311 29468 : _Fb[_qp](0, 0) = _Kb[_qp](0, 0) * _basic_def[_qp](0, 0);
312 29468 : _umax[_qp](0) = _umax_old[_qp](0);
313 29468 : _Fmax = _Fc * (1.0 + (1.0 / (_Tr * _kc)) * (1.0 - exp(-_kc * (_umax[_qp](0) - _uc))));
314 29468 : _Fcn = _Fc * (1.0 - _phi_m * (1.0 - exp(-_ac * (_umax[_qp](0) - _uc) / _uc)));
315 29468 : _ucn[_qp](0) = _Fcn / _kv;
316 : }
317 : // Cavitation un-loading phase
318 2066 : else if (_basic_def[_qp](0, 0) <= _umax_old[_qp](0))
319 : {
320 1630 : _umax[_qp](0) = _umax_old[_qp](0);
321 1630 : _Fmax = _Fc * (1.0 + (1.0 / (_Tr * _kc)) * (1.0 - exp(-_kc * (_umax[_qp](0) - _uc))));
322 1630 : _Fcn = _Fc * (1.0 - _phi_m * (1.0 - exp(-_ac * (_umax[_qp](0) - _uc) / _uc)));
323 1630 : _ucn[_qp](0) = _Fcn / _kv; // cavitation deformation
324 1630 : _Kb[_qp](0, 0) = (_Fmax - _Fcn) / (_umax[_qp](0) - _ucn[_qp](0));
325 3260 : _Fb[_qp](0, 0) = _Fcn + (_Fmax - _Fcn) / (_umax[_qp](0) - _ucn[_qp](0)) *
326 1630 : (_basic_def[_qp](0, 0) - _ucn[_qp](0));
327 : }
328 : // Cavitation reloading phase
329 : else
330 : {
331 436 : _Kb[_qp](0, 0) = (_Fc / _Tr) * exp(-_kc * (_basic_def[_qp](0, 0) - _uc));
332 872 : _Fb[_qp](0, 0) =
333 436 : _Fc * (1.0 + (1.0 / (_Tr * _kc)) * (1.0 - exp(-_kc * (_basic_def[_qp](0, 0) - _uc))));
334 :
335 : // Updating umax, Fmax, Fcn, ucn
336 436 : _umax[_qp](0) = _basic_def[_qp](0, 0);
337 436 : _Fmax = _Fc * (1.0 + (1.0 / (_Tr * _kc)) * (1.0 - exp(-_kc * (_umax[_qp](0) - _uc))));
338 436 : _Fcn = _Fc * (1.0 - _phi_m * (1.0 - exp(-_ac * (_umax[_qp](0) - _uc) / _uc)));
339 436 : _ucn[_qp](0) = _Fcn / _kv;
340 : }
341 : }
342 : else // cavitation effects not simulated
343 : {
344 0 : _Kb[_qp](0, 0) = _kv;
345 0 : _Fb[_qp](0, 0) = _Kb[_qp](0, 0) * _basic_def[_qp](0, 0);
346 : }
347 : }
348 32108 : }
349 :
350 : void
351 32108 : ComputeLRIsolatorElasticity::computeShear()
352 : {
353 : // update horizontal stiffness based on axial load and critical buckling of previous step
354 32108 : if (_horizontal_stiffness_variation)
355 : {
356 32108 : _ke = (_Gr * _A / _Tr) * (1.0 - pow(_Fb_old[_qp](0) / _Fcrn, 2));
357 :
358 32108 : if (_ke < 0)
359 : {
360 572 : _ke = 0.01 * (_Gr * _A / _Tr); // a fraction of _ke to avoid convergence issues
361 : }
362 : }
363 :
364 : // current temperature of the lead core
365 32108 : Real vel1 = _basic_vel_old[_qp](1, 0) +
366 32108 : (_gamma / _beta) * (((_basic_def[_qp](1, 0) - _basic_def_old[_qp](1, 0)) / _dt) -
367 32108 : (_basic_vel_old[_qp](1, 0)));
368 32108 : Real vel2 = _basic_vel_old[_qp](2, 0) +
369 32108 : (_gamma / _beta) * (((_basic_def[_qp](2, 0) - _basic_def_old[_qp](2, 0)) / _dt) -
370 32108 : (_basic_vel_old[_qp](2, 0)));
371 32108 : Real vel = sqrt(vel1 * vel1 + vel2 * vel2);
372 :
373 32108 : _TL_trial = calculateCurrentTemperature(_qYield, _TLC, vel);
374 :
375 : // calculating shear forces and stiffnesses in basic y and z directions
376 : // get displacement increments (trial-committed)
377 :
378 32108 : ColumnMajorMatrix delta_ub = _basic_def[_qp] - _basic_def_old[_qp];
379 :
380 32108 : if (std::sqrt(delta_ub(1) * delta_ub(1) + delta_ub(2) * delta_ub(2)) >= 0.0)
381 : {
382 : // yield displacement
383 32108 : double uy = _qYield / _k0;
384 :
385 : // calculate hysteretic evolution parameter, z, using the Newton-Raphson method
386 32108 : unsigned int iter = 0;
387 : unsigned int maxIter = 100;
388 : Real tol = 1E-9;
389 : Real beta = 0.1; // beta and gamma are as per Nagarajaiah(1991), which is opposite to from Park
390 : // et al.(1986)
391 : Real gamma = 0.9;
392 : double tmp1, tmp2, tmp3;
393 32108 : ColumnMajorMatrix f(2, 1), delta_z(2, 1), Df(2, 2);
394 : do
395 : {
396 42192 : tmp1 = beta + gamma * MathUtils::sign(_z[_qp](0) * delta_ub(1));
397 42192 : tmp2 = beta + gamma * MathUtils::sign(_z[_qp](1) * delta_ub(2));
398 42192 : tmp3 = _z[_qp](0) * delta_ub(1) * tmp1 + _z[_qp](1) * delta_ub(2) * tmp2;
399 :
400 : // function and derivative
401 42192 : f(0) = _z[_qp](0) - _z_old[_qp](0) - 1.0 / uy * (delta_ub(1) - _z[_qp](0) * tmp3);
402 42192 : f(1) = _z[_qp](1) - _z_old[_qp](1) - 1.0 / uy * (delta_ub(2) - _z[_qp](1) * tmp3);
403 :
404 84384 : Df(0, 0) = 1.0 + (1.0 / uy) *
405 42192 : (2 * _z[_qp](0) * delta_ub(1) * tmp1 + _z[_qp](1) * delta_ub(2) * tmp2);
406 42192 : Df(1, 0) = (tmp1 / uy) * _z[_qp](1) * delta_ub(1);
407 42192 : Df(0, 1) = (tmp2 / uy) * _z[_qp](0) * delta_ub(2);
408 84384 : Df(1, 1) = 1.0 + (1.0 / uy) *
409 42192 : (_z[_qp](0) * delta_ub(1) * tmp1 + 2 * _z[_qp](1) * delta_ub(2) * tmp2);
410 :
411 : // issue warning if the diagonal elements of the derivative Df is zero
412 84384 : if (MooseUtils::absoluteFuzzyLessEqual(Df(0, 0), 0.0) ||
413 42192 : MooseUtils::absoluteFuzzyLessEqual(Df(1, 1), 0.0))
414 0 : mooseError("Error in ComputeLRIsolatorElasticity block, ",
415 : name(),
416 : ". Zero Jacobian in Newton-Raphson scheme while solving ",
417 : "for the hysteretic evolution parameter, z.\n");
418 : // advance one step
419 84384 : delta_z(0) =
420 42192 : (f(0) * Df(1, 1) - f(1) * Df(0, 1)) / (Df(0, 0) * Df(1, 1) - Df(0, 1) * Df(1, 0));
421 84384 : delta_z(1) =
422 42192 : (f(0) * Df(1, 0) - f(1) * Df(0, 0)) / (Df(0, 1) * Df(1, 0) - Df(0, 0) * Df(1, 1));
423 42192 : _z[_qp] = _z[_qp] - delta_z;
424 42192 : iter++;
425 74300 : } while ((delta_z.norm() >= tol) && (iter < maxIter));
426 :
427 : // Error if Newton-Raphson scheme did not converge
428 32108 : if (iter >= maxIter)
429 0 : mooseError("Error in block, ",
430 : name(),
431 : ". Could not solve for hysteresis",
432 : " evolution parameter, z, after ",
433 : iter,
434 : " iterations and",
435 : " achieving a norm of ",
436 0 : delta_z.norm(),
437 : ".\n");
438 :
439 : // calculate derivative of hysteretic evolution parameter
440 : Real du1du2, du2du1;
441 32108 : if (delta_ub(1) * delta_ub(2) == 0)
442 : {
443 : du1du2 = 0.0;
444 : du2du1 = 0.0;
445 : }
446 : else
447 : {
448 1502 : du1du2 = delta_ub(1) / delta_ub(2);
449 1502 : du2du1 = delta_ub(2) / delta_ub(1);
450 : }
451 64216 : _dzdu(0, 0) =
452 32108 : (1.0 / uy) * (1.0 - _z[_qp](0) * (_z[_qp](0) * tmp1 + _z[_qp](1) * tmp2 * du2du1));
453 64216 : _dzdu(0, 1) =
454 32108 : (1.0 / uy) * (du1du2 - _z[_qp](0) * (_z[_qp](0) * tmp1 * du1du2 + _z[_qp](1) * tmp2));
455 64216 : _dzdu(1, 0) =
456 32108 : (1.0 / uy) * (du2du1 - _z[_qp](1) * (_z[_qp](0) * tmp1 + _z[_qp](1) * tmp2 * du2du1));
457 64216 : _dzdu(1, 1) =
458 32108 : (1.0 / uy) * (1.0 - _z[_qp](1) * (_z[_qp](0) * tmp1 * du1du2 + _z[_qp](1) * tmp2));
459 :
460 : Real dt = _t - _tC;
461 : if (dt == 0)
462 : dt = _dt;
463 :
464 : // set shear force
465 32108 : _Fb[_qp](1, 0) = _cd * vel1 + _qYield * _z[_qp](0) + _ke * _basic_def[_qp](1, 0);
466 32108 : _Fb[_qp](2, 0) = _cd * vel2 + _qYield * _z[_qp](1) + _ke * _basic_def[_qp](2, 0);
467 :
468 : // set tangent stiffness
469 32108 : _Kb[_qp](1, 1) = (_gamma / _beta) * _cd / _dt + _qYield * _dzdu(0, 0) + _ke;
470 32108 : _Kb[_qp](1, 2) = _qYield * _dzdu(0, 1);
471 32108 : _Kb[_qp](2, 1) = _qYield * _dzdu(1, 0);
472 32108 : _Kb[_qp](2, 2) = (_gamma / _beta) * _cd / _dt + _qYield * _dzdu(1, 1) + _ke;
473 :
474 32108 : if (_Kb[_qp](1, 2) * _Kb[_qp](2, 1) - _Kb[_qp](1, 1) * _Kb[_qp](2, 2) == 0)
475 0 : mooseError("Error in block, ",
476 : name(),
477 : ". Jacobian for isolator is zero due to off diagonal shear elements. Please check "
478 : "again.\n");
479 : }
480 32108 : }
481 :
482 : void
483 32108 : ComputeLRIsolatorElasticity::addPDeltaEffects()
484 : {
485 : // add geometric stiffness to local stiffness
486 32108 : _Kl[_qp](5, 1) -= 0.5 * _Fb[_qp](0, 0);
487 32108 : _Kl[_qp](5, 7) += 0.5 * _Fb[_qp](0, 0);
488 32108 : _Kl[_qp](11, 1) -= 0.5 * _Fb[_qp](0, 0);
489 32108 : _Kl[_qp](11, 7) += 0.5 * _Fb[_qp](0, 0);
490 32108 : _Kl[_qp](4, 2) += 0.5 * _Fb[_qp](0, 0);
491 32108 : _Kl[_qp](4, 8) -= 0.5 * _Fb[_qp](0, 0);
492 32108 : _Kl[_qp](10, 2) += 0.5 * _Fb[_qp](0, 0);
493 32108 : _Kl[_qp](10, 8) -= 0.5 * _Fb[_qp](0, 0);
494 32108 : _Kl[_qp](5, 5) += 0.5 * _Fb[_qp](0, 0) * _sD * _length[_qp];
495 32108 : _Kl[_qp](11, 5) -= 0.5 * _Fb[_qp](0, 0) * _sD * _length[_qp];
496 32108 : _Kl[_qp](4, 4) += 0.5 * _Fb[_qp](0, 0) * _sD * _length[_qp];
497 32108 : _Kl[_qp](10, 4) -= 0.5 * _Fb[_qp](0, 0) * _sD * _length[_qp];
498 32108 : _Kl[_qp](5, 11) -= 0.5 * _Fb[_qp](0, 0) * (1.0 - _sD) * _length[_qp];
499 32108 : _Kl[_qp](11, 11) += 0.5 * _Fb[_qp](0, 0) * (1.0 - _sD) * _length[_qp];
500 32108 : _Kl[_qp](4, 10) -= 0.5 * _Fb[_qp](0, 0) * (1.0 - _sD) * _length[_qp];
501 32108 : _Kl[_qp](10, 10) += 0.5 * _Fb[_qp](0, 0) * (1.0 - _sD) * _length[_qp];
502 32108 : }
503 :
504 : Real
505 32108 : ComputeLRIsolatorElasticity::calculateCurrentTemperature(Real _qYield, Real _TLC, Real vel)
506 : {
507 : // lead core heating
508 32108 : if (_t <= 0)
509 : return 0.0;
510 :
511 32000 : if (_t < _tC)
512 0 : _tC = 0.0;
513 32000 : Real a = _d1 / 2.0;
514 32000 : Real dt = _t - _tC;
515 : // if (dt>1) dt = 0;
516 32000 : Real a_lead = _pi * pow(a, 2);
517 32000 : Real tau = (_aS * _t) / (pow(a, 2));
518 : Real F;
519 32000 : if (tau < 0.6)
520 32000 : F = 2.0 * sqrt(tau / _pi) -
521 32000 : (tau / _pi) * (2.0 - (tau / 4.0) - pow(tau / 4.0, 2) - (15.0 / 4.0) * (pow(tau / 4.0, 3)));
522 : else
523 0 : F = 8.0 / (3.0 * _pi) - (1.0 / (2.0 * sqrt(_pi * tau))) *
524 0 : (1.0 - (1.0 / (12.0 * tau)) + (1.0 / (6.0 * pow(4.0 * tau, 2))) -
525 0 : (1.0 / (12.0 * pow(4.0 * tau, 3))));
526 : Real deltaT1 =
527 32000 : (dt / (_rhoL * _cL * _h)) *
528 32000 : ((_qYield * vel) / a_lead -
529 32000 : (_kS * _TLC / a) * (1.0 / F + 1.274 * ((_n - 1) * _ts / a) * pow(tau, -1.0 / 3.0)));
530 32000 : if (deltaT1 <= 0.0)
531 : deltaT1 = 0.0;
532 :
533 32000 : Real TL_trial1 = _TLC + deltaT1;
534 32000 : Real tCurrent2 = _t + dt;
535 32000 : tau = (_aS * tCurrent2) / (pow(a, 2));
536 32000 : if (tau < 0.6)
537 32000 : F = 2.0 * sqrt(tau / _pi) -
538 32000 : (tau / _pi) * (2.0 - (tau / 4.0) - pow(tau / 4.0, 2) - (15.0 / 4.0) * (pow(tau / 4.0, 3)));
539 : else
540 0 : F = 8.0 / (3.0 * _pi) - (1.0 / (2.0 * sqrt(_pi * tau))) *
541 0 : (1.0 - (1.0 / (12.0 * tau)) + (1.0 / (6.0 * pow(4.0 * tau, 2))) -
542 0 : (1.0 / (12.0 * pow(4.0 * tau, 3))));
543 : Real deltaT2 =
544 : (dt / (_rhoL * _cL * _h)) *
545 32000 : ((_qYield * vel) / a_lead -
546 32000 : (_kS * TL_trial1 / a) * (1.0 / F + 1.274 * ((_n - 1) * _ts / a) * pow(tau, -1.0 / 3.0)));
547 32000 : if (deltaT2 <= 0.0)
548 : deltaT2 = 0.0;
549 :
550 32000 : Real _TL_trial = _TLC + 0.5 * (deltaT1 + deltaT2);
551 :
552 32000 : return _TL_trial;
553 : }
554 :
555 : void
556 32108 : ComputeLRIsolatorElasticity::finalize()
557 :
558 : {
559 : // update lead core heating parameters
560 32108 : _TLC = _TL_trial;
561 32108 : _tC = _t;
562 32108 : if (_strength_degradation)
563 32108 : _qYield = _qYield0 * exp(-0.0069 * _TLC);
564 :
565 : // Converting forces from basic to local to global
566 32108 : _Fl[_qp] = _total_lb[_qp].transpose() * _Fb[_qp]; // local forces
567 32108 : _Fg[_qp] = _total_gl[_qp].transpose() * _Fl[_qp]; // global forces
568 :
569 : // Converting stiffness matrix from basic to local
570 32108 : _Kl[_qp] = _total_lb[_qp].transpose() * _Kb[_qp] * _total_lb[_qp];
571 :
572 : // add P-∆ effects to local stiffness
573 32108 : addPDeltaEffects();
574 :
575 : // Converting stiffness matrix from loacl to global
576 32108 : _Kg[_qp] = _total_gl[_qp].transpose() * _Kl[_qp] * _total_gl[_qp];
577 32108 : }
|
