김경종
22 KiB
Interface Modeling
fluid, evaluated at the slave node. The summation extends over all master nodes i in the vicinity of the slave node projection p \left( { \bf x } _ { \mathrm { N } } \right) (see Figure 5.2.7-1). The computation is repeated for each slave node \mathbf { x } _ { \mathrm { N } } on the surface S _ { \mathrm { f s } } \cup S _ { \mathrm { f r s } } and assembled to form the entire coupling matrix.
Similarly, the contribution to the pressure coupling term in the structural equation due to slave node \mathbf { x } _ { \mathrm { N } } is approximated by
\int_ {S _ {\mathrm{fs}} \cup S _ {\mathrm{frs}}} \delta \mathbf {u} ^ {m} \cdot \mathbf {n} ^ {-} p d S \approx p _ {N} A _ {N} \left[ \sum_ {i} \mathbf {n} ^ {-} \left(\mathbf {x} _ {\mathrm{N}}\right) \cdot \mathrm{N} ^ {i} \left(\mathbf {p} \left(\mathbf {x} _ {\mathrm{N}}\right)\right) \right],
where p _ { N } is the acoustic pressure at slave node \mathbf { x } _ { \mathrm { N } } , and the sum is again over the master nodes i in the vicinity of the slave node projection.
These expressions for the coupling terms result in matrices that are the transpose of each other. The normal vectors at the slave surface are used, so these vectors must be well-defined (see ``Surfaces: overview,'' Section 2.3.1 of the ABAQUS/Standard User's Manual).
Fluid master, solid/structural slave
If the solid medium is designated as the slave, the values on this surface are constrained to equal values interpolated from the master surface. The pointwise fluid-solid coupling condition is again enforced at the slave nodes, resulting in acoustic pressure degrees of freedom added to the solid slave surface (see Figure 5.2.7-2). These slave pressures are constrained by the master surface acoustic pressures and eliminated; the slave displacements are not constrained directly.
Figure 5.2.7-2 Fluid master.
text_image
P(x_N) i+1 p_{i+1} fluid master surface u_2 x_N x_{N+1} solid slave surface x_{N-1} p u_1 A_N n(x_N)
Hence, the contribution to the coupling term in the acoustic equation for a single slave node \mathbf { x } _ { \mathrm { N } } is approximated by
\int_ {S _ {\mathrm{fs}} \cup S _ {\mathrm{frs}}} \delta p \mathbf {n} ^ {-} \cdot \ddot {\mathbf {u}} ^ {m} d S \approx \left[ A _ {N} \mathbf {n} ^ {-} \left(\mathbf {x} _ {\mathrm{N}}\right) \cdot \ddot {\mathbf {u}} _ {N} ^ {m} \right] \left[ \sum_ {i} \mathrm{H} ^ {i} \left(\mathbf {p} \left(\mathbf {x} _ {\mathrm{N}}\right)\right) \right],
where \ddot { \mathbf { u } } _ { N } ^ { m } is the structural acceleration at the slave node and \mathrm { H } ^ { i } \left( \mathbf { p } \left( \mathbf { x } _ { \mathrm { N } } \right) \right) are the interpolants on the fluid (master) surface evaluated at the projection { \bf p } \left( { \bf x } _ { \mathrm { N } } \right) . The summation extends over the master nodes i in the vicinity of the slave node projection and is repeated for all solid/structural nodes in the
Interface Modeling
surface S _ { \mathrm { f s } } \cup S _ { \mathrm { f r s } } to compute the entire coupling matrix.
The contribution of a slave node to the coupling term in the structural equation is approximated by
\int_ {S _ {\mathrm{fs}} \cup S _ {\mathrm{frs}}} \delta \mathbf {u} ^ {m} \cdot \mathbf {n} ^ {-} p d S \approx A _ {N} \sum_ {i} \mathrm{H} ^ {i} \left(\mathbf {p} \left(\mathbf {x} _ {\mathrm{N}}\right)\right) p _ {i},
where p _ { i } is the pressure at master node i and the sum is again over the master nodes in the vicinity of the slave node projection p (xN ).
6. Loading and Constraints
6.1 Dynamic loading
6.1.1 Centrifugal, Coriolis, and rotary acceleration forces
Many of the elements in ABAQUS allow centrifugal, Coriolis, and rotary acceleration forces to be included. This section defines these load types.
It is assumed that the model (or that part of it to which these forces are applied) is described in a coordinate system that is rotating with an angular velocity, !, and/or an angular (rotary) acceleration, d \omega / d t . Let \mathbf { A } = \left[ \mathbf { a } _ { 1 } , \mathbf { a } _ { 2 } , \mathbf { a } _ { 3 } \right] be a right-handed set of unit, orthogonal vectors that form a basis for this system. Then, d { \bf A } / d t = { \pmb { \omega } } \times { \bf A } and d ^ { 2 } { \bf A } / d t ^ { 2 } = d \pmb { \omega } / d t \times { \bf A } + \pmb { \omega } \times d { \bf A } / d t .
If the angular velocity is cast as
\boldsymbol {\omega} = \omega \mathbf {n},
where \omega is the magnitude of ! and n is the unit axis of rotation, the rotary acceleration is
\frac {d \pmb {\omega}}{d t} = \frac {d \omega}{d t} \mathbf {n} + \omega \frac {d \mathbf {n}}{d t} = \alpha \mathbf {n} ^ {R O T A},
where the term !dn=dt represents the effect of the motion of the axis of rotation (precessional motion);
\alpha = d \pmb { \omega } / d t \cdot \mathbf { n } ^ { R O T A } is the magnitude of the rotary acceleration; and \mathbf { n } ^ { R O T A } is the axis of rotary acceleration. If d \mathbf { n } / d t = \mathbf { 0 } , then \alpha = d \omega / d t and \mathbf { n } ^ { R O T A } = \mathbf { n } . In component form
\omega^ {k} = \omega n ^ {k}
and
\frac {d \omega^ {k}}{d t} = \frac {d \omega}{d t} n ^ {k} + \omega \frac {d n ^ {k}}{d t}.
Let \mathbf { x } _ { 0 } be a point on the axis of rotation. The position of a material particle, \mathbf { x } , can be written
\mathbf {x} = \mathbf {x} _ {0} + y ^ {i} \mathbf {a} _ {i},
where y ^ { i } , i = 1 , 2 , 3 , are the coordinates of the point in the rotating basis system. Taking time derivatives,
{\frac {d \mathbf {x}}{d t}} = {\frac {d \mathbf {x} _ {0}}{d t}} + {\frac {d y ^ {i}}{d t}} \mathbf {a} _ {i} + y ^ {i} {\pmb {\omega}} \times \mathbf {a} _ {i}
Loading and Constraints
and
\frac {d ^ {2} \mathbf {x}}{d t ^ {2}} = \frac {d ^ {2} \mathbf {x} _ {0}}{d t ^ {2}} + \frac {d ^ {2} y ^ {i}}{d t ^ {2}} \mathbf {a} _ {i} + 2 \frac {d y ^ {i}}{d t} \pmb {\omega} \times \mathbf {a} _ {i} + y ^ {i} \pmb {\omega} \times (\pmb {\omega} \times \mathbf {a} _ {i}) + y ^ {i} \frac {d \pmb {\omega}}{d t} \times \mathbf {a} _ {i}.
We assume that the origin of the rotating system, x _ { 0 } , is fixed, so that
\frac {d \mathbf {x} _ {0}}{d t} = \frac {d ^ {2} \mathbf {x} _ {0}}{d t ^ {2}} = 0.
With this simplification
{\frac {d \mathbf {x}}{d t}} = {\frac {d y ^ {i}}{d t}} \mathbf {a} _ {i} + y ^ {i} \pmb {\omega} \times \mathbf {a} _ {i}
and
Equation 6.1.1-1
\frac {d ^ {2} \mathbf {x}}{d t ^ {2}} = \frac {d ^ {2} y ^ {i}}{d t ^ {2}} \mathbf {a} _ {i} + 2 \frac {d y ^ {i}}{d t} \pmb {\omega} \times \mathbf {a} _ {i} + y ^ {i} \pmb {\omega} \times (\pmb {\omega} \times \mathbf {a} _ {i}) + y ^ {i} \frac {d \pmb {\omega}}{d t} \times \mathbf {a} _ {i}.
The virtual work contribution from the d'Alembert forces is
\delta W ^ {A} = - \int_ {V ^ {0}} \rho^ {0} \frac {d ^ {2} \mathbf {x}}{d t ^ {2}} \cdot \delta \mathbf {x} d V ^ {0},
where \rho ^ { 0 } is the mass density of the body in its reference configuration, where its volume is V ^ { 0 } and ±x is a virtual variational field. For the part of the body described in the rotating system, the acceleration, d ^ { 2 } { \bf x } / d t ^ { 2 } , is given by Equation 6.1.1-1, while \delta \mathbf { x } = \delta y ^ { i } \mathbf { a } _ { i } only, since ! and d!=dt are prescribed and \mathbf { x } _ { 0 } is fixed. Thus,
\delta W ^ {A} = - \int_ {V ^ {0}} \rho^ {0} \left[ \frac {d ^ {2} y ^ {i}}{d t ^ {2}} \delta y ^ {i} + \left\{2 \frac {d y ^ {i}}{d t} \pmb {\omega} \times \mathbf {a} _ {i} + y ^ {i} \pmb {\omega} \times (\pmb {\omega} \times \mathbf {a} _ {i}) + y ^ {i} \frac {d \pmb {\omega}}{d t} \times \mathbf {a} _ {i} \right\} \cdot \delta y ^ {j} \mathbf {a} _ {j} \right] d V ^ {0}.
Simplifying,
\begin{array}{l} \delta W ^ {A} = - \int_ {V ^ {0}} \rho^ {0} \bigg [ \frac {d ^ {2} y ^ {i}}{d t ^ {2}} \delta y ^ {i} + 2 \frac {d y ^ {i}}{d t} \delta y ^ {j} \pmb {\omega} \cdot (\mathbf {a} _ {i} \times \mathbf {a} _ {j}) \\ \left. + y ^ {i} \delta y ^ {j} \left(\pmb {\omega} \cdot \mathbf {a} _ {i} \pmb {\omega} \cdot \mathbf {a} _ {j} - (\pmb {\omega} \cdot \pmb {\omega}) \delta_ {i j}\right) + y ^ {i} \delta y ^ {j} \frac {d \pmb {\omega}}{d t} \times \mathbf {a} _ {i} \cdot \mathbf {a} _ {j} \right] d V ^ {0}. \\ \end{array}
The terms in \delta W ^ { A } can be identified as follows. The first term,
Loading and Constraints
- \int_ {V ^ {0}} \rho^ {0} \frac {d ^ {2} y ^ {i}}{d t ^ {2}} \delta y ^ {i} d V ^ {0},
is the usual "consistent mass matrix" term associated with acceleration of the material particles with respect to the rotating system.
Writing the angular velocity of the rotating basis system as its components in that system, \omega = \omega ^ { i } \mathbf { a } _ { i } , gives the second term as
- \int_ {V ^ {0}} 2 \rho^ {0} \epsilon_ {i j k} \frac {d y ^ {i}}{d t} \delta y ^ {j} \omega^ {k} d V ^ {0},
where \epsilon _ { i j k } is the alternator tensor. This term is the Coriolis force term and arises whenever there is velocity in the rotating system, which can happen in dynamic analysis or in quasi-static cases in which a constant velocity has been introduced (for example, by the *INITIAL CONDITIONS option).
The third term is, likewise, rewritten as
- \int_ {V ^ {0}} \rho^ {0} y ^ {i} \delta y ^ {j} (\omega^ {i} \omega^ {j} - \omega^ {k} \omega^ {k} \delta_ {i j}) d V ^ {0}.
This term is the centrifugal load term.
Similarly, the fourth term is rewritten as
- \int_ {V ^ {0}} \rho^ {0} \epsilon_ {i j k} y ^ {i} \delta y ^ {j} \frac {d \omega^ {k}}{d t} d V ^ {0}.
This term is the rotary acceleration load term.
In ABAQUS/Standard the centrifugal load, Coriolis, and rotary acceleration terms contribute to the "load stiffness matrix." The centrifugal load term has a symmetric load stiffness matrix,
- \int_ {V ^ {0}} \rho^ {0} d y ^ {i} \delta y ^ {j} (\omega^ {i} \omega^ {j} - \omega^ {k} \omega^ {k} \delta_ {i j}) d V ^ {0};
the Coriolis term has an antisymmetric "load damping matrix,"
\int_ {V ^ {0}} 2 \rho^ {0} \epsilon_ {i j k} \frac {\partial}{\partial y ^ {l}} \left(\frac {d y ^ {i}}{d t}\right) d y ^ {l} \delta y ^ {j} \omega^ {k} d V ^ {0};
and the rotary acceleration term has an antisymmetric load stiffness matrix,
- \int_ {V ^ {0}} \rho^ {0} \epsilon_ {i j k} d y ^ {i} \delta y ^ {j} \frac {d \omega^ {k}}{d t} d V ^ {0}.
6.1.2 Baseline correction of accelerograms
ABAQUS/Standard offers baseline correction of acceleration records for time domain analysis. The correction technique is that proposed by Newmark (1973). An acceleration correction, a _ { 0 } ( t ) , is added to the raw data record, a ( t ) , to produce a corrected acceleration record, \textstyle a _ { c } ( t ) = a + a _ { 0 } , such that the mean square velocity over the time of the event is minimized. The acceleration correction is parabolic over any number of time intervals during the event:
a _ {0} (t) = C _ {1} + C _ {2} \left(\frac {t - T _ {1}}{T _ {2}}\right) + C _ {3} \left(\frac {t - T _ {1}}{T _ {2}}\right) ^ {2}, \qquad T _ {1} < t < T _ {2},
where T _ { 1 } , T _ { 2 } denote the limits of a time interval and C _ { k , } k = 1 , 2 , 3 , are constants obtained from the velocity minimization:
\frac {\partial}{\partial C _ {k}} \int_ {T _ {1}} ^ {T _ {2}} [ v _ {c} (t) ] ^ {2} d t = 0,
where v _ { c } ( t ) is the corrected velocity record obtained by integrating the corrected acceleration record, a _ { c } ( t ) .
It can be shown that the C _ { k } are defined for each time interval by
\left\{ \begin{array}{l} C _ {1} \\ C _ {2} \\ C _ {3} \end{array} \right\} = \left[ \begin{array}{c c c} - 3 0 0. & 9 0 0. & - 6 3 0. \\ 1 8 0 0. / \xi & - 5 7 6 0. / \xi & 4 2 0 0. / \xi \\ - 1 8 9 0. / \xi^ {2} & 6 3 0 0. / \xi^ {2} & - 4 7 2 5. / \xi^ {2} \end{array} \right] \left\{ \begin{array}{l} A _ {1} + \left(v _ {c} (T _ {1}) - v (T _ {1})\right) / (2 \Delta T) \\ A _ {2} + \left(v _ {c} (T _ {1}) - v (T _ {1})\right) / (3 \Delta T) \\ A _ {3} + \left(v _ {c} (T _ {1}) - v (T _ {1})\right) / (4 \Delta T) \end{array} \right\},
where \xi = \Delta T / T _ { 2 } ; \Delta T = T _ { 2 } - T _ { 1 } ; and A _ { 1 } , A _ { 2 } , and A _ { 3 } are defined as
A _ {1} = \frac {1}{(\Delta T) ^ {3}} \int_ {0} ^ {\Delta T} v (\tau + T _ {1}) \tau d \tau ,
A _ {2} = \frac {1}{(\Delta T) ^ {4}} \int_ {0} ^ {\Delta T} v (\tau + T _ {1}) \tau^ {2} d \tau ,
A _ {3} = \frac {1}{(\Delta T) ^ {5}} \int_ {0} ^ {\Delta T} v (\tau + T _ {1}) \tau^ {3} d \tau .
Here v ( t ) denotes the uncorrected velocity record obtained by integration of the uncorrected acceleration record, a ( t ) . These velocities are obtained by assuming the uncorrected and the corrected acceleration vary linearly over each time increment of the original acceleration history. This is not exact for the corrected acceleration record (because of the parabolic variation of the correction in time), but it is assumed that the acceleration history is discretized at sufficiently small time increments so that this is an insignificant error.
6.2 *AQUA loading
6.2.1 Drag, inertia, and buoyancy loading
For beam and truss structures immersed in fluid (e.g., offshore piping and riser problems),
ABAQUS/Standard provides a capability for introducing drag forces via Morison's equation, inertia loads, and buoyancy loads. Fluid drag is associated with velocities due to steady currents and any waves that may have been specified. Fluid inertia is associated with wave accelerations. Buoyancy has two components: the hydrostatic pressure measured to the mean fluid level and the dynamic pressure caused by the presence of waves. Partial submergence is done automatically for all fluid load types.
Drag and inertia loads are considered in two forms: distributed loads along the length of the element (distributed drag loading is further divided into a component normal to the element's axis and a component along the tangent to the element) and point drag and inertia loads where the beam changes cross-section.
Buoyancy loading is applied with a "closed-end" assumption; that is, it is assumed that the element's ends can support buoyancy loading normal to the element's cross-section. If the ends of the element are actually "open ended"--that is, the element's ends cannot support fluid pressure loads--point buoyancy forces are provided to remove the buoyancy forces at the ends of the element.
This section documents the form of these loadings. It is assumed that the fluid particle velocities and accelerations are known as functions of the current spatial location; they are defined by the steady current input, with the wave effect superposed, if the *WAVE option is used.
Continuously distributed drag and inertia loading
Distributed load types "FDD," "WDD," "FDT," and "FI" on beam, truss, or rigid beam elements provide continuously distributed drag and inertia via Morison's equations. To specify this loading, the following definitions are used:
\mathbf{v}_f is the fluid particle velocity (defined by the steady current input and possibly additional contributions from *WAVE definitions),
\mathbf{a}_f is the fluid particle acceleration when the *\mathrm{WAVE} option is used,
v_{p} is the velocity of a point on the element (nonzero during * DYNAMIC steps only),
a_{p} is the acceleration of a point on the element (nonzero during * DYNAMIC steps only),
\Delta \mathbf{v} = \mathbf{v}_f - \alpha_R\mathbf{v}_p is the relative fluid velocity,
t is a unit vector defining the axial direction at a point in the element,
\Delta \mathbf{v}_t = (\Delta \mathbf{v}\cdot \mathbf{t})\mathbf{t} is the relative tangential (axial) velocity of the fluid,
Loading and Constraints
\Delta \mathbf{v}_n = \Delta \mathbf{v} - \Delta \mathbf{v}_t is the relative transverse velocity of the fluid, C_T is the tangential (axial) drag coefficient, C_D is the transverse drag coefficient, C_M is the transverse inertia coefficient, C_A is the transverse added mass coefficient, \rho is the fluid density, \alpha_R is the structural velocity factor, and h is the exponent for tangential drag.
Then, the transverse drag force per unit length on the member is
\mathbf {F} _ {D} = \frac {1}{2} \rho C _ {D} D \Delta \mathbf {v} _ {n} \big (\Delta \mathbf {v} _ {n} \cdot \Delta \mathbf {v} _ {n} \big) ^ {\frac {1}{2}}.
The tangential drag force per unit length is
\mathbf {F} _ {T} = \frac {1}{2} \rho C _ {T} \pi D \Delta \mathbf {v} _ {t} | \Delta \mathbf {v} _ {t} | ^ {h - 1},
and the inertia force per unit length is
\mathbf {F} _ {I} = \frac {1}{4} \rho \pi D ^ {2} \big [ C _ {M} (\mathbf {a} _ {f} - \mathbf {a} _ {f} \cdot \mathbf {t t}) + C _ {A} (\mathbf {a} _ {p} - \mathbf {a} _ {p} \cdot \mathbf {t t}) \big ].
Only transverse drag is implemented for wind loading.
Point drag and inertia loading
At points where the section changes size and, thus, exposes an end area to the fluid (see Figure 6.2.1-1), additional drag and inertia forces arise.
Figure 6.2.1-1 Change in section of an immersed beam.
text_image
t Exposed end area, Δ A
For this case we need the additional definitions:
¢A
is the change in area of the section,
\mathbf { t }
is an outward normal to the exposed area (see Figure 6.2.1-1),
C _ { n }
is the drag coefficient associated with the discontinuity,
K _ { t s }
is the tangential inertia coefficient associated with the discontinuity,
L _ { t s } Lts
is the tangential added mass coefficient associated with the discontinuity, and
F _ { 1 s } , F _ { 2 s }
are fluid and structural acceleration shape factors for the tangential inertia term.
Then the drag force on the transition section is
\mathbf {F} _ {d s} = \left\{ \begin{array}{r l} \frac {1}{2} \rho C _ {n} \Delta A \Delta \mathbf {v} _ {f t} | \Delta \mathbf {v} _ {f t} | & \mathrm{for} \quad \mathbf {t} \cdot \Delta \mathbf {v} _ {f t} \leq 0, \\ \mathbf {0} & \mathrm{for} \quad \mathbf {t} \cdot \Delta \mathbf {v} _ {f t} > 0. \end{array} \right.
This force is nonzero only when the relative velocity of the fluid has a negative projection on the outward normal; i.e., when the fluid is flowing against the exposed surface. The inertia force is
\mathbf {F} _ {I s} = \rho K _ {t s} F _ {1 s} \mathbf {a} _ {f} \cdot \mathbf {t t} - \rho L _ {t s} F _ {2 s} \mathbf {a} _ {p} \cdot \mathbf {t t}.
Distributed buoyancy
ABAQUS assumes closed-end conditions when computing the distributed buoyancy loads (load type PB) for beams, pipes, rigid beams, and elbows. An open-end condition can be simulated by using load type TSB to remove the buoyancy load from the end of the element.
Loading and Constraints
For truss elements the effect of buoyancy is simply Archimedes' principle; that is, a vertical force equal to the weight of the displaced fluid is applied to the element.
For buoyancy loading the beam is assumed to have a uniform circular pipe section with the diameter specified by the user on the data lines following the *DLOAD option.
Consider a pressure field that varies with z, the vertical coordinate, where z = \mathbf { k } \cdot \mathbf { x } and k is the unit vector pointing in the vertical direction. Then
\hat {p} = p (z).
For hydrostatic pressure the dependence on the vertical coordinate is linear in z , ,
\hat {p} = \rho g (z _ {0} - z),
where z _ { 0 } is the vertical location of the free surface of the fluid (either the free surface elevation of the fluid inside the pipe or the elevation of the mean water level for external pressure), \rho is the density of the fluid, and g is the acceleration due to gravity. If the *WAVE option is used to define a wave field, the wave field generates a dynamic pressure effect. In this case the pressure field has a nonlinear variation with respect to location. Under the assumption of waves with small amplitude relative to the wave length and depth of the sea, we can assume that the dynamic pressure is slowly varying with respect to the wave direction; hence, derivatives of the pressure with respect to the coordinates in the plane parallel to the still water surface are neglected. However, the dynamic pressure field has a nonlinear dependence on z.
The pressure field results in a nodal loading contribution that can be written as follows. Let \delta \Pi _ { p } be the pressure loading contribution to the weak form of the equilibrium equations. This contribution can be written in terms of a nodal load vector, pN , as
\delta \Pi_ {p} = \mathbf {p} _ {N} \cdot \delta \mathbf {x} _ {N},
where the nodal load vector is
\mathbf {p} _ {N} = \int_ {0} ^ {l} \pi R ^ {2} \left[ \hat {p} \frac {d \mathbf {x}}{d s} \frac {d N _ {N}}{d s} + \frac {\partial \hat {p}}{\partial z} N _ {N} \mathbf {k} \right] d s,
where R is the radius of the pipe (internal or external), \hat { p } is the total pressure (including hydrostatic and dynamic pressure), N _ { N } are the shape functions associated with the nodes on the element, and l is the current length of the element.
Point buoyancy
Point buoyancy is needed to apply discrete buoyant forces on an exposed surface. Point buoyancy can be used to remove the closed-end loading condition on beam, pipe, or rigid beam elements by applying a negative load on the exposed area. Use the following definitions: