MultiPhysicsVault/.raw/FiniteElementProcedures/FiniteElementProcedures_085.md

9.20. An 11-story tapered tower is subjected to an air blast as shown. Use a computer program to solve for the response of the tower.

11 story tapered tower

bar

Component	Value
Girder properties:	E = 2.07 × 10^11 Pa
	v = 0.3
	A = 0.01 m^2
	A_s = 0.009 m^2
	I = 8.33 × 10^-5 m^4
	ρ = 7800 kg/m^3

Applied load (blast):

line

Time (msec)	Force per unit length (N/m)
0	2000
50	1000
100	500
150	200
200	0

9.5 SOLUTION OF NONLINEAR EQUATIONS IN DYNAMIC ANALYSIS

The solution of the nonlinear dynamic response of a finite element system is, in essence, obtained using the procedures already discussed: the incremental formulations presented in Chapter 6, the iterative solution procedures discussed in Section 8.4, and the time integration algorithms presented in this chapter. Hence, the major basic procedures used in a nonlinear dynamic response solution have already been presented, and we only need to briefly summarize in the following how these procedures are employed together in a nonlinear dynamic analysis.

9.5.1 Explicit Integration

The most common explicit time integration operator used in nonlinear dynamic analysis is probably the central difference operator. As in linear analysis (see Section 9.2), the equilibrium of the finite element system is considered at time t in order to calculate the displacements at time t + \Delta t . Neglecting the effect of a damping matrix, we operate for each discrete time step solution on the equations,

\mathbf {M} ^ {\prime} \ddot {\mathbf {U}} = ^ {\prime} \mathbf {R} - ^ {\prime} \mathbf {F} \tag {9.103}

where the nodal point force vector {}^{t}F is evaluated as discussed in Section 6.3. The solution for the nodal point displacements at time t + \Delta t is obtained using the central difference approximation for the accelerations [given in (9.3)] to substitute for {}^{t}\ddot{U} in (9.103). Thus, as in linear analysis, if we know {}^{t-\Delta t}U and {}^{t}U , the relations in (9.3) and (9.103) are employed to calculate {}^{t+\Delta t}U . The solution therefore simply corresponds to a forward marching in time; the main advantage of the method is that with M a diagonal matrix the solution of {}^{t+\Delta t}U does not involve a triangular factorization of a coefficient matrix.

The shortcoming in the use of the central difference method lies in the severe time step restriction: for stability, the time step size \Delta t must be smaller than a critical time step \Delta t_{cr} , which is equal to T_{n}/\pi , where T_{n} is the smallest period in the finite element system. This time step restriction was derived considering a linear system (see Example 9.12), but the result is also applicable to nonlinear analysis, since for each time step the nonlinear response calculation may be thought of—in an approximate way—as a linear analysis. However, whereas in a linear analysis the stiffness properties remain constant, in a nonlinear analysis these properties change during the response calculations. These changes in the material and/or geometric conditions enter into the evaluation of the force vector {}^{\prime}\mathbf{F} as discussed in Chapter 6. Since therefore the value of T_{n} is not constant during the response calculation, the time step \Delta t needs to be decreased if the system stiffens, and this time step adjustment must be performed in a conservative manner so that the condition \Delta t \leq T_{n}/\pi is satisfied with certainty at all times.

To emphasize this point, consider an analysis in which the time step is always smaller than the critical time step except for a few successive solution steps, and for these solution steps the time step \Delta t is just slightly larger than the critical time step. In such a case, the analysis results may not show an “obvious” solution instability but instead a significant solution error is accumulated over the solution steps for which the time step size was larger than the critical value for stability. The situation is quite different from what is observed in linear analysis, where the solution quickly “blows up” if the time step is larger than the critical time step size for stability. This phenomenon is demonstrated somewhat in the response predicted for the simple one degree of freedom spring-mass system shown in Fig. 9.8. In the solution the time step \Delta t is slightly larger than the critical time step for stability in the stiff region of the spring. Since the time step corresponds to a stable time step for small spring displacements, the response calculations are partly stable and partly unstable. The calculated response is shown in Fig 9.8, and it is observed that although the predicted displacements are grossly in error, the solution does not blow up. Hence, if this single degree of freedom system corresponded to a higher frequency in a large finite element model, a significant error accumulation could take place without an obvious blow-up of the solution.

