15 Solids
15.1 Kinematics
In the flex representation method, the kinematic description is written in terms of the reference and current trimmed U-spline manifolds. Consider a reference U-spline manifold with a right-handed orthogonal coordinate system and a motion or deformation of the reference U-spline manifold defined by the mapping where the mapped manifold is called the current or deformed manifold. We can measure the change in position of each coordinate under the mapping through a displacement field where The motion can also be written in terms of the displacement field as The deformation gradient is a linear operator defined as where and is the identity mapping on .
15.2 Local Mesh Size
Terms in many practical discrete formulations depend on the local element size. Examples include penalty methods, stabilized formulations, and shock-capturing viscosities. There is some ambiguity about what this may mean in a trimmed spline, but it is usually sufficient to consider the size of elements in the untrimmed spline. Even without considering trimming there is some ambiguity in what is meant by "element size" when a spline element may be distorted anisotropically and/or non-uniformly. A widely-used (e.g., [11]) approach to expressing element size in this setting is through the following tensor: where the untrimmed element parent coordinates are assumed to be in the parametric domain , in which is the parametric dimension of the element. The factor of 4 is to maintain consistency with many references from the stabilized fluid dynamics literature where the parent domain is instead assumed to be . For a uniform, isotropic, axis-aligned mesh with an element size of in each direction, we have Following the design of element size estimates for stabilized methods in [11], we can estimate the element size squared given only the tensor : This could be made exact for the isotropic mesh with a factor of , but this is typically omitted from isotropic element size estimates in the literature, as only the right asymptotics are typically needed, and this factor is effectively absorbed into other scaling factors in the stabilization terms under consideration.
A simpler estimate might initially appear to be , which would also give element size squared to within a constant factor in the isotropic case, or, even more simply, , which would estimate element size directly. However, a desired property of an element size estimate is that it will coverge to the smallest directional element size under anisotropic refinement. For example, if we have directional element sizes and in a two-dimensional mesh, equation would estimate which is a well-known regularization of . Thus, if a two-dimensional mesh were refined in only the direction while holding fixed, we would estimate as . This property would not hold for a determinant-based estimate of element size. The ability to recover the same asymptotic element size estimate from anisotropic refinement of a higher-dimensional mesh as from a true lower-dimensional mesh is another reason for omitting the factor of from the estimate equation .
For fluid dynamics applications, it is also common to estimate the element size along the direction of a vector . This is enabled by through the formula which is commonly used in calculating optimal stabilizing streamline diffusion or wall-normal element sizes in boundary layer models.
15.3 Weak Form
In the flex representation method, the physical domain of the problem is a trimmed U-spline manifold with boundary denoted by . The weak form of the equations of motion in the current configuration can be written as: Find such that, for all , where is the mass density of the material, is the time-dependent body force, is the time-dependent traction force, and is the Cauchy stress, which depends on the displacement and material constitutive model.
15.4 Semi-Discrete Form
The semi-discrete equations of motion are a system of coupled second order ordinary differential equations given as follows where is the consistent mass matrix, and denote the internal and the external forces, respectively, and damping has been ignored. The displacement, velocity, and acceleration solution vectors , and contain the solution coefficients with respect to the U-spline basis.
Equation can also be written in residual form
For the consistent mass matrix, the submatrix components are For each force vector, the subvector components are where is the Cauchy stress in Voigt form and are the standard strain-displacement matrices in the current configuration defined as
Unless otherwise indicated, to solve this discrete system of equations we use Newton-Raphson iteration to handle the nonlinearities within classical time-stepping schemes such as generalized-.
15.5 Penalty Enforcement of Dirichlet Boundary Conditions
This section discusses how a penalty regularization is used to approximate Dirichlet boundary conditions.
15.5.1 Displacement Dirichlet Conditions
To enforce Dirichlet boundary conditions on displacement, we augment the weak form of the problem with additional terms to enforce boundary conditions through a penalty regularization. In particular, we add to the left-hand side of equation , where is a penalty parameter, is the prescribed boundary displacement, is the test function of the weak formulation, and is the set whose displacement is constrained. The constrained set is often a surface, as suggested by notation, but the penalty formulation in this section can be applied to volumes, curves, or points as well.
15.5.2 Penalty Parameter Estimation
If the penalty parameter is too small, then it will allow excessive error in the boundary condition, whereas, if it is too large, it can cause a number of problems, including ill-conditioned algebraic equations, small critical time steps in explicit dynamics, and diminished solution quality in immersed discretzations, especially when the immersed boundary is skew to the background mesh. It is therefore important to be able to estimate the penalty parameter automatically and robustly. For displacement boundary conditions in solid mechanics, we use the formula where is a dimensionless scaling factor, is an estimate of the material’s Young’s modulus, as discussed further in Scalar Stiffness Estimation, is an estimate of the local element size, computed using equation , is the dimension of the body (e.g., three for a volume) and is the dimension of the set constrained by the boundary condition (e.g., two for a surface). For the dimensionless scaling factor, a value of is a reasonable default for practical engineering calculations across a wide range of problems, but higher values may be necessary for applications like patch tests, material model calibration, or fundamental convergence studies.
15.5.3 Interpretation of the Penalty Parameter
To interpret the meaning of the dimensionless penalty parameter in practical terms, it is helpful to consider the case of a one-dimensional linear-elastic bar in axial tension, with one of its ends constrained (implying ). The traction applied by the penalty term needs to balance the stress in the bar, and we can therefore relate the error in the boundary condition to the change in length of a single element. Force balance implies that the tensile stress in the bar must be , while the change in length of a single element is given by . Taken together with equation , this implies that the boundary condition error is times smaller than the change in length of an element, which converges to zero as and the overall specimen’s change in length is distributed over more elements.
It may at first be alarming that for (e.g., constraints on points or lines in a volumetric body) the penalty does not increase as . However, these are cases where we expect a zero reaction force in the continuous limit, because the solution to any nonzero concentrated load would be singular. Such constraints should only be used for situations where a zero (or in practice small) reaction force is expected, e.g., when removing rigid modes.
15.5.4 Velocity and Acceleration Dirichlet Conditions
In some applications, it may be more natural to apply Dirichlet boundary conditions to velocity or acceleration. Rather than formulating penalty terms involving velocity and acceleration, we convert velocity and acceleration boundary data into displacement data through numerical integration, then apply the displacement penalty formulation of equation . The reason for this is to be able to reuse the robust penalty estimation formula of equation and have a consistent interpretation of the boundary condition regularization. An implicit Crank–Nicolson scheme is used to accumulate the current prescribed displacement, regardless of what time integration method is used to solve the semi-discrete problem. If velocity data is given as a function of time at each boundary material point, then the displacement data at each point is accumulated like and if acceleration data is given, then velocity and displacement are accumulated like where and are indices of consecutive time points and is the time step size. The time points in this Crank–Nicolson scheme for boundary data correspond to times at which the boundary condition residual is assembled. These times are not necessarily the same time points considered to be integer-indexed time levels in the overall problem’s integrator. E.g., if an implicit midpoint scheme is applied to the overall problem, then the time points and in the above formulas for boundary data integration will correspond to midpoints of consecutive time steps used for the semi-discrete problem.
15.6 Time Stepping Algorithms
The implicit and explicit time stepping methods used in Coreform Flex are all variations on the generalized- method of [19] for the temporal discretization of the initial value problem given by equation . The base method is stated as follows: given find such that where is the time step and the real-valued parameters , , and define the method. The above statement of the method is generally implicit, but may produce explicit formulas for the -level unknowns for certain choices of parameters, as will be demonstrated in the following sections.
15.6.1 Initial Acceleration
One awkward implementational aspect of the generalized- method is that it requires an initial condition for acceleration to compute the first step. This is not part of the continuous problem statement, which only includes initial conditions for displacement and velocity, with the initial acceleration being uniquely determined by the physics of the problem. For explicit versions of the method, we compute this initial acceleration as where uses a lumped diagonal approximation for efficiency. For implicit versions of the formulation, we use an implementation from PETSc, which estimates initial acceleration by first performing two half-steps of size using the implicit first-order method given by generalized- parameters and , which has no dependence on . The two half-steps produce approximations and of the velocity solution at times and , in addition to the initial velocity given at . Using these approximations, the initial acceleration is then approximated using the following one-sided finite difference: where is subsequently recomputed using the full step size and original generalized- parameters, now that an initial acceleration is available.
15.6.2 Implicit Dynamics
In problems not requiring high-fidelity approximations of stress wave propagation, it is most efficient to use implicit members of the generalized- family, where has some nontrivial dependence on the -level solution, requiring the solution of some algebraic system of equations.
15.6.2.1 Implicit Dynamics With High-Frequency Damping
The most widely-used family of generalized- methods used for imlicit dynamics is parameterized by the single parameter . It was shown in [19] that by choosing the parameters we can achieve second order accuracy and unconditional stability for linear problems. The parameter has the interpretation of controlling numerical damping for high-frequency modes. More specifically, a linear elastodynamic problem can be expressed as a superposition of vibrational modes at different frequencies. Applying the generalized- scheme to the single degree-of-freedom ordinary differential equation governing one of these modes with angular frequency , the parameter is the spectral radius of the method’s amplification matrix (mapping the solution from one time step to the next), in the limit of . In other words, expresses how rapdily under-resolved vibrational modes decay. Choosing implies no dissipation and is equivalent to the implicit midpoint rule, while choosing produces maximum high-frequency dissipation.
The reason for introducing numerical dissipation in high-frequency modes is that, in semi-discrete problems that have been discretized in space, the highest-frequency modes are essentially artifacts of the spatial discretization. In nonlinear problems, these spurious modes can have nontrivial interactions with the lower-frequency modes, so it is desirable to eliminate them from the analysis. Note that for all values of , the limit of has no dissipation, and integrators using all choices of converge to the same undamped solution as . This is what distinguishes numerical damping from physical damping (e.g., mass- or stiffness-proproportional damping) that modifies the governing equations. The use of numerical dissipation therefore involves a tacit assumption that the time step is not significantly smaller than the periods of vibrational modes that are well-resolved by the spatial discretization. However, that is normally the case by a wide margin in problems where implicit analysis is cost-effective.
15.6.2.2 Backward Euler
The backward Euler method for the full dynamic problem cannot be derived from the generalized- scheme. However, we primarily use the backward Euler method for quasi-static problems with no dependence on the velocity or acceleration. This special case can still be formally derived from the generalized- formulas by taking and only using the first relation of , where the other parameters are not relevant. However, in practice, we use a distinct code path for the backward Euler integration of quasi-static problems.
15.6.2.3 Adaptive Time Stepping
In complicated nonlinear problems, the algebraic equations in each step of an implicit method may be very difficult to solve. They typically become easier to solve as the time step is reduced, for two reasons: first, because the converged solution of the previous step is closer to the converged solution of the next step and serves as a good starting point for nonlinear iteration and, second, because the tangent matrix approaches a solution-independent mass matrix which effectively reduces the system to -projection as . In practice, the ability to solve algebraic equations is often what limits time step size rather than accuracy considerations, and a suitable step size is difficult to estimate beforehand. As such, we include adaptive time stepping to re-run steps with reduced step sizes if the nonlinear solver is unable to converge. Our time step adaptivity is based solely on whether the nonlinear solve is able to converge, rather than using a local error estimator to control accuracy as is sometimes seen in the academic literature. The only control on accuracy is the maximum time step size, , which is specified in the input.
The adaptivity algorithm can be summarized as follows. If the nonlinear solve for a step fails to converge, the step is repeated with for until either it succeeds or the step size dips below a minimum threshold . Each time a step is successfully solved, the time step is updated by for .
15.6.3 Explicit Dynamics
Several different explicit schemes can be formally derived from the generalized- framework.
15.6.3.1 Explicit Central Difference
A type of explicit central difference scheme referred to as the Verlet method can be recovered from the generalized- method by setting and simplifying the solution of the implicit relation into the explicit formula where is replaced by a lumped diagonal approximation for efficiency. However, a modification of this known as the leapfrog method, where velocity is computed at half-integer steps, is easier to adapt to variable time step sizes. Specifically, and the same explicit formula equation for , where can vary from step to step.
Given an estimate for the maximum vibrational mode frequency of the semi-discrete problem, the computation of which is discussed further in Maximum Frequency Estimation, the critical time step used for the leapfrog method is given by which serves as an upper bound on . Rather than saturating on this upper bound, the time step is usually selected slightly lower than the critical value, to account for approximations in the computation of . The reduction of the time step relative to the theoretical critical time step is specified using a safety factor such that
15.6.3.2 Explicit Dynamics With High-Frequency Damping
The central difference method of Explicit Central Difference contains no numerical dissipation, and therefore carries some risk of polluting the computed structural response with effects of spurious high-frequency modes. Following [47], a family of explicit schemes with tunable high-frequency dissipation can be derived from the implicit generalized- formulas by selecting the parameters where is the spectral radius of the amplification matrix for the highest-frequency mode at a time step size where it is maximally damped. The explicit update formula for the -level acceleration is given by where is a diagonal lumped-mass approximation, while the displacement and velocity are updated following the usual generalized- formulas.
Time step estimation is more involved than in the central difference scheme, because it depends on the choice of . The maximum time step for which the method is stable is given by which is in general slightly greater than . This critical time step only coincides with when , where it is also equal to the critical time step for the central difference method. The time step at which high-frequency damping is controlled becomes moderately smaller when some numerical damping is added. The spectral radius of the amplification matrix for the highest-frequency mode goes from to 1 over the interval between and , so saturating on the critical time step would defeat the purpose of using a method with high-frequency dissipation. As in Explicit Central Difference, the critical time step is reduced by a user-specified saftey factor . The step size used for time integration is then computed as
15.6.3.3 Maximum Frequency Estimation
Our critical time step formulas for explicit integrators depend on being able to estimate the highest frequency of any vibrational mode in the semi-discrete problem. For linear elastodynamic problems, we can define a stiffness matrix such that Then the vibrational mode with angular frequency is the solution to the generalized eigenvalue problem The largest eigenvalue this problem can be found efficiently via power iteration when is symmetric, as is typical in linear elasticity. Power iteration consists of repeatedly computing starting from a randomly-oriented unit vector and continuing until reaching a fixed point where for some dimensionless tolerance . As the power iteration approaches a fixed point, the quantity takes on the interpretation of a Rayleigh quotient that can be used to compute the maximum vibrational frequency: To avoid forming the stiffness matrix in explicit calculations, we approximate its action using a finite difference: where is the current displacement. This will produce a tangent stiffness in the current configuration when applied to nonlinear problems. Although nonlinear problems are technically outside the scope of the modal stability analysis for time integrators, the resulting estimate is nonetheless effective in practice when combined with a reasonable safety factor for the critical time step.
15.6.3.4 Mass-Proportional Damping
The amount of numerical damping that can be incorporated into explicit calculations while remaining consistent with the original governing equations is limited by the upper bound on time step size imposed by stability considerations. If more dissipation is needed, one typically includes it as a regularization of the problem itself, by modifying the governing equations to include an artificial dissipative term. The simplest choice of dissipation is mass-proportional damping. At the semi-discrete level, this can be formulated as adding a term to in equation , where is the mass damping coefficient, with units of frequency. This internal force contribution is calculated using in the method of Explicit Dynamics With High-Frequency Damping, which is consistent with the choice of . In the leapfrog method of Explicit Central Difference, we only have velocity directly available at half steps, so the damping contribution is calculated using ; this is inconsistent with , but the method would become implicit if we attempted to use .
Mass damping can be used to define a damping ratio for a vibrational mode of frequency , A nonzero damping ratio impacts critical time step calculations. The critical time step for the method of Explicit Dynamics With High-Frequency Damping becomes where are intermediate variables introduced to simplify the formula, and the damping ratio is computed for the frequency . Note that the critical time step with nonzero damping ratio goes to zero if there is no numerical dissipation (), so it is recommended to use in combination with mass damping to avoid excessive reduction of the critical time step.
In the leapfrog central difference scheme of Explicit Central Difference, we do not adjust critical time step for mass damping, relying instead on the safety factor to absorb effects of artificial damping on stability.
15.7 Near-incompressibility
We employ a stabilized formulation to allow for locking-free equal-order U-spline interpolation of both displacement and pressure, avoiding the discrete saddle point problem that would necessitate careful choice of inf-sup stable interpolations in the incompressible limit. Although degree elevation of a U-spline basis resolves “locking” in the sense of artifically-stiff response, stresses are still very innaccurate for most problems.
Traditional low-order FEA (i.e., linear elements) often makes use of specialized quadrature schemes for the nearly-incompressible limit that are simple to implement and emulate the use of a discontinuous, elementwise-constant pressure space. In this case, the pressure-dependence of the problem can be condensed out at the element level, leaving a pure displacement-based formulation. For high-order, smooth U-splines, no such element-level condensation can be performed, and the design of reduced quadrature schemes that alleviate locking for trimmed splines remains an open research problem, so formulations must be developed that handle the pressure dependence explicitly.
15.7.1 Stabilized Static Pressure Equation
For static and implicit problems, we use a simplified version of the the pressure-stabilized formulation from [39], which can be stated as: Find such that, for all , where , is an estimate of the bulk modulus, as discussed further in Scalar Stiffness Estimation, and is an expression for the pressure, given entirely in terms of (derivatives of) the displacement field, in which defines the Cauchy stress for a given constitutive model entirely through kinematics, without invoking the unknown pressure field . The symbol , without a subscript “”, as it appears in the weak form, refers to an alternate expression for Cauchy stress, which defines the deviatoric part in terms of displacement, but uses the unknown field to define the hydrostatic part: where the deviatoric part of stress is defined in terms of displacement only as
The stabilization parameter is given by where is an estimate of the material’s shear modulus, as discussed further in Scalar Stiffness Estimation, is the degree of the basis, is a local proxy for element diameter whose square is estimated by equation , and is a dimensionless constant (with being a good default value).
The pressure stabilization perturbs the problem slightly, but the scaling of by maintains second-order accuracy of the displacement field. Following [39], it is technically possible to formulate higher-than-second-order-accurate stabilization by replacing with the full residual of the momentum-balance equation, but this imposes stability limits on and the time step that may be difficult to satisfy in practice and it requires costly computation of second-order spatial derivatives of the displacement field. The consistency error of the simplified formulation is usually much smaller than other sources of error in practical calculations, and higher-than-second-order convergence is often still seen empirically with high-order splines.
15.7.2 Stabilized Pressure Rate Equation
The absence of a time derivative in the static pressure equation makes it difficult to solve using lumped-mass explicit time integration. Thus, in explicit dynamic problems, we follow [87], modifying the pressure equation to be in rate form, viz., The pressure rate equation is also scaled by the estimated bulk modulus to obtain a more balanced system of equations. The rate-form pressure-stabilized problem can be stated as: Find such that, for all , where . For the rate-form pressure equation, the stabilization parameter is changed to Unlike in the static formulation, we retain inertial and body-force contributions to the momentum balance residual in the stabilization term, as these may be a significant source of pressure gradients in highly-dynamic problems and the scaling of stabilization by instead of makes the explicit formulation more sensitive to consistency error introduced by stabilization. This formulation could in principle be rendered strongly-consistent by including the divergence of shear stress alongside the pressure gradient, but this is left out for efficiency and robustness. (The cited reference [87] focuses exclusively on linear simplicial elements, and also omits the shear stress divergence.)
In many constitutive models (e.g., metals), can be expressed as a function of only the deformation gradient determinant, in which case where is the so-called “tangent bulk modulus” of the material.
This formulation softens the bulk response to prevent locking, but this softening can sometimes result in element inversion on coarse discretizations of problems with extreme deformations. In large-deformation explicit analyses, we include an additional nonlinear stabilization term in the formulation: where activates this term smoothly when is less than the prescribed threshold , which defaults to 1/2. Because in the exact solution, this term remains formally high-order-consistent, but it may stiffen the structural response on coarse meshes. (Note that taking would completely cancel the effect of the pressure stabilization, reverting to a pure-displacement formulation that locks in the nearly-incompressible limit.)
15.7.2.1 Explicit Time Integration
The leapfrog central difference scheme used for explicit displacement-based solid dynamics in Coreform IGA is very efficient for second-order-in-time problems, but would be unstable for the first-order pressure rate equation. This can be seen by noticing that, for fixed , the pressure is effectively governed by the heat equation with a source term from and a conductivity modulated by . Leapfrog is well-known to be unconditionally-unstable for the heat equation, and first-order systems generally.
There are two options for explicit time integration of pressure-stabilized problems. We recommend using a specialized explicit integrator for this system. It is a predictor-corrector based on the implicit generalized- method with that updates displacement first, followed by pressure in the correction step. Stabilized problems may also be integrated with a more general explicit predictor-multicorrector scheme using a monolithic displacement/pressure update and arbitrary generalized- parameters, but this typically requires multiple correction passes for stability and relies on implicit time step adaptivity rather than computing a critical time step.