김경종
23 KiB
q _ {i} = e _ {i j k} ^ {\varphi} \varepsilon_ {j k} + D _ {i j} ^ {\varphi} E _ {j}.
Equation 2.10.1-4
A choice of either stress coefficients or strain coefficients for the piezoelectric material matrix is available in the input definitions for ABAQUS.
Kinematics
For the piezoelectric elements both displacements and electric potentials exist at the nodal locations. The displacements and electrical potentials are approximated within the element as
\mathbf {u} = \mathbf {N} ^ {N} u ^ {N}
and
\varphi = \mathbf {N} ^ {N} \varphi^ {N},
where { \bf N } ^ { N } is the array of interpolating functions and u ^ { N } and \varphi ^ { N } are nodal quantities. The body forces and charges as well as the surface forces and charges are interpolated in a similar manner.
The strains and electrical potential gradients are given as
\pmb {\varepsilon} = \mathbf {B} _ {u} ^ {N} \boldsymbol {u} ^ {N}
and
\mathbf {E} = - \mathbf {B} _ {\varphi} ^ {N} \varphi^ {N},
where \mathbf { B } _ { u } ^ { N } and \mathbf { B } _ { \varphi } ^ { N } are the spatial derivatives of { { \bf N } ^ { N } } .
System equations
With these approximate fields and the constitutive properties given above, in conjunction with the equilibrium and conservation equations, the following system of equations is derived in terms of nodal quantities:
M ^ {M N} \ddot {u} ^ {N} + K _ {u u} ^ {M N} u ^ {N} + K _ {\varphi u} ^ {M N} \varphi^ {N} = P ^ {M}
Equation 2.10.1-5
and
K _ {\varphi u} ^ {M N} u ^ {N} - K _ {\varphi \varphi} ^ {M N} \varphi^ {N} = - Q ^ {M},
Equation 2.10.1-6
where
Procedures
M ^ {M N} = \int_ {V} \rho \mathbf {N} ^ {M} \cdot \mathbf {N} ^ {N} d V
is the mass matrix (no inertia terms exist for the electrical flux conservation equation), ½ is the mass density,
K _ {u u} ^ {M N} = \int_ {V} \mathbf {B} _ {u} ^ {M}: \mathbf {D} _ {m}: \mathbf {B} _ {u} ^ {N} d V
is the displacement stiffness matrix,
K _ {\varphi \varphi} ^ {M N} = \int_ {V} \mathbf {B} _ {\varphi} ^ {M} \cdot \mathbf {D} _ {\varphi} \cdot \mathbf {B} _ {\varphi} ^ {N} d V
is the dielectric "stiffness" matrix,
K _ {\varphi u} ^ {M N} = \int_ {V} \mathbf {B} _ {\varphi} ^ {M} \cdot \mathbf {e}: \mathbf {B} _ {u} ^ {N} d V
is the piezoelectric coupling matrix,
P ^ {M} = \int_ {V} \mathbf {N} ^ {M} \cdot \mathbf {P} _ {v} d V + \int_ {S} \mathbf {N} ^ {M} \cdot \mathbf {P} _ {s} d S + P _ {c} ^ {M}
is the mechanical force vector, and
Q ^ {M} = \int_ {V} \mathbf {N} ^ {M} \cdot \mathbf {Q} _ {v} d V + \int_ {S} \mathbf {N} ^ {M} \cdot \mathbf {Q} _ {s} d S + Q _ {c} ^ {M}
is the electrical charge vector. In these expressions the constitutive properties are specified in a matrix form where \mathbf { D } _ { m } is the mechanical relationship, \mathbf { D } _ { \varphi } is the electrical relationship, and e is the piezoelectrical relationship. The "load" vectors include the body, surface, and concentrated quantities, as shown. The unknowns are the nodal displacements and potentials. Once these are determined, the strains and potential gradients can be computed using the expressions given above. The stresses and electrical flux densities are computed by means of constitutive relationships.
2.11 Heat transfer
2.11.1 Uncoupled heat transfer analysis
The ABAQUS/Standard capability for uncoupled heat transfer analysis is intended to model solid body heat conduction with general, temperature-dependent conductivity; internal energy (including latent heat effects); and quite general convection and radiation boundary conditions. This section describes the basic energy balance, constitutive models, boundary conditions, finite element discretization, and time integration procedures used.
Heat transfer in flowing materials (convection) is discussed in ``Convection/diffusion,'' Section 2.11.3. Radiation heat transfer in cavities is discussed in ``Cavity radiation,'' Section 2.11.4. All such heat transfer mechanisms can be present in a model.
Energy balance
The basic energy balance is (Green and Naghdi)
Equation 2.11.1-1
\int_ {V} \rho \dot {U} d V = \int_ {S} q d S + \int_ {V} r d V,
where V is a volume of solid material, with surface area S ; \rho is the density of the material; \dot { U } is the material time rate of the internal energy; q is the heat flux per unit area of the body, flowing into the body; and r is the heat supplied externally into the body per unit volume.
It is assumed that the thermal and mechanical problems are uncoupled in the sense that U = U ( \theta ) only, where µ is the temperature of the material, and q and r do not depend on the strains or displacements of the body. For simplicity a Lagrangian description is assumed, so "volume" and "surface" mean the volume and surface in the reference configuration.
Constitutive definition
This relationship is usually written in terms of a specific heat, neglecting coupling between mechanical and thermal problems:
c (\theta) = \frac {d U}{d \theta},
except for latent heat effects at phase changes, which are given separately in terms of solidus and liquidus temperatures (the lower and upper temperature bounds of the phase change range) and the total internal energy associated with the phase change, called the latent heat. When latent heat is given, it is assumed to be in addition to the specific heat effect (see Figure 2.11.1-1). For many cases it is reasonable to assume that the phase change occurs within a known temperature range, and then the *LATENT HEAT option is used. However, in some cases it may be necessary to include a kinetic theory for the phase change to model the effect accurately (an example would be the prediction of crystallization in a polymer casting process). For such cases the user can model the process in considerable detail, using the solution-dependent state variable feature in ABAQUS, with user subroutine HETVAL.
Figure 2.11.1-1 Specific heat, latent heat definition.
Procedures
line
| θ, temperature | U, internal energy | c(θ) = dU/dθ |
|---|---|---|
| 0 | Low | 0 |
| Low | Medium | 0 |
| High | High | 0 |
| High | Very High | 0 |
Heat conduction is assumed to be governed by the Fourier law,
Equation 2.11.1-2
\mathbf {f} = - \mathbf {k} \frac {\partial \theta}{\partial \mathbf {x}},
where k is the conductivity matrix, \mathbf { k } = \mathbf { k } ( \theta ) ; f is the heat flux; and x is position. The conductivity k can be fully anisotropic, orthotropic, or isotropic.
Boundary conditions
Boundary conditions can be specified as prescribed temperature, \theta = \theta ( \mathbf { x } , t ) ; prescribed surface heat flux, q = q ( \mathbf { x } , t ) per area; prescribed volumetric heat flux, q = r ( \mathbf { x } , t ) per volume; surface convection: q = h ( \theta - \theta ^ { 0 } ) , where h = h ( \mathbf { x } , t ) is the film coefficient and \theta ^ { 0 } = \theta ^ { 0 } ( \mathbf { x } , t ) is the sink temperature; and
radiation: q = A \bigg ( ( \theta - \theta ^ { Z } ) ^ { 4 } - ( \theta ^ { 0 } - \theta ^ { Z } ) ^ { 4 } \bigg ) , where A is the radiation constant (emissivity times the Stefan-Boltzmann constant) and \theta ^ { Z } is the absolute zero on the temperature scale used. Surfaces can also participate in cavity radiation effects. The cavity radiation formulation in ABAQUS is described in ``Cavity radiation,'' Section 2.11.4.
Spatial discretization
A variational statement of the energy balance, Equation 2.11.1-1, together with the Fourier law, Equation 2.11.1-2, is obtained directly by the standard Galerkin approach as
Equation 2.11.1-3
\int_ {V} \rho \dot {U} \delta \theta d V + \int_ {V} \frac {\partial \delta \theta}{\partial \mathbf {x}} \cdot \cdot \mathbf {k} \frac {\partial \theta}{\partial \mathbf {x}} d V = \int_ {V} \delta \theta r d V + \int_ {s _ {q}} \delta \theta q d S,
where \delta \theta is an arbitrary variational field satisfying the essential boundary conditions. The body is approximated geometrically with finite elements, so the temperature is interpolated as
\theta = N ^ {N} (\mathbf {x}) \theta^ {N}, \quad N = 1, 2, \dots ,
where \theta ^ { N } are nodal temperatures. The Galerkin approach assumes that \delta \theta , the variational field, is interpolated by the same functions:
\delta \theta = N ^ {N} \delta \theta^ {N}.
First- and second-order polynomials in one, two, and three dimensions are used for the N ^ { N } . With these interpolations the variational statement, Equation 2.11.1-3, becomes
\begin{array}{l} \delta \theta^ {N} \left\{\int_ {V} N ^ {N} \rho \dot {U} d V + \int_ {V} \frac {\partial N ^ {N}}{\partial \mathbf {x}} \cdot \mathbf {k} \cdot \frac {\partial \theta}{\partial \mathbf {x}} d V \right. \\ = \int_ {V} N ^ {N} r d V + \int_ {S _ {q}} N ^ {N} q d S \Biggr \}, \\ \end{array}
and since the \delta \theta ^ { N } are arbitrarily chosen, this gives the system of equations
Equation 2.11.1-4
\begin{array}{l} \int_ {V} N ^ {N} \rho \dot {U} d V + \int_ {V} \frac {\partial N ^ {N}}{\partial \mathbf {x}} \cdot \mathbf {k} \cdot \frac {\partial \theta}{\partial \mathbf {x}} d V \\ = \int_ {V} N ^ {N} r d V + \int_ {S _ {q}} N ^ {N} q d S. \\ \end{array}
This set of equations is the "continuous time description" of the geometric approximation.
Time integration
ABAQUS/Standard uses the backward difference algorithm:
\dot {U} _ {t + \Delta t} = (U _ {t + \Delta t} - U _ {t}) (1 / \Delta t).
Equation 2.11.1-5
This operator is chosen for a number of reasons. First of all, we choose from one-step operators of the form
f _ {t + \Delta t} = f _ {t} + \big ((1 - \gamma) \dot {f} _ {t} + \gamma \dot {f} _ {t + \Delta t} \big) \Delta t
because of their simplicity in implementation (for example, no special starting procedures are needed) and well-understood behavior. For \gamma < 1 / 2 such operators are only conditionally stable for linear heat transfer problems. We prefer to work with unconditionally stable methods, because ABAQUS is most commonly applied to problems where the solution is sought over very long time periods (compared to the stability limit for the explicit form of the operator, \gamma = 0 ) , and so choose \gamma \geq 1 / 2 . Of these operators the central difference method, \gamma = 1 / 2 , has the highest accuracy. However, that form of the operator tends to produce oscillations in the early time solution that are not present in the backward difference form. Thus, we use \gamma = 1 \colon backward difference. Introducing the operator, Equation 2.11.1-5, into the energy balance Equation 2.11.1-4 gives
Equation 2.11.1-6
\begin{array}{l} \frac {1}{\Delta t} \int_ {V} N ^ {N} \rho (U _ {t + \Delta t} - U _ {t}) d V + \int_ {V} \frac {\partial N ^ {N}}{\partial \mathbf {x}} \cdot \mathbf {k} \cdot \frac {\partial \theta}{\partial \mathbf {x}} d V \\ - \int_ {V} N ^ {N} r d V - \int_ {S _ {q}} N ^ {N} q d S = 0. \\ \end{array}
This nonlinear system is solved by a modified Newton method. The method is modified Newton because the tangent matrix (the Jacobian matrix)--that is, the rate of change of the left-hand side of Equation 2.11.1-6 with respect to \theta _ { t + \Delta t } ^ { N } --is not formed exactly. The formation of the terms in this tangent matrix is now described.
The internal energy term gives a Jacobian contribution:
\frac {1}{\Delta t} \int_ {V} N ^ {N} \rho \frac {d U}{d \theta} \bigg | _ {t + \Delta t} N ^ {M} d V.
( d U / d \theta ) | _ { t + \Delta t } is the specific heat, c ( \theta ) , outside the latent heat range, and is c + L / ( \theta _ { L } - \theta _ { S } ) if \theta _ { L } > \theta _ { t + \Delta t } > \theta _ { S } at the integration point, where \theta _ { L } and \theta _ { S } are the liquidus and solidus temperatures and L is the latent heat associated with this phase change.
In severe latent heat cases this term can result in numerical instabilities, as the stiffness term d U / d \theta is small outside the solidus-liquidus temperature range and is very stiff inside that rather narrow range. To avoid such instabilities in those cases this term is modified to a secant term during the early iterations of the solution to a time step. Since the modification occurs only in cases involving latent heat, it affects only those problems.
Procedures
The conductivity term gives a Jacobian contribution:
\int_ {V} \frac {\partial N ^ {N}}{\partial \mathbf {x}} \cdot \mathbf {k} \Big | _ {t + \Delta t} \cdot \frac {\partial N ^ {M}}{\partial \mathbf {x}} d V + \int_ {V} \frac {\partial N ^ {N}}{\partial \mathbf {x}} \cdot \frac {\partial \mathbf {k}}{\partial \theta} \Big | _ {t + \Delta t} \cdot \frac {\partial \theta}{\partial \mathbf {x}} \Big | _ {t + \Delta t} N ^ {M} d V.
The second of these terms is typically small, since the conductivity usually varies only slowly with temperature. Because of this, and because the term is not symmetric, it is usually more efficient to omit it. This term is omitted unless the unsymmetric solver is chosen. Prescribed surface fluxes and body fluxes can also be temperature dependent and will then give rise to Jacobian contributions.
With film and radiation conditions, the surface flux term gives a Jacobian contribution:
\int_ {S} N ^ {N} \frac {\partial q}{\partial \theta} \Big | _ {t + \Delta t} N ^ {M} d S.
For film conditions, q = h ( \theta ) ( \theta - \theta ^ { o } ) ;
\frac {\partial q}{\partial \theta} = \frac {\partial h}{\partial \theta} (\theta - \theta^ {o}) + h,
while for radiation, q = A ( \theta ^ { 4 } - \theta ^ { o 4 } ) ;
\frac {\partial q}{\partial \theta} = 4 A \theta^ {3}.
These terms are included in exactly this form in the Jacobian. The modified Newton method is then
Equation 2.11.1-7
\begin{array}{l} \left[ \frac {1}{\Delta t} \int_ {V} N ^ {N} \rho \frac {d U}{d \theta} \right| _ {t + \Delta t} N ^ {M} d V + \int_ {V} \frac {\partial N ^ {N}}{\partial \mathbf {x}} \cdot \mathbf {k} \big | _ {t + \Delta t} \cdot \frac {\partial N ^ {M}}{\partial \mathbf {x}} d V \\ + \int_ {S} N ^ {N} \left(\frac {\partial h}{\partial \theta} (\theta - \theta^ {o}) + h + 4 A \theta^ {3}\right) \left. N ^ {M} d S \right] \overline {{{c}}} ^ {M} \\ = \int_ {V} N ^ {N} r d V + \int_ {S _ {q}} N ^ {N} q d S - \frac {1}{\Delta t} \int_ {V} N ^ {N} \rho (U _ {t + \Delta t} - U _ {t}) d V \\ - \int_ {V} \frac {\partial N ^ {N}}{\partial \mathbf {x}} \cdot \mathbf {k} \cdot \frac {\partial \theta}{\partial \mathbf {x}} d V, \\ \end{array}
\theta _ { t + \Delta t , i + 1 } ^ { N } = \theta _ { t + \Delta t , i } ^ { N } + \overline { { c } } ^ { N } ; i = iteration number.
For purely linear systems Equation 2.11.1-7 is linear in \overline { { c } } ^ { M } and, hence, in \theta _ { t + \Delta t } ^ { N } , so a single equation solution provides the \theta _ { t + \Delta t } ^ { N } . Since the method usually is only a minor modification of Newton's method, convergence is rapid.
ABAQUS/Standard uses an automatic (self-adaptive) time stepping algorithm to choose \Delta t . . This is based on a user-supplied tolerance on the maximum temperature change allowed in a time increment, and the increment is adjusted according to this parameter, as well as the convergence rate of Equation
2.11.1-7 in nonlinear cases.
The first-order heat transfer elements (such as 2-node link, 4-node quadrilateral, and 8-node brick) use a numerical integration rule with the integration stations located at the corners of the element for the heat capacitance terms. This means that the Jacobian term associated with the internal energy rate is diagonal. This approach is especially effective when strong latent heat effects are present. The second-order elements use conventional Gaussian integration. Thus, second-order elements are to be preferred for problems when the solution will be smooth (without latent heat effects), whereas the first-order elements should be used in nonsmooth cases (with latent heat).
The HEATCAP element is available for modeling lumped heat capacitance at a point. The associated concentrated film and concentrated radiation loading options are specified with the *CFILM and *CRADIATE options, respectively. These features are also allowed in coupled temperature-displacement and coupled thermal-electrical analysis.
2.11.2 Shell heat conduction
This section describes the formulation used in the shell heat conduction elements in ABAQUS/Standard. The basis of the elements is a combination of piecewise quadratic interpolation of temperature through the thickness of the shell and either linear interpolation (in elements DS3 and DS4) or quadratic interpolation (in elements DS6 and DS8) on the reference surface of the shell. The isoparametric interpolation functions for the shell reference surface are identical in form to those used for the solid quadrilateral and triangular elements and can be found in ``Solid isoparametric quadrilaterals and hexahedra,'' Section 3.2.4, and ``Triangular, tetrahedral, and wedge elements,'' Section 3.2.6, respectively. Nodal temperature values are stored at a set of points through the thickness (points P below) at each node of the element (nodes N below). For the purpose of numerical integration of the finite element equations, \mathbf { a } \ 2 \times 2 Gauss integration scheme with a 2 \times 2 nodal integration scheme for the internal energy and specific heat term is used for the quadrilateral element DS4 and a 3 \times 3 Gauss integration scheme is used for the quadrilateral element DS8. Three- and six-point integration schemes are used for the triangular elements DS3 and DS6, respectively, the details of which can be found in ``Triangular, tetrahedral, and wedge elements,'' Section 3.2.6.
Let \left( \phi _ { 1 } , \phi _ { 2 } \right) be material coordinates of a point in the reference surface of the shell, and let s _ { 3 } measure position through the thickness of the shell so that - h / 2 - z _ { 0 } \leq s _ { 3 } \leq h / 2 - z 0 , where h is the thickness of the shell, z _ { 0 } is the offset of the reference surface from the midsurface as discussed in ``Transverse shear stiffness in composite shells and offsets from the midsurface, '' Section 3.6.8. The position of any point in the shell is given by
\mathbf {x} = \mathbf {x} ^ {o} (\phi_ {1}, \phi_ {2}) + s _ {3} \mathbf {n} (\phi_ {1}, \phi_ {2}),
where
xo \mathbf { x } ^ { o }
is the position of a point in the reference surface, and
n
is the unit normal to the reference surface of the shell.
The temperature interpolation can be written as
\theta = N ^ {N} (\phi_ {1}, \phi_ {2}) M ^ {P} (s _ {3}) \overline {{\theta}} ^ {N P},
Equation 2.11.2-1
where
M ^ {P} (s _ {3})
is a piecewise parabolic interpolation,
N ^ {N} (\phi_ {1}, \phi_ {2})
is an interpolator in the reference surface, and
\overline {{\theta}} ^ {N P}
are nodal temperature values (at node N , point P through the thickness).
The basic heat energy balance is Equation 2.11.1-3, with the approximate Jacobian matrix for the Newton method based on
\int_ {V} \rho \delta \theta \frac {d \dot {U}}{d \theta_ {t + \Delta t}} d \theta d V + \int_ {V} \frac {\partial \delta \theta}{\partial \mathbf {x}} \cdot \mathbf {k} \cdot \frac {\partial d \theta}{\partial \mathbf {x}} d V - \int_ {S} \delta \theta \frac {\partial q}{\partial \theta_ {t + \Delta t}} d \theta d S,
where d \theta is the correction to the temperature solution at time t + \Delta t . The derivation of this form is discussed in ``Uncoupled heat transfer analysis,'' Section 2.11.1. The form of these terms for shell heat conduction elements is obtained by introducing the interpolator, Equation 2.11.2-1, and neglecting the change in area, with respect to s _ { 3 } , of surfaces parallel to the reference surface.
The internal energy rate term (the first term in Equation 2.11.1-3) contributes, to the residual,
\int_ {A} \left[ N ^ {N} \int_ {- h / 2} ^ {h / 2} M ^ {P} \rho \dot {U} _ {t + \Delta t} d s _ {3} \right] d A,
and to the Jacobian,
\int_ {A} \left[ N ^ {N} \int_ {- h / 2} ^ {h / 2} M ^ {P} \rho \frac {d \dot {U}}{d \theta_ {t + \Delta t}} M ^ {Q} d s _ {3} \right] N ^ {M} d A.
For the second term the temperature derivatives are taken with respect to a local orthogonal system \left( { { s _ { 1 } } , { s _ { 2 } } , { s _ { 3 } } } \right) , where s _ { 1 } and s _ { 2 } measure distance along local base vectors \mathbf { e } _ { 1 } and \mathbf { e } _ { 2 } , in the reference surface of the shell, set up according to the standard convention in ABAQUS for such local systems in shells. The term is formed by first introducing an intermediate set of temperature values,
Procedures
\left\lfloor \theta^ {P}, \frac {\partial \theta^ {P}}{\partial s _ {1}}, \frac {\partial \theta^ {P}}{\partial s _ {2}} \right\rfloor ,
corresponding to each temperature value point through the thickness, at each section where integration through the thickness is performed. Since the number of temperature values on the section is the same as the number of integration points, M ^ { P } is unity in the appropriate locations and zero everywhere else. Then we can interpolate to the section by
\left\{ \begin{array}{l} \theta^ {P} \\ \frac {\partial \theta^ {P}}{\partial s _ {1}} \\ \frac {\partial \theta^ {P}}{\partial s _ {2}} \end{array} \right\} = \left\{\beta^ {M} \right\} \overline {{\theta}} ^ {M P},
where
\left\{\beta^ {M} \right\} = \left\{ \begin{array}{l} N ^ {M} \\ \frac {\partial N ^ {M}}{\partial s _ {1}} \\ \frac {\partial N ^ {M}}{\partial s _ {2}} \end{array} \right\}.
The piecewise quadratic interpolation through the thickness then gives
\left\{ \begin{array}{l} \frac {\partial \theta}{\partial s _ {1}} \\ \frac {\partial \theta}{\partial s _ {2}} \\ \frac {\partial \theta}{\partial s _ {3}} \end{array} \right\} = \left[ \gamma^ {P} \right] \left\{ \begin{array}{l} \theta^ {P} \\ \frac {\partial \theta^ {P}}{\partial s _ {1}} \\ \frac {\partial \theta^ {P}}{\partial s _ {2}} \end{array} \right\},
where
\left[ \gamma^ {P} \right] = \left[ \begin{array}{c c c} 0 & M ^ {P} & 0 \\ 0 & 0 & M ^ {P} \\ d M ^ {P} / d s _ {3} & 0 & 0 \end{array} \right].
The conductivity term in the Jacobian then is
\int_ {A} \left\lfloor \beta^ {N} \right\rfloor \left[ \int_ {- h / 2} ^ {h / 2} \left[ \gamma^ {P} \right] ^ {T} [ k ] \left[ \gamma^ {Q} \right] d s _ {3} \right] \left\{\beta^ {M} \right\} d A,
where [ k ] is the local conductivity matrix and the same term multiplied by \theta _ { t + \Delta t } ^ { M Q } appears in the residual.