The proper choice of the time step \Delta t is therefore a most important factor. Guidelines for the choice of \Delta t have been given in Section 9.4.4.

^{0}x=0.0
^{0}\dot{x}=0.955
x_{T} = 0.0025

Figure 9.8 Response of bilinear elastic system as predicted using the central difference method; \Delta t_{cr} = 0.001061027 ; the accurate response with displacement \leqslant 0.1 was calculated with \Delta t = 0.000106103 ; the “unstable” response was calculated with \Delta t = 0.00106103 .

line

Time	Displacement
0.0	0.0
0.1	0.0
0.2	0.0
0.3	0.0
0.4	0.0
0.5	0.0
0.6	0.0

Figure 9.8 (continued)

9.5.2 Implicit Integration

All the implicit time integration schemes discussed previously for linear dynamic analysis can also be employed in nonlinear dynamic response calculations. Although unstable in certain analyses, see K. J. Bathe [F], a very common technique used is the trapezoidal rule and we use this method to demonstrate the basic additional considerations involved in a nonlinear analysis.

As in linear analysis, using implicit time integration, we consider the equilibrium of the system at time t + \Delta t . This requires in nonlinear analysis that an iteration be performed. In practice, a full Newton-Raphson is usually best (see Section 8.4), but for illustration we consider here the modified Newton-Raphson iteration (and also neglect damping effects):

\mathbf {M} ^ {t + \Delta t} \ddot {\mathbf {U}} ^ {(k)} + ^ {t} \mathbf {K} \Delta \mathbf {U} ^ {(k)} = ^ {t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} ^ {(k - 1)} \tag {9.104}

{ } ^ { t + \Delta t } \mathbf { U } ^ { ( k ) } = { } ^ { t + \Delta t } \mathbf { U } ^ { ( k - 1 ) } + \Delta \mathbf { U } ^ { ( k ) } \tag {9.105}

