Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "ComputeIncrementalBeamStrain.h"
11 : #include "MooseMesh.h"
12 : #include "Assembly.h"
13 : #include "NonlinearSystem.h"
14 : #include "MooseVariable.h"
15 : #include "Function.h"
16 :
17 : #include "libmesh/quadrature.h"
18 : #include "libmesh/utility.h"
19 :
20 : registerMooseObject("TensorMechanicsApp", ComputeIncrementalBeamStrain);
21 :
22 : InputParameters
23 698 : ComputeIncrementalBeamStrain::validParams()
24 : {
25 698 : InputParameters params = Material::validParams();
26 698 : params.addClassDescription("Compute a infinitesimal/large strain increment for the beam.");
27 1396 : params.addRequiredCoupledVar(
28 : "rotations", "The rotations appropriate for the simulation geometry and coordinate system");
29 1396 : params.addRequiredCoupledVar(
30 : "displacements",
31 : "The displacements appropriate for the simulation geometry and coordinate system");
32 1396 : params.addRequiredParam<RealGradient>("y_orientation",
33 : "Orientation of the y direction along "
34 : "with Iyy is provided. This should be "
35 : "perpendicular to the axis of the beam.");
36 1396 : params.addRequiredCoupledVar(
37 : "area",
38 : "Cross-section area of the beam. Can be supplied as either a number or a variable name.");
39 1396 : params.addCoupledVar("Ay",
40 : 0.0,
41 : "First moment of area of the beam about y axis. Can be supplied "
42 : "as either a number or a variable name.");
43 1396 : params.addCoupledVar("Az",
44 : 0.0,
45 : "First moment of area of the beam about z axis. Can be supplied "
46 : "as either a number or a variable name.");
47 1396 : params.addCoupledVar("Ix",
48 : "Second moment of area of the beam about x axis. Can be "
49 : "supplied as either a number or a variable name. Defaults to Iy+Iz.");
50 1396 : params.addRequiredCoupledVar("Iy",
51 : "Second moment of area of the beam about y axis. Can be "
52 : "supplied as either a number or a variable name.");
53 1396 : params.addRequiredCoupledVar("Iz",
54 : "Second moment of area of the beam about z axis. Can be "
55 : "supplied as either a number or a variable name.");
56 1396 : params.addParam<bool>("large_strain", false, "Set to true if large strain are to be calculated.");
57 1396 : params.addParam<std::vector<MaterialPropertyName>>(
58 : "eigenstrain_names",
59 : {},
60 : "List of beam eigenstrains to be applied in this strain calculation.");
61 1396 : params.addParam<FunctionName>(
62 : "elasticity_prefactor",
63 : "Optional function to use as a scalar prefactor on the elasticity vector for the beam.");
64 698 : return params;
65 0 : }
66 :
67 524 : ComputeIncrementalBeamStrain::ComputeIncrementalBeamStrain(const InputParameters & parameters)
68 : : Material(parameters),
69 524 : _has_Ix(isParamValid("Ix")),
70 524 : _nrot(coupledComponents("rotations")),
71 524 : _ndisp(coupledComponents("displacements")),
72 524 : _rot_num(_nrot),
73 524 : _disp_num(_ndisp),
74 524 : _area(coupledValue("area")),
75 524 : _Ay(coupledValue("Ay")),
76 524 : _Az(coupledValue("Az")),
77 524 : _Iy(coupledValue("Iy")),
78 524 : _Iz(coupledValue("Iz")),
79 524 : _Ix(_has_Ix ? coupledValue("Ix") : _zero),
80 524 : _original_length(declareProperty<Real>("original_length")),
81 524 : _total_rotation(declareProperty<RankTwoTensor>("total_rotation")),
82 524 : _total_disp_strain(declareProperty<RealVectorValue>("total_disp_strain")),
83 524 : _total_rot_strain(declareProperty<RealVectorValue>("total_rot_strain")),
84 1048 : _total_disp_strain_old(getMaterialPropertyOld<RealVectorValue>("total_disp_strain")),
85 1048 : _total_rot_strain_old(getMaterialPropertyOld<RealVectorValue>("total_rot_strain")),
86 524 : _mech_disp_strain_increment(declareProperty<RealVectorValue>("mech_disp_strain_increment")),
87 524 : _mech_rot_strain_increment(declareProperty<RealVectorValue>("mech_rot_strain_increment")),
88 1048 : _material_stiffness(getMaterialPropertyByName<RealVectorValue>("material_stiffness")),
89 524 : _K11(declareProperty<RankTwoTensor>("Jacobian_11")),
90 524 : _K21_cross(declareProperty<RankTwoTensor>("Jacobian_12")),
91 524 : _K21(declareProperty<RankTwoTensor>("Jacobian_21")),
92 524 : _K22(declareProperty<RankTwoTensor>("Jacobian_22")),
93 524 : _K22_cross(declareProperty<RankTwoTensor>("Jacobian_22_cross")),
94 1048 : _large_strain(getParam<bool>("large_strain")),
95 1048 : _eigenstrain_names(getParam<std::vector<MaterialPropertyName>>("eigenstrain_names")),
96 524 : _disp_eigenstrain(_eigenstrain_names.size()),
97 524 : _rot_eigenstrain(_eigenstrain_names.size()),
98 524 : _disp_eigenstrain_old(_eigenstrain_names.size()),
99 524 : _rot_eigenstrain_old(_eigenstrain_names.size()),
100 524 : _nonlinear_sys(_fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0)),
101 524 : _soln_disp_index_0(_ndisp),
102 524 : _soln_disp_index_1(_ndisp),
103 524 : _soln_rot_index_0(_ndisp),
104 524 : _soln_rot_index_1(_ndisp),
105 524 : _initial_rotation(declareProperty<RankTwoTensor>("initial_rotation")),
106 524 : _effective_stiffness(declareProperty<Real>("effective_stiffness")),
107 1048 : _prefactor_function(isParamValid("elasticity_prefactor") ? &getFunction("elasticity_prefactor")
108 1048 : : nullptr)
109 : {
110 : // Checking for consistency between length of the provided displacements and rotations vector
111 524 : if (_ndisp != _nrot)
112 1 : mooseError("ComputeIncrementalBeamStrain: The number of variables supplied in 'displacements' "
113 : "and 'rotations' must match.");
114 :
115 : // fetch coupled variables and gradients (as stateful properties if necessary)
116 2092 : for (unsigned int i = 0; i < _ndisp; ++i)
117 : {
118 3138 : MooseVariable * disp_variable = getVar("displacements", i);
119 1569 : _disp_num[i] = disp_variable->number();
120 :
121 3138 : MooseVariable * rot_variable = getVar("rotations", i);
122 1569 : _rot_num[i] = rot_variable->number();
123 : }
124 :
125 523 : if (_large_strain && (_Ay[0] > 0.0 || _Ay[1] > 0.0 || _Az[0] > 0.0 || _Az[1] > 0.0))
126 1 : mooseError("ComputeIncrementalBeamStrain: Large strain calculation does not currently "
127 : "support asymmetric beam configurations with non-zero first or third moments of "
128 : "area.");
129 :
130 543 : for (unsigned int i = 0; i < _eigenstrain_names.size(); ++i)
131 : {
132 42 : _disp_eigenstrain[i] = &getMaterialProperty<RealVectorValue>("disp_" + _eigenstrain_names[i]);
133 42 : _rot_eigenstrain[i] = &getMaterialProperty<RealVectorValue>("rot_" + _eigenstrain_names[i]);
134 21 : _disp_eigenstrain_old[i] =
135 42 : &getMaterialPropertyOld<RealVectorValue>("disp_" + _eigenstrain_names[i]);
136 21 : _rot_eigenstrain_old[i] =
137 42 : &getMaterialPropertyOld<RealVectorValue>("rot_" + _eigenstrain_names[i]);
138 : }
139 522 : }
140 :
141 : void
142 3881 : ComputeIncrementalBeamStrain::initQpStatefulProperties()
143 : {
144 : // compute initial orientation of the beam for calculating initial rotation matrix
145 : const std::vector<RealGradient> * orientation =
146 3881 : &_subproblem.assembly(_tid, _nonlinear_sys.number()).getFE(FEType(), 1)->get_dxyzdxi();
147 3881 : RealGradient x_orientation = (*orientation)[0];
148 3881 : x_orientation /= x_orientation.norm();
149 :
150 7762 : RealGradient y_orientation = getParam<RealGradient>("y_orientation");
151 3881 : y_orientation /= y_orientation.norm();
152 3881 : Real sum = x_orientation(0) * y_orientation(0) + x_orientation(1) * y_orientation(1) +
153 3881 : x_orientation(2) * y_orientation(2);
154 :
155 3881 : if (std::abs(sum) > 1e-4)
156 1 : mooseError("ComputeIncrementalBeamStrain: y_orientation should be perpendicular to "
157 : "the axis of the beam.");
158 :
159 : // Calculate z orientation as a cross product of the x and y orientations
160 : RealGradient z_orientation;
161 3880 : z_orientation(0) = (x_orientation(1) * y_orientation(2) - x_orientation(2) * y_orientation(1));
162 3880 : z_orientation(1) = (x_orientation(2) * y_orientation(0) - x_orientation(0) * y_orientation(2));
163 3880 : z_orientation(2) = (x_orientation(0) * y_orientation(1) - x_orientation(1) * y_orientation(0));
164 :
165 : // Rotation matrix from global to original beam local configuration
166 3880 : _original_local_config(0, 0) = x_orientation(0);
167 3880 : _original_local_config(0, 1) = x_orientation(1);
168 3880 : _original_local_config(0, 2) = x_orientation(2);
169 3880 : _original_local_config(1, 0) = y_orientation(0);
170 3880 : _original_local_config(1, 1) = y_orientation(1);
171 3880 : _original_local_config(1, 2) = y_orientation(2);
172 3880 : _original_local_config(2, 0) = z_orientation(0);
173 3880 : _original_local_config(2, 1) = z_orientation(1);
174 3880 : _original_local_config(2, 2) = z_orientation(2);
175 :
176 3880 : _total_rotation[_qp] = _original_local_config;
177 :
178 : RealVectorValue temp;
179 3880 : _total_disp_strain[_qp] = temp;
180 3880 : _total_rot_strain[_qp] = temp;
181 3880 : }
182 :
183 : void
184 204144 : ComputeIncrementalBeamStrain::computeProperties()
185 : {
186 : // fetch the two end nodes for current element
187 : std::vector<const Node *> node;
188 612432 : for (unsigned int i = 0; i < 2; ++i)
189 408288 : node.push_back(_current_elem->node_ptr(i));
190 :
191 : // calculate original length of a beam element
192 : // Nodal positions do not change with time as undisplaced mesh is used by material classes by
193 : // default
194 : RealGradient dxyz;
195 816576 : for (unsigned int i = 0; i < _ndisp; ++i)
196 612432 : dxyz(i) = (*node[1])(i) - (*node[0])(i);
197 :
198 204144 : _original_length[0] = dxyz.norm();
199 :
200 : // Fetch the solution for the two end nodes at time t
201 204144 : const NumericVector<Number> & sol = *_nonlinear_sys.currentSolution();
202 204144 : const NumericVector<Number> & sol_old = _nonlinear_sys.solutionOld();
203 :
204 816576 : for (unsigned int i = 0; i < _ndisp; ++i)
205 : {
206 612432 : _soln_disp_index_0[i] = node[0]->dof_number(_nonlinear_sys.number(), _disp_num[i], 0);
207 612432 : _soln_disp_index_1[i] = node[1]->dof_number(_nonlinear_sys.number(), _disp_num[i], 0);
208 612432 : _soln_rot_index_0[i] = node[0]->dof_number(_nonlinear_sys.number(), _rot_num[i], 0);
209 612432 : _soln_rot_index_1[i] = node[1]->dof_number(_nonlinear_sys.number(), _rot_num[i], 0);
210 :
211 612432 : _disp0(i) = sol(_soln_disp_index_0[i]) - sol_old(_soln_disp_index_0[i]);
212 612432 : _disp1(i) = sol(_soln_disp_index_1[i]) - sol_old(_soln_disp_index_1[i]);
213 612432 : _rot0(i) = sol(_soln_rot_index_0[i]) - sol_old(_soln_rot_index_0[i]);
214 612432 : _rot1(i) = sol(_soln_rot_index_1[i]) - sol_old(_soln_rot_index_1[i]);
215 : }
216 :
217 : // For small rotation problems, the rotation matrix is essentially the transformation from the
218 : // global to original beam local configuration and is never updated. This method has to be
219 : // overriden for scenarios with finite rotation
220 204144 : computeRotation();
221 204144 : _initial_rotation[0] = _original_local_config;
222 :
223 612432 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
224 408288 : computeQpStrain();
225 :
226 204144 : if (_fe_problem.currentlyComputingJacobian())
227 48262 : computeStiffnessMatrix();
228 204144 : }
229 :
230 : void
231 408288 : ComputeIncrementalBeamStrain::computeQpStrain()
232 : {
233 408288 : const Real A_avg = (_area[0] + _area[1]) / 2.0;
234 408288 : const Real Iz_avg = (_Iz[0] + _Iz[1]) / 2.0;
235 408288 : Real Ix = _Ix[_qp];
236 408288 : if (!_has_Ix)
237 407968 : Ix = _Iy[_qp] + _Iz[_qp];
238 :
239 : // Rotate the gradient of displacements and rotations at t+delta t from global coordinate
240 : // frame to beam local coordinate frame
241 408288 : const RealVectorValue grad_disp_0(1.0 / _original_length[0] * (_disp1 - _disp0));
242 : const RealVectorValue grad_rot_0(1.0 / _original_length[0] * (_rot1 - _rot0));
243 : const RealVectorValue avg_rot(
244 408288 : 0.5 * (_rot0(0) + _rot1(0)), 0.5 * (_rot0(1) + _rot1(1)), 0.5 * (_rot0(2) + _rot1(2)));
245 :
246 408288 : _grad_disp_0_local_t = _total_rotation[0] * grad_disp_0;
247 408288 : _grad_rot_0_local_t = _total_rotation[0] * grad_rot_0;
248 408288 : _avg_rot_local_t = _total_rotation[0] * avg_rot;
249 :
250 : // displacement at any location on beam in local coordinate system at t
251 : // u_1 = u_n1 - rot_3 * y + rot_2 * z
252 : // u_2 = u_n2 - rot_1 * z
253 : // u_3 = u_n3 + rot_1 * y
254 : // where u_n1, u_n2, u_n3 are displacements at neutral axis
255 :
256 : // small strain
257 : // e_11 = u_1,1 = u_n1, 1 - rot_3, 1 * y + rot_2, 1 * z
258 : // e_12 = 2 * 0.5 * (u_1,2 + u_2,1) = (- rot_3 + u_n2,1 - rot_1,1 * z)
259 : // e_13 = 2 * 0.5 * (u_1,3 + u_3,1) = (rot_2 + u_n3,1 + rot_1,1 * y)
260 :
261 : // axial and shearing strains at each qp along the length of the beam
262 408288 : _mech_disp_strain_increment[_qp](0) = _grad_disp_0_local_t(0) * _area[_qp] -
263 408288 : _grad_rot_0_local_t(2) * _Ay[_qp] +
264 408288 : _grad_rot_0_local_t(1) * _Az[_qp];
265 408288 : _mech_disp_strain_increment[_qp](1) = -_avg_rot_local_t(2) * _area[_qp] +
266 408288 : _grad_disp_0_local_t(1) * _area[_qp] -
267 408288 : _grad_rot_0_local_t(0) * _Az[_qp];
268 408288 : _mech_disp_strain_increment[_qp](2) = _avg_rot_local_t(1) * _area[_qp] +
269 408288 : _grad_disp_0_local_t(2) * _area[_qp] +
270 408288 : _grad_rot_0_local_t(0) * _Ay[_qp];
271 :
272 : // rotational strains at each qp along the length of the beam
273 : // rot_strain_1 = integral(e_13 * y - e_12 * z) dA
274 : // rot_strain_2 = integral(e_11 * z) dA
275 : // rot_strain_3 = integral(e_11 * -y) dA
276 : // Iyz is the product moment of inertia which is zero for most cross-sections so it is assumed to
277 : // be zero for this analysis
278 : const Real Iyz = 0;
279 408288 : _mech_rot_strain_increment[_qp](0) =
280 408288 : _avg_rot_local_t(1) * _Ay[_qp] + _grad_disp_0_local_t(2) * _Ay[_qp] +
281 408288 : _grad_rot_0_local_t(0) * Ix + _avg_rot_local_t(2) * _Az[_qp] -
282 408288 : _grad_disp_0_local_t(1) * _Az[_qp];
283 408288 : _mech_rot_strain_increment[_qp](1) = _grad_disp_0_local_t(0) * _Az[_qp] -
284 408288 : _grad_rot_0_local_t(2) * Iyz +
285 408288 : _grad_rot_0_local_t(1) * _Iz[_qp];
286 408288 : _mech_rot_strain_increment[_qp](2) = -_grad_disp_0_local_t(0) * _Ay[_qp] +
287 408288 : _grad_rot_0_local_t(2) * _Iy[_qp] -
288 408288 : _grad_rot_0_local_t(1) * Iyz;
289 :
290 408288 : if (_large_strain)
291 : {
292 88960 : _mech_disp_strain_increment[_qp](0) +=
293 88960 : 0.5 *
294 88960 : ((Utility::pow<2>(_grad_disp_0_local_t(0)) + Utility::pow<2>(_grad_disp_0_local_t(1)) +
295 88960 : Utility::pow<2>(_grad_disp_0_local_t(2))) *
296 88960 : _area[_qp] +
297 88960 : Utility::pow<2>(_grad_rot_0_local_t(2)) * _Iy[_qp] +
298 88960 : Utility::pow<2>(_grad_rot_0_local_t(1)) * _Iz[_qp] +
299 88960 : Utility::pow<2>(_grad_rot_0_local_t(0)) * Ix);
300 88960 : _mech_disp_strain_increment[_qp](1) += (-_avg_rot_local_t(2) * _grad_disp_0_local_t(0) +
301 88960 : _avg_rot_local_t(0) * _grad_disp_0_local_t(2)) *
302 : _area[_qp];
303 88960 : _mech_disp_strain_increment[_qp](2) += (_avg_rot_local_t(1) * _grad_disp_0_local_t(0) -
304 88960 : _avg_rot_local_t(0) * _grad_disp_0_local_t(1)) *
305 : _area[_qp];
306 :
307 88960 : _mech_rot_strain_increment[_qp](0) += -_avg_rot_local_t(1) * _grad_rot_0_local_t(2) * _Iy[_qp] +
308 88960 : _avg_rot_local_t(2) * _grad_rot_0_local_t(1) * _Iz[_qp];
309 88960 : _mech_rot_strain_increment[_qp](1) += (_grad_disp_0_local_t(0) * _grad_rot_0_local_t(1) -
310 88960 : _grad_disp_0_local_t(1) * _grad_rot_0_local_t(0)) *
311 : _Iz[_qp];
312 88960 : _mech_rot_strain_increment[_qp](2) += -(_grad_disp_0_local_t(2) * _grad_rot_0_local_t(0) -
313 88960 : _grad_disp_0_local_t(0) * _grad_rot_0_local_t(2)) *
314 : _Iy[_qp];
315 : }
316 :
317 408288 : _total_disp_strain[_qp] = _total_rotation[0].transpose() * _mech_disp_strain_increment[_qp] +
318 408288 : _total_disp_strain_old[_qp];
319 408288 : _total_rot_strain[_qp] = _total_rotation[0].transpose() * _mech_rot_strain_increment[_qp] +
320 408288 : _total_disp_strain_old[_qp];
321 :
322 : // Convert eigenstrain increment from global to beam local coordinate system and remove eigen
323 : // strain increment
324 408768 : for (unsigned int i = 0; i < _eigenstrain_names.size(); ++i)
325 : {
326 480 : _mech_disp_strain_increment[_qp] -=
327 960 : _total_rotation[0] * ((*_disp_eigenstrain[i])[_qp] - (*_disp_eigenstrain_old[i])[_qp]) *
328 480 : _area[_qp];
329 480 : _mech_rot_strain_increment[_qp] -=
330 960 : _total_rotation[0] * ((*_rot_eigenstrain[i])[_qp] - (*_rot_eigenstrain_old[i])[_qp]);
331 : }
332 :
333 408288 : Real c1_paper = std::sqrt(_material_stiffness[0](0));
334 408288 : Real c2_paper = std::sqrt(_material_stiffness[0](1));
335 :
336 408288 : Real effec_stiff_1 = std::max(c1_paper, c2_paper);
337 :
338 408288 : Real effec_stiff_2 = 2 / (c2_paper * std::sqrt(A_avg / Iz_avg));
339 :
340 655416 : _effective_stiffness[_qp] = std::max(effec_stiff_1, _original_length[0] / effec_stiff_2);
341 :
342 408288 : if (_prefactor_function)
343 0 : _effective_stiffness[_qp] *= std::sqrt(_prefactor_function->value(_t, _q_point[_qp]));
344 408288 : }
345 :
346 : void
347 48262 : ComputeIncrementalBeamStrain::computeStiffnessMatrix()
348 : {
349 48262 : const Real youngs_modulus = _material_stiffness[0](0);
350 48262 : const Real shear_modulus = _material_stiffness[0](1);
351 :
352 48262 : const Real A_avg = (_area[0] + _area[1]) / 2.0;
353 48262 : const Real Iy_avg = (_Iy[0] + _Iy[1]) / 2.0;
354 48262 : const Real Iz_avg = (_Iz[0] + _Iz[1]) / 2.0;
355 48262 : Real Ix_avg = (_Ix[0] + _Ix[1]) / 2.0;
356 48262 : if (!_has_Ix)
357 48242 : Ix_avg = Iy_avg + Iz_avg;
358 :
359 : // K = |K11 K12|
360 : // |K21 K22|
361 :
362 : // relation between translational displacements at node 0 and translational forces at node 0
363 48262 : RankTwoTensor K11_local;
364 : K11_local.zero();
365 48262 : K11_local(0, 0) = youngs_modulus * A_avg / _original_length[0];
366 48262 : K11_local(1, 1) = shear_modulus * A_avg / _original_length[0];
367 48262 : K11_local(2, 2) = shear_modulus * A_avg / _original_length[0];
368 48262 : _K11[0] = _total_rotation[0].transpose() * K11_local * _total_rotation[0];
369 :
370 : // relation between displacements at node 0 and rotational moments at node 0
371 48262 : RankTwoTensor K21_local;
372 : K21_local.zero();
373 48262 : K21_local(2, 1) = shear_modulus * A_avg * 0.5;
374 48262 : K21_local(1, 2) = -shear_modulus * A_avg * 0.5;
375 48262 : _K21[0] = _total_rotation[0].transpose() * K21_local * _total_rotation[0];
376 :
377 : // relation between rotations at node 0 and rotational moments at node 0
378 48262 : RankTwoTensor K22_local;
379 : K22_local.zero();
380 48262 : K22_local(0, 0) = shear_modulus * Ix_avg / _original_length[0];
381 48262 : K22_local(1, 1) = youngs_modulus * Iz_avg / _original_length[0] +
382 48262 : shear_modulus * A_avg * _original_length[0] / 4.0;
383 48262 : K22_local(2, 2) = youngs_modulus * Iy_avg / _original_length[0] +
384 48262 : shear_modulus * A_avg * _original_length[0] / 4.0;
385 48262 : _K22[0] = _total_rotation[0].transpose() * K22_local * _total_rotation[0];
386 :
387 : // relation between rotations at node 0 and rotational moments at node 1
388 48262 : RankTwoTensor K22_local_cross = -K22_local;
389 48262 : K22_local_cross(1, 1) += 2.0 * shear_modulus * A_avg * _original_length[0] / 4.0;
390 48262 : K22_local_cross(2, 2) += 2.0 * shear_modulus * A_avg * _original_length[0] / 4.0;
391 48262 : _K22_cross[0] = _total_rotation[0].transpose() * K22_local_cross * _total_rotation[0];
392 :
393 : // relation between displacements at node 0 and rotational moments at node 1
394 48262 : _K21_cross[0] = -_K21[0];
395 :
396 : // stiffness matrix for large strain
397 48262 : if (_large_strain)
398 : {
399 : // k1_large is the stiffness matrix obtained from sigma_xx * d(epsilon_xx)
400 8948 : RankTwoTensor k1_large_11;
401 : // row 1
402 8948 : k1_large_11(0, 0) = Utility::pow<2>(_grad_disp_0_local_t(0)) +
403 8948 : 1.5 * Utility::pow<2>(_grad_rot_0_local_t(2)) * Iy_avg +
404 8948 : 1.5 * Utility::pow<2>(_grad_rot_0_local_t(1)) * Iz_avg +
405 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(1)) +
406 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(2)) +
407 8948 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(0)) * Ix_avg;
408 8948 : k1_large_11(1, 0) = 0.5 * _grad_disp_0_local_t(0) * _grad_disp_0_local_t(1) -
409 8948 : 1.0 / 3.0 * _grad_rot_0_local_t(0) * _grad_rot_0_local_t(1) * Iz_avg;
410 8948 : k1_large_11(2, 0) = 0.5 * _grad_disp_0_local_t(0) * _grad_disp_0_local_t(2) -
411 8948 : 1.0 / 3.0 * _grad_rot_0_local_t(0) * _grad_rot_0_local_t(2) * Iy_avg;
412 :
413 : // row 2
414 8948 : k1_large_11(0, 1) = k1_large_11(1, 0);
415 8948 : k1_large_11(1, 1) = Utility::pow<2>(_grad_disp_0_local_t(1)) +
416 8948 : 1.5 * Utility::pow<2>(_grad_rot_0_local_t(0)) * Iz_avg +
417 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(0)) +
418 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(2)) +
419 8948 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(2)) * Iy_avg +
420 8948 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(1)) * Iz_avg +
421 8948 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(0)) * Iy_avg;
422 8948 : k1_large_11(2, 1) = 0.5 * _grad_disp_0_local_t(1) * _grad_disp_0_local_t(2);
423 :
424 : // row 3
425 8948 : k1_large_11(0, 2) = k1_large_11(2, 0);
426 8948 : k1_large_11(1, 2) = k1_large_11(2, 1);
427 8948 : k1_large_11(2, 2) = Utility::pow<2>(_grad_disp_0_local_t(2)) +
428 8948 : 1.5 * Utility::pow<2>(_grad_rot_0_local_t(0)) * Iy_avg +
429 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(0)) +
430 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(1)) +
431 8948 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(0)) * Iz_avg +
432 8948 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(2)) * Iy_avg +
433 8948 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(2)) * Iz_avg;
434 :
435 8948 : k1_large_11 *= 1.0 / 4.0 / Utility::pow<2>(_original_length[0]);
436 :
437 8948 : RankTwoTensor k1_large_21;
438 : // row 1
439 8948 : k1_large_21(0, 0) = 0.5 * _grad_disp_0_local_t(0) * _grad_rot_0_local_t(0) * (Ix_avg)-1.0 /
440 8948 : 3.0 * _grad_disp_0_local_t(1) * _grad_rot_0_local_t(1) * Iz_avg -
441 8948 : 1.0 / 3.0 * _grad_disp_0_local_t(2) * _grad_rot_0_local_t(2);
442 8948 : k1_large_21(1, 0) = 1.5 * _grad_disp_0_local_t(0) * _grad_rot_0_local_t(1) * Iz_avg -
443 8948 : 1.0 / 3.0 * _grad_disp_0_local_t(1) * _grad_rot_0_local_t(0) * Iz_avg;
444 8948 : k1_large_21(2, 0) = 1.5 * _grad_disp_0_local_t(0) * _grad_rot_0_local_t(2) * Iy_avg -
445 8948 : 1.0 / 3.0 * _grad_disp_0_local_t(2) * _grad_rot_0_local_t(0) * Iy_avg;
446 :
447 : // row 2
448 8948 : k1_large_21(0, 1) = k1_large_21(1, 0);
449 8948 : k1_large_21(1, 1) = 0.5 * _grad_disp_0_local_t(1) * _grad_rot_0_local_t(1) * Iz_avg -
450 8948 : 1.0 / 3.0 * _grad_disp_0_local_t(0) * _grad_rot_0_local_t(0) * Iz_avg;
451 8948 : k1_large_21(2, 1) = 0.5 * _grad_disp_0_local_t(1) * _grad_rot_0_local_t(2) * Iy_avg;
452 :
453 : // row 3
454 8948 : k1_large_21(0, 2) = k1_large_21(2, 0);
455 8948 : k1_large_21(1, 2) = k1_large_21(2, 1);
456 8948 : k1_large_21(2, 2) = 0.5 * _grad_disp_0_local_t(2) * _grad_rot_0_local_t(2) * Iy_avg -
457 8948 : 1.0 / 3.0 * _grad_disp_0_local_t(0) * _grad_rot_0_local_t(0) * Iy_avg;
458 8948 : k1_large_21 *= 1.0 / 4.0 / Utility::pow<2>(_original_length[0]);
459 :
460 8948 : RankTwoTensor k1_large_22;
461 : // row 1
462 8948 : k1_large_22(0, 0) = Utility::pow<2>(_grad_rot_0_local_t(0)) * Utility::pow<2>(Ix_avg) +
463 8948 : 1.5 * Utility::pow<2>(_grad_disp_0_local_t(1)) * Iz_avg +
464 8948 : 1.5 * Utility::pow<2>(_grad_disp_0_local_t(2)) * Iy_avg +
465 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(0)) * Ix_avg +
466 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(2)) * Iz_avg +
467 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(1)) * Iy_avg +
468 8948 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(2)) * Iy_avg * Ix_avg +
469 8948 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(1)) * Iz_avg * Ix_avg;
470 8948 : k1_large_22(1, 0) = 0.5 * _grad_rot_0_local_t(0) * _grad_rot_0_local_t(1) * Iz_avg * Ix_avg -
471 8948 : 1.0 / 3.0 * _grad_disp_0_local_t(0) * _grad_disp_0_local_t(1) * Iz_avg;
472 8948 : k1_large_22(2, 0) = 0.5 * _grad_rot_0_local_t(0) * _grad_rot_0_local_t(2) * Iy_avg * Ix_avg -
473 8948 : 1.0 / 3.0 * _grad_disp_0_local_t(0) * _grad_disp_0_local_t(2) * Iy_avg;
474 :
475 : // row 2
476 8948 : k1_large_22(0, 1) = k1_large_22(1, 0);
477 8948 : k1_large_22(1, 1) = Utility::pow<2>(_grad_rot_0_local_t(1)) * Iz_avg * Iz_avg +
478 8948 : 1.5 * Utility::pow<2>(_grad_disp_0_local_t(0)) * Iz_avg +
479 8948 : 1.5 * Utility::pow<2>(_grad_rot_0_local_t(2)) * Iy_avg * Iz_avg +
480 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(1)) * Iz_avg +
481 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(2)) * Iz_avg +
482 8948 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(0)) * Iz_avg * Ix_avg;
483 8948 : k1_large_22(2, 1) = 1.5 * _grad_rot_0_local_t(1) * _grad_rot_0_local_t(2) * Iy_avg * Iz_avg;
484 :
485 : // row 3
486 8948 : k1_large_22(0, 2) = k1_large_22(2, 0);
487 8948 : k1_large_22(1, 2) = k1_large_22(2, 1);
488 8948 : k1_large_22(2, 2) = Utility::pow<2>(_grad_rot_0_local_t(2)) * Iy_avg * Iy_avg +
489 8948 : 1.5 * Utility::pow<2>(_grad_disp_0_local_t(0)) * Iy_avg +
490 8948 : 1.5 * Utility::pow<2>(_grad_rot_0_local_t(1)) * Iy_avg * Iz_avg +
491 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(1)) * Iy_avg +
492 8948 : 0.5 * Utility::pow<2>(_grad_disp_0_local_t(2)) * Iy_avg +
493 : 0.5 * Utility::pow<2>(_grad_rot_0_local_t(0)) * Iz_avg * Ix_avg;
494 :
495 8948 : k1_large_22 *= 1.0 / 4.0 / Utility::pow<2>(_original_length[0]);
496 :
497 : // k2_large and k3_large are contributions from tau_xy * d(gamma_xy) and tau_xz * d(gamma_xz)
498 : // k2_large for node 1 is negative of that for node 0
499 8948 : RankTwoTensor k2_large_11;
500 : // col 1
501 8948 : k2_large_11(0, 0) =
502 8948 : 0.25 * Utility::pow<2>(_avg_rot_local_t(2)) + 0.25 * Utility::pow<2>(_avg_rot_local_t(1));
503 8948 : k2_large_11(1, 0) = -1.0 / 6.0 * _avg_rot_local_t(0) * _avg_rot_local_t(1);
504 8948 : k2_large_11(2, 0) = -1.0 / 6.0 * _avg_rot_local_t(0) * _avg_rot_local_t(2);
505 :
506 : // col 2
507 8948 : k2_large_11(0, 1) = k2_large_11(1, 0);
508 8948 : k2_large_11(1, 1) = 0.25 * _avg_rot_local_t(0);
509 :
510 : // col 3
511 8948 : k2_large_11(0, 2) = k2_large_11(2, 0);
512 8948 : k2_large_11(2, 2) = 0.25 * Utility::pow<2>(_avg_rot_local_t(0));
513 :
514 8948 : k2_large_11 *= 1.0 / 4.0 / Utility::pow<2>(_original_length[0]);
515 :
516 8948 : RankTwoTensor k2_large_22;
517 : // col1
518 8948 : k2_large_22(0, 0) = 0.25 * Utility::pow<2>(_avg_rot_local_t(0)) * Ix_avg;
519 8948 : k2_large_22(1, 0) = 1.0 / 6.0 * _avg_rot_local_t(0) * _avg_rot_local_t(1) * Iz_avg;
520 8948 : k2_large_22(2, 0) = 1.0 / 6.0 * _avg_rot_local_t(0) * _avg_rot_local_t(2) * Iy_avg;
521 :
522 : // col2
523 8948 : k2_large_22(0, 1) = k2_large_22(1, 0);
524 8948 : k2_large_22(1, 1) = 0.25 * Utility::pow<2>(_avg_rot_local_t(2)) * Iz_avg +
525 8948 : 0.25 * Utility::pow<2>(_avg_rot_local_t(1)) * Iz_avg;
526 :
527 : // col3
528 8948 : k2_large_22(0, 2) = k2_large_22(2, 0);
529 8948 : k2_large_22(2, 2) = 0.25 * Utility::pow<2>(_avg_rot_local_t(2)) * Iy_avg +
530 8948 : 0.25 * Utility::pow<2>(_avg_rot_local_t(1)) * Iy_avg;
531 :
532 8948 : k2_large_22 *= 1.0 / 4.0 / Utility::pow<2>(_original_length[0]);
533 :
534 : // k3_large for node 1 is same as that for node 0
535 8948 : RankTwoTensor k3_large_22;
536 : // col1
537 8948 : k3_large_22(0, 0) = 0.25 * Utility::pow<2>(_grad_disp_0_local_t(2)) +
538 8948 : 0.25 * _grad_rot_0_local_t(0) * Ix_avg +
539 8948 : 0.25 * Utility::pow<2>(_grad_disp_0_local_t(1));
540 8948 : k3_large_22(1, 0) = -1.0 / 6.0 * _grad_disp_0_local_t(0) * _grad_disp_0_local_t(1) +
541 8948 : 1.0 / 6.0 * _grad_rot_0_local_t(0) * _grad_rot_0_local_t(1) * Iz_avg;
542 8948 : k3_large_22(2, 0) = -1.0 / 6.0 * _grad_disp_0_local_t(0) * _grad_disp_0_local_t(2) +
543 8948 : 1.0 / 6.0 * _grad_rot_0_local_t(0) * _grad_rot_0_local_t(2) * Iy_avg;
544 :
545 : // col2
546 8948 : k3_large_22(0, 1) = k3_large_22(1, 0);
547 8948 : k3_large_22(2, 2) = 0.25 * Utility::pow<2>(_grad_disp_0_local_t(0)) +
548 8948 : 0.25 * _grad_rot_0_local_t(2) * Iy_avg +
549 8948 : 0.25 * _grad_rot_0_local_t(1) * Iz_avg;
550 :
551 : // col3
552 8948 : k3_large_22(0, 2) = k3_large_22(2, 0);
553 : k3_large_22(2, 2) = 0.25 * Utility::pow<2>(_grad_disp_0_local_t(0)) +
554 : 0.25 * _grad_rot_0_local_t(2) * Iy_avg +
555 : 0.25 * _grad_rot_0_local_t(1) * Iz_avg;
556 :
557 8948 : k3_large_22 *= 1.0 / 16.0;
558 :
559 8948 : RankTwoTensor k3_large_21;
560 : // col1
561 8948 : k3_large_21(0, 0) = -1.0 / 6.0 *
562 8948 : (_grad_disp_0_local_t(2) * _avg_rot_local_t(2) +
563 8948 : _grad_disp_0_local_t(1) * _avg_rot_local_t(1));
564 8948 : k3_large_21(1, 0) = 0.25 * _grad_disp_0_local_t(0) * _avg_rot_local_t(1) -
565 8948 : 1.0 / 6.0 * _grad_disp_0_local_t(1) * _avg_rot_local_t(0);
566 8948 : k3_large_21(2, 0) = 0.25 * _grad_disp_0_local_t(0) * _avg_rot_local_t(2) -
567 8948 : 1.0 / 6.0 * _grad_disp_0_local_t(2) * _avg_rot_local_t(0);
568 :
569 : // col2
570 8948 : k3_large_21(0, 1) = 0.25 * _grad_disp_0_local_t(1) * _avg_rot_local_t(0) -
571 8948 : 1.0 / 6.0 * _grad_disp_0_local_t(0) * _avg_rot_local_t(1);
572 8948 : k3_large_21(1, 1) = -1.0 / 6.0 * _grad_disp_0_local_t(0) * _avg_rot_local_t(0);
573 :
574 : // col3
575 8948 : k3_large_21(0, 2) = 0.25 * _grad_disp_0_local_t(2) * _avg_rot_local_t(0) -
576 8948 : 1.0 / 6.0 * _grad_disp_0_local_t(0) * _avg_rot_local_t(2);
577 8948 : k3_large_21(2, 2) = -1.0 / 6.0 * _grad_disp_0_local_t(0) * _avg_rot_local_t(0);
578 :
579 8948 : k3_large_21 *= 1.0 / 8.0 / _original_length[0];
580 :
581 8948 : RankTwoTensor k4_large_22;
582 : // col 1
583 8948 : k4_large_22(0, 0) = 0.25 * _grad_rot_0_local_t(0) * _avg_rot_local_t(0) * Ix_avg +
584 8948 : 1.0 / 6.0 * _grad_rot_0_local_t(2) * _avg_rot_local_t(2) * Iy_avg +
585 8948 : 1.0 / 6.0 * _grad_rot_0_local_t(1) * _avg_rot_local_t(1) * Iz_avg;
586 8948 : k4_large_22(1, 0) = 1.0 / 6.0 * _grad_rot_0_local_t(1) * _avg_rot_local_t(0) * Iz_avg;
587 8948 : k4_large_22(2, 0) = 1.0 / 6.0 * _grad_rot_0_local_t(2) * _avg_rot_local_t(0) * Iy_avg;
588 :
589 : // col2
590 8948 : k4_large_22(0, 1) = 1.0 / 6.0 * _grad_rot_0_local_t(0) * _avg_rot_local_t(1) * Iz_avg;
591 8948 : k4_large_22(1, 1) = 0.25 * _grad_rot_0_local_t(1) * _avg_rot_local_t(1) * Iz_avg +
592 8948 : 1.0 / 6.0 * _grad_rot_0_local_t(0) * _avg_rot_local_t(0) * Iz_avg;
593 8948 : k4_large_22(2, 1) = 0.25 * _grad_rot_0_local_t(1) * _avg_rot_local_t(2) * Iz_avg;
594 :
595 : // col 3
596 8948 : k4_large_22(0, 2) = 1.0 / 6.0 * _grad_rot_0_local_t(0) * _avg_rot_local_t(2) * Iy_avg;
597 8948 : k4_large_22(1, 2) = 0.25 * _grad_rot_0_local_t(2) * _avg_rot_local_t(1) * Iy_avg;
598 8948 : k4_large_22(2, 2) = 0.25 * _grad_rot_0_local_t(2) * _avg_rot_local_t(2) * Iy_avg +
599 8948 : 1.0 / 6.0 * _grad_rot_0_local_t(0) * _avg_rot_local_t(0) * Iy_avg;
600 :
601 8948 : k3_large_22 += 1.0 / 8.0 / _original_length[0] * (k4_large_22 + k4_large_22.transpose());
602 :
603 : // Assembling final matrix
604 8948 : _K11[0] += _total_rotation[0].transpose() * (k1_large_11 + k2_large_11) * _total_rotation[0];
605 8948 : _K22[0] += _total_rotation[0].transpose() * (k1_large_22 + k2_large_22 + k3_large_22) *
606 8948 : _total_rotation[0];
607 8948 : _K21[0] += _total_rotation[0].transpose() * (k1_large_21 + k3_large_21) * _total_rotation[0];
608 8948 : _K21_cross[0] +=
609 8948 : _total_rotation[0].transpose() * (-k1_large_21 + k3_large_21) * _total_rotation[0];
610 8948 : _K22_cross[0] += _total_rotation[0].transpose() * (-k1_large_22 - k2_large_22 + k3_large_22) *
611 17896 : _total_rotation[0];
612 : }
613 48262 : }
614 :
615 : void
616 159664 : ComputeIncrementalBeamStrain::computeRotation()
617 : {
618 159664 : _total_rotation[0] = _original_local_config;
619 159664 : }
|
