Files
김경종 b7f84e1c0f
Tests / Hermetic test suite (push) Has been cancelled
Tests / Skill frontmatter validation (push) Has been cancelled
add documents
2026-05-29 15:59:56 +09:00

422 lines
25 KiB
Markdown
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<!-- source-page: 701 -->
On combining the lumped-mass matrices Eqs. (16.2.12) and (16.6.1) for the bar and beam, respectively, the resulting local axes plane frame lumped-mass matrix is
$$
\underline {{{\hat {m}}}} = \frac {\rho A L}{2} \left[ \begin{array}{c c c c c c} \hat {d} _ {1 x} & \hat {d} _ {1 y} & \hat {\phi} _ {1} & \hat {d} _ {2 x} & \hat {d} _ {2 y} & \hat {\phi} _ {2} \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{array} \right] \tag {16.7.7}
$$
The global mass matrix m for a plane frame element arbitrarily oriented in x-y coordinates is transformed according to Eq. (16.7.1), where the transformation matrix T is now given by Eq. (5.1.10) and either Eq. (16.7.6) for consistent-mass or (16.7.7) for lumped-mass matrices.
Because a longhand solution of the time-dependent plane frame problem is quite lengthy, only a computer program solution will be presented in Section 16.9.
# Plane Stress=Strain Element
The plane stress, plane strain, constant-strain triangle element (Figure 1621) consistent-mass matrix is obtained by using the shape functions from Eq. (6.2.18) and the shape function matrix given by substituting
$$
\underline {{N}} = \left[ \begin{array}{c c c c c c} N _ {1} & 0 & N _ {2} & 0 & N _ {3} & 0 \\ 0 & N _ {1} & 0 & N _ {2} & 0 & N _ {3} \end{array} \right]
$$
into Eq. (16.2.19) to obtain
$$
[ m ] = \rho \int_ {V} \underline {{N}} ^ {T} \underline {{N}} d V \tag {16.7.8}
$$
![](images/page-701_e80a58c9aeaca469b553c3c0f452593315007fef96174bff2dd2e7e4cc967a34.jpg)
<details>
<summary>text_image</summary>
y
v₃
3
u₃
v₂
v₁
1
u₁
2
u₂
x
</details>
Figure 1621 CST element with nodal degrees of freedom
<!-- source-page: 702 -->
Letting $d V = t d A$ and noting that $\begin{array} { r } { \int _ { A } N _ { 1 } ^ { 2 } d A = \frac { 1 } { 6 } A , \int _ { A } N _ { 1 } N _ { 2 } d A = \frac { 1 } { 1 2 } A } \end{array}$ , and so on, we obtain the CST global consistent-mass matrix as
$$
\underline {{m}} = \frac {\rho t A}{1 2} \left[ \begin{array}{c c c c c c} 2 & 0 & 1 & 0 & 1 & 0 \\ & 2 & 0 & 1 & 0 & 1 \\ & & 2 & 0 & 1 & 0 \\ & & & 2 & 0 & 1 \\ & & & & 2 & 0 \\ \text { Symmetry } & & & & & 2 \end{array} \right] \tag {16.7.9}
$$
For the isoparametric quadrilateral element for plane stress and plane strain considered in Chapter 10, we use the shape functions given by Eq. (10.3.5) with the shape function matrix given in Eq. (10.3.4) substituted into Eq. (16.7.10). This yields the quadrilateral element consistent-mass matrix as
$$
\underline {{m}} = \rho t \int_ {- 1} ^ {1} \int_ {- 1} ^ {1} \underline {{N}} ^ {T} \underline {{N}} \det \underline {{J}} d s d t \tag {16.7.10}
$$
The integral in Eq. (16.7.10) is evaluated best by numerical integration as described in Section 10.5.
# Axisymmetric Element
The axisymmetric triangular element (considered in Chapter 9 and shown in Figure 1622) consistent-mass matrix is given by
$$
\underline {{m}} = \int_ {V} \rho \underline {{N}} ^ {T} \underline {{N}} d V = \int_ {A} \rho \underline {{N}} ^ {T} \underline {{N}} 2 \pi r d A \tag {16.7.11}
$$
Since $r = N _ { 1 } r _ { 1 } + N _ { 2 } r _ { 2 } + N _ { 3 } r _ { 3 }$ , we have
$$
\underline {{m}} = 2 \pi \rho \int_ {A} (N _ {1} r _ {1} + N _ {2} r _ {2} + N _ {3} r _ {3}) \underline {{N}} ^ {T} \underline {{N}} d A \tag {16.7.12}
$$
![](images/page-702_7eaa158d505530fdc803a8b9722fe5780e5f58fbb643a007639d77a112cc6de0.jpg)
<details>
<summary>text_image</summary>
z
w3
3
u3
w2
2
u2
w1
1
u1
r
</details>
Figure 1622 Axisymmetric triangular element showing nodal degrees of freedom
<!-- source-page: 703 -->
Noting that
$$
\int_ {A} N _ {1} ^ {3} d A = \frac {2 A}{2 0} \quad \int_ {A} N _ {1} ^ {2} N _ {2} d A = \frac {2 A}{6 0} \tag {16.7.13}
$$
$$
\int_ {A} N _ {1} N _ {2} N _ {3} d A = \frac {2 A}{1 2 0} \quad \text { and so on }
$$
we obtain
$$
\underline {{m}} = \frac {\pi \rho A}{1 0} \left[ \begin{array}{c c c c c c} \frac {4}{3} r _ {1} + 2 \bar {r} & 0 & 2 \bar {r} - \frac {r _ {3}}{3} & 0 & 2 \bar {r} - \frac {r _ {2}}{3} & 0 \\ & \frac {4}{3} r _ {1} + 2 \bar {r} & 0 & 2 \bar {r} - \frac {r _ {3}}{3} & 0 & 2 \bar {r} - \frac {r _ {2}}{3} \\ & & \frac {4}{3} r _ {2} + 2 \bar {r} & 0 & 2 \bar {r} - \frac {r _ {1}}{3} & 0 \\ & & & \frac {4}{3} r _ {2} + 2 \bar {r} & 0 & 2 \bar {r} - \frac {r _ {1}}{3} \\ & & & & \frac {4}{3} r _ {3} + 2 \bar {r} & 0 \\ \text {Symmetry} & & & & & \frac {4}{3} r _ {3} + 2 \bar {r} \end{array} \right] \tag {16.7.14}
$$
where r ¼ r1 þ r2 þ r3 $\bar { r } = \frac { r _ { 1 } + r _ { 2 } + r _ { 3 } } { 3 }$
# Tetrahedral Solid Element
Finally, the tetrahedral solid element (considered in Chapter 11) consistent-mass matrix is obtained by substituting the shape function matrix Eq. (11.2.9) with shape functions defined in Eq. (11.2.10) into Eq. (16.2.19) and performing the integration to obtain
$$
\underline {{m}} = \frac {\rho V}{2 0} \left[ \begin{array}{c c c c c c c c c c c c} 2 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\ & 2 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ & & 2 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 \\ & & & 2 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\ & & & & & 2 & 0 & 0 & 1 & 0 & 0 & 1 \\ & & & & & & 2 & 0 & 0 & 1 & 0 & 0 \\ & & & & & & & 2 & 0 & 0 & 1 & 0 \\ & & & & & & & & 2 & 0 & 0 & 1 \\ & & & & & & & & & 2 & 0 & 0 \\ & & & & & & & & & & 2 & 0 \\ \text {Symmetry} & & & & & & & & & & & 2 \end{array} \right] \tag {16.7.15}
$$
<!-- source-page: 704 -->
# d 16.8 Time-Dependent Heat Transfer
In this section, we consider the time-dependent heat transfer problem in one dimension only. The basic differential equation for time-dependent heat transfer in one dimension was given previously by Eq. (13.1.7) with the boundary conditions given by Eqs. (13.1.10) and (13.1.11).
The finite element formulation of the equations can be obtained by minimization of the following functional:
$$
\begin{array}{l} \pi_ {h} = \frac {1}{2} \iiint_ {V} \left[ K _ {x x} \left(\frac {\partial T}{\partial x}\right) ^ {2} - 2 (Q - c \rho \dot {T}) T \right] d V \\ - \iint_ {S _ {2}} q ^ {*} T d S + \frac {1}{2} \iint_ {S _ {3}} h (T - T _ {\infty}) ^ {2} d s \tag {16.8.1} \\ \end{array}
$$
Equation (16.8.1) is similar to Eq. (13.4.10) with definitions given by Eq. (13.4.11) except that the Q term is now replaced by
$$
Q - c \rho \dot {T} \tag {16.8.2}
$$
where, again, c is the specific heat of the material, and the dot over the variable T denotes differentiation with respect to time. Again, Eq. (13.4.22) obtained in Section 13.4 for the conductivity or stiffness matrix and Eqs. (13.4.23)(13.4.25) for the force matrix terms are applicable here.
The term given by Eq. (16.8.2) yields an additional contribution to the basic element equations previously obtained for the time-independent problem as follows:
$$
\Omega_ {Q} = - \iiint_ {V} T (Q - c \rho \dot {T}) d V \tag {16.8.3}
$$
Again, the temperature function is given by
$$
\{T \} = [ N ] \{\hat {t} \} \tag {16.8.4}
$$
where ½N
is the shape function matrix given by Eq. (13.4.3) or Eqs. (16.2.3) for the simple one-dimensional element, and f^tg is the nodal temperature matrix. Substituting Eq. (16.8.4) into Eq. (16.8.3) and differentiating with respect to time where indicated yield
$$
\Omega_ {Q} = - \iint_ {V} \left([ N ] \{\hat {t} \} Q - c \rho [ N ] \{\hat {t} \} [ N ] \{\dot {\hat {t}} \}\right) d V \tag {16.8.5}
$$
where the fact that ½N
is a function only of the coordinate system has been taken into account. Equation (16.8.5) must be minimized with respect to the nodal temperatures as follows:
$$
\frac {\partial \Omega_ {Q}}{\partial \{\hat {t} \}} = - \iiint_ {V} [ N ] ^ {T} Q d V + \iiint_ {V} c \rho [ N ] ^ {T} [ N ] d V \{\dot {\hat {t}} \} \tag {16.8.6}
$$
<!-- source-page: 705 -->
where we have assumed that $\{ \dot { \hat { t } } \}$ remains constant during the differentiation with respect to f^tg. Equation (16.8.6) results in the additional time-dependent term added to Eq. (13.4.18). Hence, using previous definitions for the stiffness and force matrices, we obtain the element equations as
$$
\{\hat {f} \} = [ \hat {k} ] \{\hat {t} \} + [ \hat {m} ] \{\dot {\hat {t}} \} \tag {16.8.7}
$$
where now
$$
[ \hat {m} ] = \iint_ {V} c \rho [ N ] ^ {T} [ N ] d V \tag {16.8.8}
$$
For an element with constant cross-sectional area A, the differential volume is $d V = A d { \hat { x } }$ . Substituting the one-dimensional shape function matrix Eq. (13.4.3) into Eq. (16.8.8) yields
$$
[ \hat {m} ] = c \rho A \int_ {0} ^ {L} \left\{ \begin{array}{c} 1 - \frac {\hat {x}}{L} \\ \frac {\hat {x}}{L} \end{array} \right\} \left[ \begin{array}{c c} 1 - \frac {\hat {x}}{L} & \frac {\hat {x}}{L} \end{array} \right] d \hat {x}
$$
or
$$
[ \hat {m} ] = \frac {c \rho A L}{6} \left[ \begin{array}{l l} 2 & 1 \\ 1 & 2 \end{array} \right] \tag {16.8.9}
$$
Equation (16.8.9) is analogous to the consistent-mass matrix Eq. (16.2.23). The lumped-mass matrix for the heat conduction problem is then
$$
[ \hat {m} ] = \frac {c \rho A L}{2} \left[ \begin{array}{l l} 1 & 0 \\ 0 & 1 \end{array} \right] \tag {16.8.10}
$$
which is analogous to Eq. (16.2.12) for the one-dimensional stress element.
The time-dependent heat-transfer problem can now be solved in a manner analogous to that for the stress analysis problem. We present the numerical time integration scheme.
# Numerical Time Integration
The numerical time integration method described here is similar to Newmarks method used for structural dynamics analysis and can be used to solve time-dependent or transient heat-transfer problems.
We begin by assuming that two temperature states $\underline { { T } } _ { i }$ at time $t _ { i }$ and $\underline { { T } } _ { i + 1 }$ at time $t _ { i + 1 }$ are related by
$$
\underline {{T}} _ {i + 1} = \underline {{T}} _ {i} + [ (1 - \beta) \dot {\underline {{T}}} _ {i} + \beta \dot {\underline {{T}}} _ {i + 1} ] (\Delta t) \tag {16.8.11}
$$
Equation (16.8.11) is known as the generalized trapezoid rule. Much like Newmarks method for numerical time integration of the second-order equations of structural dynamics, Eq. (16.8.11) includes a parameter $\beta$ that is chosen by the user.
Next we express Eq. (16.8.7) in global form as
$$
\{F \} = [ K ] \{T \} + [ M ] \{\dot {T} \} \tag {16.8.12}
$$
<!-- source-page: 706 -->
We now write Eq. (16.8.12) for time $t _ { i }$ and then for time $t _ { i + 1 }$ . We then multiply the first of these two equations by $1 - \beta$ and the second by $\beta$ to obtain
$$
(1 - \beta) (\underline {{K}} \underline {{T}} _ {i} + \underline {{M}} \dot {\underline {{T}}} _ {i}) = (1 - \beta) \underline {{F}} _ {i} \tag {16.8.13a}
$$
$$
\beta (\underline {{K}} \underline {{T}} _ {i + 1} + \underline {{M}} \dot {\underline {{T}}} _ {i + 1}) = \beta \underline {{F}} _ {i + 1} \tag {16.8.13b}
$$
Next we add Eqs. (16.8.13a and b) together to obtain
$$
\underline {{M}} [ (1 - \beta) \dot {\underline {{T}}} _ {i} + \beta \dot {\underline {{T}}} _ {i + 1} ] + \underline {{K}} [ (1 - \beta) \underline {{T}} _ {i} + \beta \underline {{T}} _ {i + 1} ] = (1 - \beta) \underline {{F}} _ {i} + \beta \underline {{F}} _ {i + 1} \tag {16.8.14}
$$
Now, using Eq. (16.8.11), we can eliminate the time derivative terms from $\operatorname { E q }$ (16.8.14) to write
$$
\frac {\underline {{M}} (\underline {{T}} _ {i + 1} - \underline {{T}} _ {i})}{\Delta t} + \underline {{K}} [ (1 - \beta) \underline {{T}} _ {i} + \beta \underline {{T}} _ {i + 1} ] = (1 - \beta) \underline {{F}} _ {i} + \beta \underline {{F}} _ {i + 1} \tag {16.8.15}
$$
Rewriting Eq. (16.8.15) by grouping the $\underline { { T } } _ { i + 1 }$ terms on the left side, we have
$$
\left(\frac {1}{\Delta t} \underline {{M}} + \beta \underline {{K}}\right) \underline {{T}} _ {i + 1} = \left[ \frac {1}{\Delta t} \underline {{M}} - (1 - \beta) \underline {{K}} \right] \underline {{T}} _ {i} + (1 - \beta) \underline {{F}} _ {i} + \beta \underline {{F}} _ {i + 1} \tag {16.8.16}
$$
The time integration to solve for $\underline { T }$ begins as follows. Given a known initial temperature $\underline { { T } } _ { 0 }$ at time $t = 0$ and a time step $\Delta t ,$ we solve Eq. (16.8.16) for $\underline { { T } } _ { 1 }$ at $t = \Delta t$ . Then, using $\underline { { T } } _ { 1 }$ , we determine $\underline { { T } } _ { 2 }$ at $t = 2 ( \Delta t )$ , and so on. For a constant $\Delta t ,$ the left-side coefficient of $T _ { i + 1 }$ need be evaluated only one time (assuming $\underline { { M } }$ and $\underline { { K } }$ do not vary with time). The matrix Eq. (16.8.16) can then be solved in the usual manner, such as by Gauss elimination. For a one-dimensional heat-transfer analysis, element k is given by Eqs. (13.4.22) and (13.4.28), whereas f is given by Eqs. (13.4.26) and (13.4.29).
It has been shown that depending on the value of $\beta ,$ the time step $\Delta t$ may have an upper limit for the numerical analysis to be stable. If $\begin{array} { r } { \beta < \frac { 1 } { \bar { \gamma } } . } \end{array}$ , the largest $\Delta t$ for stability as shown in Reference [12] is
$$
\Delta t = \frac {2}{(1 - 2 \beta) \lambda_ {\max}} \tag {16.8.17}
$$
where $\lambda _ { \mathrm { m a x } }$ is the largest eigenvalue of
$$
(\underline {{K}} - \lambda \underline {{M}}) \underline {{T}} ^ {\prime} = 0 \tag {16.8.18}
$$
in which, as in Eq. (16.4.2), we have
$$
\underline {{T}} (t) = \underline {{T}} ^ {\prime} e ^ {i \lambda t} \tag {16.8.19}
$$
with $\underline { T } ^ { \prime }$ representing the natural modes. If $\begin{array} { r } { \beta \geqslant \frac { 1 } { 2 } , } \end{array}$ , the numerical analysis is unconditionally stable; that is, stability of solution (but not accuracy) is guaranteed for Dt greater than that given by Eq. (16.8.17), or as $\Delta t$ becomes indefinitely large. Various numerical integration methods result, depending on specific values of $\beta \colon$
$\beta = 0 :$ Forward difference, or Euler [3], which is said to be conditionally stable (that is, Dt must be no greater than that given by Eq. (16.8.17) to obtain a stable solution).
<!-- source-page: 707 -->
$\begin{array} { r } { \beta = \frac { 1 } { 2 } : } \end{array}$ Crank-Nicolson, or trapezoid, rule, which is unconditionally stable.
$\begin{array} { r } { \beta = \frac { 2 } { 3 } \colon } \end{array}$ Galerkin, which is unconditionally stable.
$\beta = 1$ : Backward difference, which is unconditionally stable.
If $\beta = 0 _ { ; }$ , the numerical integration method is called explicit; that is, we can solve for $\underline { { T } } _ { i + 1 }$ directly at time $\Delta t$ knowing only previous information at $t = \underline { { T } } _ { i } . \ : \mathrm { I f } \beta > 0$ , the method is called implicit. If a diagonal mass-type matrix M exists and $\beta = 0$ , the computational effort for each time step is small (see Example 16.4, where a lumped-mass matrix was used), but so must be Dt. The choice of $\beta > \frac { 1 } { 2 }$ is often used. However, if $\begin{array} { r } { \beta = \frac { 1 } { 2 } } \end{array}$ and sharp transients exist, the method generates spurious oscillations in the solution. Using $\begin{array} { r } { \beta > \frac { 1 } { 2 } , } \end{array}$ along with smaller $\Delta t$ [12], is probably better. Example 16.7 illustrates the solution of a one-dimensional time-dependent heat-transfer problem using the numerical time integration scheme [Eq. (16.8.16)].
# Example 16.7
A circular fin (Figure 1623) is made of pure copper with a thermal conductivity of $K _ { x x } = 4 0 0 ~ \mathrm { { W / ( m \cdot ^ { \circ } C ) } } , h = 1 5 0 ~ \mathrm { { W / ( m ^ { 2 } \cdot ^ { \circ } C ) } }$ , mass density $\rho = 8 9 0 0 ~ \mathrm { k g } / \mathrm { m } ^ { 3 }$ , and specific heat $c = 3 7 5 ~ \mathrm { J / ( k g \cdot ^ { \circ } C ) } ~ ( 1 ~ \mathrm { J } = 1 ~ \mathrm { W } \cdot \mathrm { s } )$ . The initial temperature of the fin is $2 5 ^ { \circ } \mathrm { C }$ . The fin length is 2 cm, and the diameter is 0.4 cm. The right tip of the fin is insulated. The base of the fin is then suddenly increased to a temperature of $8 5 ^ { \circ } \mathrm { C }$ and maintained at this temperature. Use the consistent form of the capacitance matrix, a time step of 0.1 s, and $\begin{array} { r } { \beta = \frac { 2 } { 3 } . } \end{array}$ . Use two elements of equal length. Determine the temperature distribution up to 3 s.
Using Eq. (13.4.22), the stiffness matrix is
$$
\begin{array}{l} \begin{array}{c c c c c} 1 & 2 & & 1 & 2 \\ 2 & 3 & & 2 & 3 \end{array} \\ \underline {{k}} ^ {(1)} = \underline {{k}} ^ {(2)} = \frac {A K _ {x x}}{L} \left[ \begin{array}{c c} 1 & - 1 \\ - 1 & 1 \end{array} \right] + \frac {h P L}{6} \left[ \begin{array}{c c} 2 & 1 \\ 1 & 2 \end{array} \right] \\ \underline {{k}} ^ {(1)} = \underline {{k}} ^ {(2)} = \frac {\pi (0 . 0 0 4) ^ {2} (4 0 0)}{4 (0 . 0 1)} \left[ \begin{array}{c c} 1 & - 1 \\ - 1 & 1 \end{array} \right] + \frac {1 5 0 (2 \pi) (0 . 0 0 2) (0 . 0 1)}{6} \left[ \begin{array}{c c} 2 & 1 \\ 1 & 2 \end{array} \right] \tag {16.8.20} \\ \end{array}
$$
![](images/page-707_d79c41082554a992e104ecda4b39ff26b8addbcb797794c0544d4b106b92d8ae.jpg)
<details>
<summary>text_image</summary>
T∞ = 25°C
1
2
3
Insulated tip
85°C
2 cm
</details>
Figure 1623 Rod subjected to time-dependent temperature
<!-- source-page: 708 -->
Assembling the element stiffness matrices, Eq. (16.8.20), we obtain the global stiffness matrix as
$$
\underline {{{K}}} = \left[ \begin{array}{c c c} 1 & 2 & 3 \\ 0. 5 0 8 9 4 & - 0. 4 9 9 5 1 & 0 \\ - 0. 4 9 9 5 1 & 1. 0 1 7 8 8 & - 0. 4 9 9 5 1 \\ 0 & - 0. 4 9 9 5 1 & 0. 5 0 8 9 4 \end{array} \right] \frac {\mathrm{W}}{\mathrm {^ {\circ} C}} \tag {16.8.21}
$$
Using Eq. (13.4.25), we obtain each element force matrix as
$$
\{f _ {h} ^ {(1)} \} = \{f _ {h} ^ {(2)} \} = \frac {h T _ {\infty} P L}{2} \left\{ \begin{array}{l} 1 \\ 1 \end{array} \right\} = \frac {(1 5 0) (2 5 ^ {\circ} \mathrm{C}) (2 \pi) (0 . 0 0 2) (0 . 0 1)}{2} \left\{ \begin{array}{l} 1 \\ 1 \end{array} \right\}
$$
$$
\underline {{f}} _ {h} ^ {(1)} = \underline {{f}} _ {h} ^ {(2)} = \left\{ \begin{array}{l} 0. 2 3 5 6 1 \\ 0. 2 3 5 6 1 \end{array} \right\} \tag {16.8.22}
$$
Using Eq. (16.8.22), we find that the assembled global force matrix is
$$
\{F \} = \left\{ \begin{array}{l} 0. 2 3 5 6 1 \\ 0. 4 7 1 2 2 \\ 0. 2 3 5 6 1 \end{array} \right\} \mathrm{W} \tag {16.8.23}
$$
Next using Eq. (16.8.9), we obtain each element mass (capacitance) matrix as
$$
[ m ] = \frac {c \rho A L}{6} \left[ \begin{array}{c c} 2 & 1 \\ 1 & 2 \end{array} \right]
$$
$$
\underline {{m}} ^ {(1)} = \underline {{m}} ^ {(2)} = \frac {(3 7 5) (8 9 0 0) \frac {\pi (0 . 0 0 4) ^ {2}}{4} (0 . 0 1)}{6} \left[ \begin{array}{c c} 2 & 1 \\ 1 & 2 \end{array} \right]
$$
$$
= 0. 0 6 9 9 0 \left[ \begin{array}{l l} 2 & 1 \\ 1 & 2 \end{array} \right] \mathbf {W} \cdot \mathrm{s} / ^ {\circ} \mathbf {C} \tag {16.8.24}
$$
Using Eq. (16.8.24), the assembled capacitance matrix is
$$
\underline {{M}} = \left[ \begin{array}{l l l} 1 & 2 & 3 \\ 0. 1 3 9 8 0 & 0. 0 6 9 9 0 & 0 \\ 0. 0 6 9 9 0 & 0. 2 7 9 6 0 & 0. 0 6 9 9 0 \\ 0 & 0. 0 6 9 9 0 & 0. 1 3 9 8 0 \end{array} \right] \frac {\mathrm{W} \cdot \mathrm{s}}{\circ \mathrm{C}} \tag {16.8.25}
$$
Using Eq. (16.8.16) and Eqs. (16.8.21) and (16.8.25), we obtain
$$
\left(\frac {1}{\Delta t} \underline {{M}} + \beta \underline {{K}}\right) = \left[ \begin{array}{l l l} 1. 7 3 7 4 & 0. 3 6 6 0 3 & 0 \\ 0. 3 6 6 0 3 & 3. 4 7 4 7 & 0. 3 6 6 0 3 \\ 0 & 0. 3 6 6 0 3 & 1. 7 3 7 4 \end{array} \right] \frac {\mathrm{W}}{\circ \mathrm{C}} \tag {16.8.26}
$$
<!-- source-page: 709 -->
Table 164 Nodal temperatures at various times for Example 16.7
<table><tr><td rowspan="2">Time (s)</td><td colspan="3">Temperature of Node Numbers (°C)</td></tr><tr><td>1</td><td>2</td><td>3</td></tr><tr><td>0.1</td><td>85</td><td>18.534</td><td>26.371</td></tr><tr><td>0.2</td><td>85</td><td>29.732</td><td>21.752</td></tr><tr><td>0.3</td><td>85</td><td>36.404</td><td>22.662</td></tr><tr><td>0.4</td><td>85</td><td>41.032</td><td>25.655</td></tr><tr><td>0.5</td><td>85</td><td>44.665</td><td>29.312</td></tr><tr><td>0.6</td><td>85</td><td>47.749</td><td>33.059</td></tr><tr><td>0.7</td><td>85</td><td>50.482</td><td>36.669</td></tr><tr><td>0.8</td><td>85</td><td>52.956</td><td>40.062</td></tr><tr><td>0.9</td><td>85</td><td>55.218</td><td>43.218</td></tr><tr><td>1.0</td><td>85</td><td>57.296</td><td>46.139</td></tr><tr><td>1.1</td><td>85</td><td>59.208</td><td>48.837</td></tr><tr><td>1.2</td><td>85</td><td>60.969</td><td>51.327</td></tr><tr><td>1.3</td><td>85</td><td>62.593</td><td>53.623</td></tr><tr><td>1.4</td><td>85</td><td>64.089</td><td>55.741</td></tr><tr><td>1.5</td><td>85</td><td>65.469</td><td>57.693</td></tr><tr><td>1.6</td><td>85</td><td>66.742</td><td>59.493</td></tr><tr><td>1.7</td><td>85</td><td>67.915</td><td>61.152</td></tr><tr><td>1.8</td><td>85</td><td>68.996</td><td>62.683</td></tr><tr><td>1.9</td><td>85</td><td>69.993</td><td>64.094</td></tr><tr><td>2.0</td><td>85</td><td>70.912</td><td>65.395</td></tr><tr><td>2.1</td><td>85</td><td>71.760</td><td>66.594</td></tr><tr><td>2.2</td><td>85</td><td>72.542</td><td>67.700</td></tr><tr><td>2.3</td><td>85</td><td>73.262</td><td>68.720</td></tr><tr><td>2.4</td><td>85</td><td>73.926</td><td>69.660</td></tr><tr><td>2.5</td><td>85</td><td>74.539</td><td>70.527</td></tr><tr><td>2.6</td><td>85</td><td>75.104</td><td>71.326</td></tr><tr><td>2.7</td><td>85</td><td>75.624</td><td>72.063</td></tr><tr><td>2.8</td><td>85</td><td>76.104</td><td>72.742</td></tr><tr><td>2.9</td><td>85</td><td>76.547</td><td>73.368</td></tr><tr><td>3.0</td><td>85</td><td>76.955</td><td>73.946</td></tr></table>$$
\text { and } \quad \left[ \frac {1}{\Delta t} \underline {{M}} - (1 - \beta) \underline {{K}} \right] = \left[ \begin{array}{l l l} 1. 2 2 8 0 & 0. 8 6 5 5 & 0 \\ 0. 8 6 5 5 & 2. 4 5 7 & 0. 8 6 5 5 \\ 0 & 0. 8 6 5 5 & 1. 2 2 8 0 \end{array} \right] \frac {\mathrm{W}}{\mathrm {^ {\circ} C}} \tag {16.8.27}
$$
where $\begin{array} { r } { \beta = \frac { 2 } { 3 } } \end{array}$ and $\Delta t = 0 .$ 1 s have been used to obtain Eqs. (16.8.26) and (16.8.27). For the first time step, $t = 0 . 1 \ { \mathrm { s } } ,$ , we then use Eqs. (16.8.23), (16.8.27), and (16.8.26) in
<!-- source-page: 710 -->
Eq. (16.8.16) to obtain
$$
\begin{array}{l} \left[ \begin{array}{c c c} 1. 7 3 7 4 & 0. 3 6 6 0 3 & 0 \\ 0. 3 6 6 0 3 & 3. 4 7 4 7 & 0. 3 6 6 0 3 \\ 0 & 0. 3 6 6 0 3 & 1. 7 3 7 4 \end{array} \right] \left\{ \begin{array}{c} 8 5 ^ {\circ} \mathrm{C} \\ t _ {2} \\ t _ {3} \end{array} \right\} \\ = \left[ \begin{array}{l l l} 1. 2 2 8 0 & 0. 8 6 5 5 & 0 \\ 0. 8 6 5 5 & 2. 4 5 7 & 0. 8 6 5 5 \\ 0 & 0. 8 6 5 5 & 1. 2 2 8 0 \end{array} \right] \left\{ \begin{array}{l} 2 5 ^ {\circ} \mathrm{C} \\ 2 5 ^ {\circ} \mathrm{C} \\ 2 5 ^ {\circ} \mathrm{C} \end{array} \right\} + \left\{ \begin{array}{l} 0. 2 3 5 6 1 \\ 0. 4 7 1 2 2 \\ 0. 2 3 5 6 1 \end{array} \right\} \tag {16.8.28} \\ \end{array}
$$
In Eq. (16.8.28), we should note that because $\underline { { F } } _ { i } = \underline { { F } } _ { i + 1 }$ for all time, the sum of the terms is $( 1 - \beta ) \underline { { F } } _ { i } + \beta \underline { { F } } _ { i + 1 } = \underline { { F } } _ { i }$ for all time. This is the column matrix on the right side of Eq. (16.8.28). We now solve Eq. (16.8.28) in the usual manner by partitioning the second and third equations of Eq. (16.8.28) from the first equation and solving the second and third equations simultaneously for $t _ { 2 }$ and $t _ { 3 }$ . The results are
$$
t _ {2} = 1 8. 5 3 4 ^ {\circ} \mathrm{C} \quad t _ {3} = 2 6. 3 7 1 ^ {\circ} \mathrm{C}
$$
At time $t = 0 . 2 ~ \mathrm { s } ,$ Eq. (16.8.28) becomes
$$
\begin{array}{l} \left[ \begin{array}{c c c} 1. 7 3 7 4 & 0. 3 6 6 0 3 & 0 \\ 0. 3 6 6 0 3 & 3. 4 7 4 7 & 0. 3 6 6 0 3 \\ 0 & 0. 3 6 6 0 3 & 1. 7 3 7 4 \end{array} \right] \left\{ \begin{array}{c} 8 5 ^ {\circ} \mathrm{C} \\ t _ {2} \\ t _ {3} \end{array} \right\} \\ = \left[ \begin{array}{l l l} 1. 2 2 8 0 & 0. 8 6 5 5 & 0 \\ 0. 8 6 5 5 & 2. 4 5 7 & 0. 8 6 5 5 \\ 0 & 0. 8 6 5 5 & 1. 2 2 8 0 \end{array} \right] \left\{ \begin{array}{c} 8 5 ^ {\circ} \mathrm{C} \\ 1 8. 5 3 4 ^ {\circ} \mathrm{C} \\ 2 6. 3 7 1 ^ {\circ} \mathrm{C} \end{array} \right\} + \left\{ \begin{array}{l} 0. 2 3 5 6 1 \\ 0. 4 7 1 2 2 \\ 0. 2 3 5 6 1 \end{array} \right\} \tag {16.8.29} \\ \end{array}
$$
![](images/page-710_0604522f0ea37e64bb8c25d575501c902124a5597fa6e6ed069bc58683c9739e.jpg)
<details>
<summary>line</summary>
| Time, s | Node 2 Temperature, °C | Node 3 Temperature, °C |
| ------- | ---------------------- | ---------------------- |
| 0.0 | 25.0 | 25.0 |
| 0.2 | 20.0 | 20.0 |
| 0.4 | 30.0 | 25.0 |
| 0.6 | 40.0 | 35.0 |
| 0.8 | 50.0 | 45.0 |
| 1.0 | 55.0 | 50.0 |
| 1.2 | 60.0 | 55.0 |
| 1.4 | 65.0 | 60.0 |
| 1.6 | 68.0 | 63.0 |
| 1.8 | 70.0 | 65.0 |
| 2.0 | 72.0 | 67.0 |
| 2.2 | 73.0 | 68.0 |
| 2.4 | 74.0 | 69.0 |
| 2.6 | 75.0 | 70.0 |
| 2.8 | 76.0 | 71.0 |
| 3.0 | 77.0 | 72.0 |
</details>
Figure 1624 Temperature as a function of time for nodes 2 and 3 of Example 16.7