Using the trapezoidal rule of time integration (Newmark's method with \delta = \frac{1}{2} and \alpha = \frac{1}{4} ), the following assumptions are employed:

{ } ^ { t + \Delta t } \mathbf { U } = { } ^ { t } \mathbf { U } + \frac { \Delta t } { 2 } ( { } ^ { t } \dot { \mathbf { U } } + { } ^ { t + \Delta t } \dot { \mathbf { U } } ) \tag {9.106}

^ {t + \Delta t} \dot {\mathbf {U}} = ^ {t} \dot {\mathbf {U}} + \frac {\Delta t}{2} (^ {t} \ddot {\mathbf {U}} + ^ {t + \Delta t} \ddot {\mathbf {U}}) \tag {9.107}

Using the relations in (9.105) to (9.107), we thus obtain

{ } ^ { t + \Delta t } \ddot { \mathbf { U } } ^ { ( k ) } = \frac { 4 } { \Delta t ^ { 2 } } ( { } ^ { t + \Delta t } \mathbf { U } ^ { ( k - 1 ) } - { } ^ { t } \mathbf { U } + \Delta \mathbf { U } ^ { ( k ) } ) - \frac { 4 } { \Delta t } { } ^ { t } \dot { \mathbf { U } } - { } ^ { t } \ddot { \mathbf { U } } \tag {9.108}

and substituting into (9.104), we have

{ } ^ { t } \hat { \mathbf { K } } \Delta \mathbf { U } ^ { ( k ) } = { } ^ { t + \Delta t } \mathbf { R } - { } ^ { t + \Delta t } \mathbf { F } ^ { ( k - 1 ) } - \mathbf { M } \left( \frac { 4 } { \Delta t ^ { 2 } } ( { } ^ { t + \Delta t } \mathbf { U } ^ { ( k - 1 ) } - { } ^ { t } \mathbf { U } ) - \frac { 4 } { \Delta t } { } ^ { t } \dot { \mathbf { U } } - { } ^ { t } \ddot { \mathbf { U } } \right) \tag {9.109}

where

^ \prime \hat {\mathbf {K}} = ^ {\prime} \mathbf {K} + \frac {4}{\Delta t ^ {2}} \mathbf {M} \tag {9.110}

We now notice that the iterative equations in dynamic nonlinear analysis using implicit time integration are of the same form as the equations that we considered in static nonlinear analysis, except that both the coefficient matrix and the nodal point force vector contain contributions from the inertia of the system. We can therefore directly conclude that all iterative solution strategies discussed in Section 8.4 for static analysis are also directly applicable to the solution of (9.109). However, since the inertia of the system renders its dynamic response, in general, “more smooth” than its static response, convergence of the iteration can, in general, be expected to be more rapid than in static analysis, and the convergence behavior can be improved by decreasing \Delta t . The numerical reason for the better convergence characteristics in a dynamic analysis as \Delta t decreases lies in the contribution of the mass matrix to the coefficient matrix. This contribution increases and ultimately becomes dominant as the time step decreases (see Section 8.4.1).

It is interesting to note that in the first solutions of nonlinear dynamic finite element response, equilibrium iterations were not performed in the step-by-step incremental analysis; i.e., the relation in (9.109) was simply solved for k = 1 and the incremental displacement \Delta\mathbf{U}^{(1)} was accepted as an accurate approximation to the actual displacement increment from time t to time t + \Delta t . However, it was then recognized that the iteration can actually be of utmost importance (see K. J. Bathe and E. L. Wilson [B]) since any error admitted in the incremental solution at a particular time directly affects in a path-dependent manner the solution at any subsequent time. Indeed, because any nonlinear dynamic response is highly path-dependent, the analysis of a nonlinear dynamic problem requires iteration at each time step more stringently than does a static analysis.

A simple demonstration of this observation is given in Fig. 9.9. This figure shows the results obtained in the analysis of a simple pendulum that was idealized as a truss element with a concentrated mass at its free end. The pendulum was released from a horizontal position, and the response was calculated for about one period of oscillation. In the analysis the convergence tolerances already discussed in Section 8.4.4 were used, but including the

line

| t (sec) | Angle θ (degrees) - tight tolerances | Angle θ (degrees) - loose tolerances | | ------- | ------------------------------------- | -------------------------------------- | | 0 | 90 | 90 | | 1 | 0 | 0 | | 2 | -60 | -60 | | 3 | 0 | 0 | | 4 | 60 | 60 |

Figure 9.9 Analysis of simple pendulum using trapezoidal rule, RNORM = mg

effect of inertia, i.e., convergence is reached when the following conditions are satisfied:

\frac {\left\| ^ {t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} ^ {(i - 1)} - \mathbf {M} ^ {t + \Delta t} \ddot {\mathbf {U}} ^ {(i - 1)} \right\| _ {2}}{\text { RNORM }} \leq \text { RTOL } \tag {9.111}

and \frac{\Delta\mathbf{U}^{(i)T}(t^{+\Delta t}\mathbf{R}-t^{+\Delta t}\mathbf{F}^{(i-1)}-\mathbf{M}^{t+\Delta t}\ddot{\mathbf{U}}^{(i-1)})}{\Delta\mathbf{U}^{(1)T}(t^{+\Delta t}\mathbf{R}-t^{\prime}\mathbf{F}-\mathbf{M}^{\prime}\ddot{\mathbf{U}})}\leq\mathrm{ETOL} (9.112)

in which RTOL is a force tolerance and ETOL is an energy tolerance. Figure 9.9 demonstrates the importance of iterating and doing so with a sufficiently tight convergence tolerance. In this analysis energy is lost if the convergence tolerance is not tight enough, but depending on the problem being considered, the predicted response may also blow up. We used in this example the Newmark trapezoidal rule, and obtained good results. However, it is known that this method can actually become unstable in nonlinear dynamic solutions (even though very tight convergence tolerances for the iterations are used). For such analyses, the Bathe method is effective for the reasons discussed in Section 9.4.3, because no parameters need to be set and only a reasonable time step size needs to be chosen, like in any direct time integration solution, see K.J. Bathe [F] and Z. Kazanci and K.J. Bathe [A].

9.5.3 Solution Using Mode Superposition

In considering linear analysis, we discussed in Section 9.3 that the essence of mode superposition is a transformation from the element nodal point degrees of freedom to the generalized degrees of freedom of the vibration mode shapes. Since the dynamic equilibrium equations in the basis of the mode shape vectors decouple (assuming proportional damping), mode superposition analysis can be very effective in linear analysis if only some vibration modes are excited by the loading. The same basic principles are also applicable in nonlinear analysis; however, in this case the vibration mode shapes and frequencies change, and to transform the coefficient matrix in (9.109) into diagonal form the free-vibration mode shapes of the system at time t need to be used in the transformation. The calculation of the vibration mode shapes and frequencies at time t, when these quantities have been calculated at a previous time, could be achieved economically using the subspace iteration method (see Section 11.6). However, the complete mode superposition analysis of nonlinear dynamic response is generally effective only when the solution can be obtained without updating the stiffness matrix too frequently. In this case, the governing finite element equilibrium equations for the solution of the response at time t + \Delta t are

\mathbf {M} ^ {t + \Delta t} \ddot {\mathbf {U}} ^ {(k)} + ^ {r} \mathbf {K} \Delta \mathbf {U} ^ {(k)} = ^ {t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} ^ {(k - 1)} \quad k = 1, 2, \dots \tag {9.113}

where ^{7}K is the stiffness matrix corresponding to the configuration at some previous time \tau . In the mode superposition analysis we now use

{ } ^ { t + \Delta t } \mathbf { U } = \sum _ { i = r } ^ { s } \boldsymbol { \phi } _ { i } { } ^ { t + \Delta t } x _ { i } \tag {9.114}

where t+\Delta t x_{i} is the ith generalized modal displacement at time t + \Delta t and

{ } ^ { \tau } \mathbf { K } \phi _ { i } = \omega _ { i } ^ { 2 } \mathbf { M } \phi _ { i } ; \quad i = r , \dots , s \tag {9.115}

that is, \omega_{i} , \phi_{i} are free-vibration frequencies (radians/second) and mode shape vectors of the system at time \tau . Using (9.114) in the usual way, the equations in (9.113) are transformed to

{ } ^ { t + \Delta t } \ddot { \mathbf { X } } ^ { ( k ) } + \mathbf { \Omega } ^ { 2 } \Delta \mathbf { X } ^ { ( k ) } = \mathbf { \Phi } ^ { T } ( { } ^ { t + \Delta t } \mathbf { R } - { } ^ { t + \Delta t } \mathbf { F } ^ { ( k - 1 ) } ) \quad k = 1 , 2 , \dots \tag {9.116}

where

\boldsymbol {\Omega} ^ {2} = \left[ \begin{array}{c c c c} \omega_ {r} ^ {2} & & & \\ & \cdot & & \\ & & \cdot & \omega_ {s} ^ {2} \end{array} \right]; \quad \boldsymbol {\Phi} = \left[ \begin{array}{c c c c} \boldsymbol {\Phi} _ {r}, & \dots , & \boldsymbol {\Phi} _ {s} \end{array} \right]; \quad \quad^ {t + \Delta t} \mathbf {X} = \left[ \begin{array}{c} ^ {t + \Delta t} x _ {r} \\ \vdots \\ ^ {t + \Delta t} x _ {s} \end{array} \right] \tag {9.117}

The relations in (9.116) are the equilibrium equations at time t + \Delta t in the generalized modal displacements of time \tau ; the corresponding mass matrix is an identity matrix, the stiffness matrix is \Omega^{2} , the external load vector is \Phi^{T t+\Delta t}R , and the force vector corresponding to the element stresses at the end of iteration (k - 1) is \Phi^{T t+\Delta t}F^{(k-1)} . The solution of (9.116) can be obtained using, for example, the trapezoidal rule of time integration (see Section 9.5.2).

In general, the use of mode superposition in nonlinear dynamic analysis can be effective if only a relatively few mode shapes need to be considered in the analysis. Such conditions may be encountered, for example, in the analysis of earthquake response and vibration excitation, and it is in these areas that the technique has been employed.

9.5.4 Exercises

9.21. Consider the simple pendulum idealization shown. Use a finite element program to solve for the response of the system (see Fig. 9.9), for 100 cycles.

text_image

g = 980 cm/sec² Initial conditions: 0θ = 90° 0θ̂ = 0 Length = 304.43 cm Tip mass = 10 kg EA = 10¹⁰ (kg·cm)/sec²

One truss element with tip concentrated mass is employed

9.22. Consider the cantilever beam shown. The beam is initially at rest when the tip load P is suddenly applied. Use a finite element program to solve for the dynamic response of the beam, allowing for the large displacement effects. Use the trapezoidal rule, the central difference method and, if available, also mode superposition to solve for the response.

text_image

10 in P 0.1 in P 0.5 lb Time

Young's modulus E = 1.2 \times 10^{6} lb/in ^{2} Poisson's ratio \nu = 0.2 Mass density \rho = 1.0 \times 10^{-4} (lb·sec ^{2} )/in ^{4} Width = 1 in

9.23. Use a computer program to solve for the dynamic buckling load of the arch considered in Fig. 6.23.
9.24. Use a computer program to analyze the pipe whip problem described in the figure. You can perform a direct integration solution or a mode superposition solution. (These types of problems are of importance in the analysis of postulated accident conditions; see, for example, S. M. Ma and K. J. Bathe [A].)

text_image

D₀ = 30.0 in a = 3.0 in D₁ = 27.75 in b = 21.0 in L = 360.0 in d = 5.75 in Diameter d a b D₀ D₁ D₂ Restraint

line

t	P(t) (lb)
0	6.57 × 10⁵

line

ε	σ_y
0	0
ε	σ_y

Pipe material:

E = 2. 6 9 8 \times 1 0 ^ {7} \text {   psi   }

v = 0. 3

\sigma_ {Y} = 2. 9 1 4 \times 1 0 ^ {4} \text { psi }

\rho = 7. 1 8 \times 1 0 ^ {- 4} \text {   slug / in } ^ {3}

line

| Parameter | Value | | --------------- | --------------- | | E | 2.99 × 10⁷ psi | | σ_y | 3.80 × 10⁴ psi |

9.6 SOLUTION OF NONSTRUCTURAL PROBLEMS; HEAT TRANSFER AND FLUID FLOWS

Although we considered in the previous sections the solution of the dynamic response of structures and solids, it should be recognized that many of the basic concepts discussed are also directly applicable to the analysis of other types of problems. Namely, in the solution of a nonstructural problem, the choice lies again between the use of explicit or implicit time integration, mode superposition analysis needs to be considered, and it may also be advantageous to use different time integration schemes for different domains of the complete element assemblage (see Section 9.2.5). The stability and accuracy properties of the time integration schemes used are analyzed in basically the same way as in structural analysis, and the important basic observations made concerning nonlinear structural analysis are also applicable to the analysis of nonlinear nonstructural problems.

9.6.1 The \alpha -Method of Time Integration

The nonstructural problems that we have in mind are heat transfer, field problems, and fluid flow (see Chapter 7). The major difference in the time integration of the governing equations of these problems, when compared to structural analysis, is that we now deal only with first derivatives in time. Therefore, time integration operators different from those discussed in the previous sections are employed.

Based on the discussion in Section 9.4 we can present a time integration scheme used in the analysis of heat transfer and fluid flows by considering a typical one degree of freedom equilibrium equation,

\dot {\eta} + \lambda \eta = r; \quad \eta | _ {t = 0} = ^ {0} \eta \tag {9.118}

where, for example, considering a heat transfer problem, \eta is the unknown temperature, \lambda is the diffusivity, and r is the heat input to the system. The \alpha -method of time integration can be employed effectively for the solution of (9.118) and is given by the following assumptions:

{ } ^ { t + \alpha \Delta t } \dot { \eta } = ( { } ^ { t + \Delta t } \eta - { } ^ { t } \eta ) / \Delta t \tag {9.119}

{ } ^ { t + \alpha \Delta t } \eta = ( 1 - \alpha ) ^ { t } \eta + \alpha ^ { t + \Delta t } \eta

where \alpha is a constant that is chosen to yield optimum stability and accuracy properties. To solve for t^{+\Delta t}\eta we proceed as described in Section 9.2. Namely, if \eta is known, we can use (9.118) at time t + \alpha \Delta t and (9.119) to solve for t^{+\Delta t}\eta , and so on. This \alpha -method was already used in Section 6.6.3 in the solution of the inelastic response of finite element systems.

The properties of the integration procedure depend on the value of \alpha that is employed. The following procedures are in frequent use (see, for example, L. Collatz [A]).

\alpha = 0 , explicit Euler forward method, stable provided \Delta t \leq 2 / \lambda , first-order accurate in \Delta t

\alpha = \frac{1}{2} , implicit trapezoidal rule, unconditionally stable, second-order accurate in \Delta t

\alpha = 1 , implicit Euler backward method, unconditionally stable, first-order accurate in \Delta t .

To evaluate these stability properties we proceed as in Section 9.4.2. Now we use (9.118), with r = 0 and \lambda constant, at time t + \alpha \Delta t , substitute from (9.119), and solve for the variable n at time t + \Delta t in terms of all known quantities.

{ } ^ { t + \Delta t } \eta = \frac { 1 - ( 1 - \alpha ) \lambda \Delta t } { 1 + \alpha \lambda \Delta t } { } ^ { t } \eta \tag {9.120}

Therefore, for stability we need

\left| \frac {1 - (1 - \alpha) \lambda \Delta t}{1 + \alpha \lambda \Delta t} \right| \leq 1 \tag {9.121}

which shows that the \alpha -method is unconditionally stable if \alpha \geq \frac{1}{2} . In the case \alpha < \frac{1}{2} , the method is only stable provided

\Delta t \leq \frac {2}{(1 - 2 \alpha) \lambda} \tag {9.122}

These results have been used in the summary given above for the cases \alpha = 0, \frac{1}{2} , and 1.

To evaluate the accuracy properties we proceed in a similar manner as in Section 9.4.3 but now consider the initial value problem

\dot {\eta} + \lambda \eta = 0; \quad {} ^ {0} \eta = 1 \tag {9.123}

Assume that we perform the numerical solution of (9.123) for a time period 1/\lambda , with time step magnitude 1/n\lambda , where n is the number of time steps used for the time period 1/\lambda . Then we can define as our error measure the percentage in absolute difference between the numerical and exact solutions at time 1/\lambda . Figure 9.10 shows this error measure for the Euler forward and backward methods and the trapezoidal rule.

The information in Fig. 9.10 is useful for a direct integration solution if the largest value of \lambda in the finite element mesh, for which the response is to be accurately integrated, is known. Namely, the relation in (9.118) can be considered as the governing differential equation in time for the mode shape corresponding to the value \lambda (see the discussion in Sections 9.3.2 and 9.4.4).

line

Time	Analytical	Numerical
0.0	1.0	1.0
1/λ	0.0	0.0

line

Number of steps used per 1/λ time period	α = 0	α = 0.5	α = 1
3	20.0	1.0	14.0
4	14.0	1.0	11.0
5	12.0	1.0	9.0
6	10.0	1.0	8.0
7	8.0	1.0	7.0
8	7.0	1.0	6.0
9	6.0	1.0	5.0
10	5.0	1.0	4.0

Figure 9.10 Error in numerical solution of (9.123) using the \alpha method

However, in actual analysis, we need to select the finite element discretization, the time integration scheme, and the time step \Delta t . The following considerations are then useful.

Consider a one-dimensional heat transfer condition in a continuous medium, as shown in Fig. 9.11. The uniform initial temperature is \theta_{i} , when suddenly the free surface at x = 0 is subjected to a temperature \theta_{0} . The governing differential equation of equilibrium of this problem was derived in Example 3.16.

The exact solution of the mathematical model gives the temperature distributions shown schematically in Fig. 9.11(b). This figure also shows the penetration depth \gamma

text_image

θ(x = 0, t = 0) = θi θ(x = 0, t = 0⁺) = θ0 x ∞ ∞ Initial temperature θi θ(x, 0) = θi, x ≥ 0

(a) Problem considered

line

x	Time t1	Time t2	Time t3
0	θ	θ	θ
γ₁	~0.8	~0.6	~0.4
γ₂	~0.6	~0.4	~0.2
γ₃	~0.4	~0.2	~0.1

(b) Penetration depths at three different times (schematic presentation)
Figure 9.11 Analysis of continuous medium

defined as

\gamma = 4 \sqrt {a t} \tag {9.124}

where a is the thermal diffusity, ^{2} a = k/ρc. At this distance,

\frac {\theta (\gamma) - \theta_ {i}}{\theta_ {0} - \theta_ {i}} <   0. 0 1 \tag {9.125}

(see Fig.9.11).

This penetration depth is also used for the problem when instead of an imposed temperature, suddenly a heat flux is imposed. In this case we have

\frac {\theta (\gamma) - \theta_ {i}}{\theta^ {s} - \theta_ {i}} <   0. 0 1 \tag {9.126}

where \theta^s is the surface temperature at time t .

The finite element discretization is chosen by use of \gamma . Assume that t_{\mathrm{min}} is the minimum time at which temperature results are desired. Then, if N is the number of