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