-
Notifications
You must be signed in to change notification settings - Fork 59
/
Copy pathfd1d_heat_steady.html
468 lines (413 loc) · 13.9 KB
/
fd1d_heat_steady.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
<html>
<head>
<title>
FD1D_HEAT_STEADY - Finite Difference Solution of a 1D Steady State Heat Equation
</title>
</head>
<body bgcolor="#EEEEEE" link="#CC0000" alink="#FF3300" vlink="#000055">
<h1 align = "center">
FD1D_HEAT_STEADY <br> Finite Difference Solution of a 1D Steady State Heat Equation
</h1>
<hr>
<p>
<b>FD1D_HEAT_STEADY</b>
is a FORTRAN90 program which
applies the finite difference method to estimate the solution of
the steady state heat equation over a one dimensional region, which
can be thought of as a thin metal rod.
</p>
<p>
We will assume the rod extends over the range A <= X <= B.
</p>
<p>
The quantity of interest is the temperature U(X) at each point in the rod.
</p>
<p>
We will assume that the temperature of the rod is fixed at known values
at the two endpoints. Symbolically, we write the <i>boundary conditions</i>:
<pre>
U(A) = UA
U(B) = UB
</pre>
</p>
<p>
Inside the rod, we assume that a version of the steady (time independent)
heat equation applies. This assumes that the situation in the rod has
"settled down", so that the temperature configuration has no further
tendency to change. The equation we will consider is
<pre>
- d/dx ( K(X) * d/dx U(X) ) = F(X)
</pre>
Here, the right hand side term F(X) allows us to consider internal
heat sources in the metal - perhaps a portion of the rod is sitting
above a blow torch. The coefficient K(X) is a measure of heat conductivity.
It measures the rate at which the heat from a local hot spot will spread out.
</p>
<p>
If the heat source function F(X) is zero everywhere, and if K(X) is
constant, then the solution U(X) will be the straight line function
that matches the two endpoint values. Making F(X) positive over a small
interval will "heat up" that portion. You can simulate a rod that is
divided into regions of different materials by setting the function K(X)
to have a given value K1 over some subinteral of [A,B] and value K2 over
the rest of the region.
</p>
<p>
To estimate the value of the solution, we will pick a uniform mesh
of N points X(1) through X(N), from A to B. At the I-th point, we will
compute an estimated temperature U(I). To do this, we will need to
use the boundary conditions and the differential equation.
</p>
<p>
Since X(1) = A and X(N) = B, we can use the boundary conditions to
set U(1) = UA and U(N) = UB.
</p>
<p>
For the points in the interior, we need to approximate the differential
equation in a way that allows us to determine the solution values.
We will do this using a finite difference approximation.
</p>
<p>
Suppose we are working at node X(I), which is associated with U(I).
Then a centered difference approximation to
<pre>
- d/dx ( K(X) * d/dx U(X) )
</pre>
is
<pre>
- ( + K(X(I)+DX/2) * ( U(X(I+1)) - U(X(I)) ) / DX )
- K(X(I)-DX/2) * ( U(X(I)) - U(X(I-1)) ) / DX ) / DX
</pre>
</p>
<p>
If we rearrange the terms in this approximation, and set it equal to F(X(I)),
we get the finite difference approximation to the differential equation at X(I):
<pre>
- K(X(I)-DX/2) * U(X(I-1)
+ ( K(X(I)-DX/2) + K(X(I)+DX(2)) ) * U(X(I))
- K(X(I)+DX(2)) * U(X(I+1)) = DX * DX * F(X(I))
</pre>
</p>
<p>
This means that we have N-2 equations, each of which involves
an unknown U(I) and its left and right neighbors, plus the two boundary
conditions, which give us a total of N equations in N unknowns.
</p>
<p>
We can set up and solve these linear equations using a matrix A for
the coefficients, and a vector RHS for the right hand side terms,
and calling a function to solve the system A*U=RHS will give us the solutioni.
</p>
<p>
Because finite differences are only an approximation to derivatives, this
process only produces estimates of the solution. Usually, we can reduce
this error by decreasing DX.
</p>
<p>
This program assumes that the user will provide a calling program, and
functions to evaluate K(X) and F(X).
</p>
<h3 align = "center">
Usage:
</h3>
<p>
<pre>
call <b>fd1d_steady_heat</b> ( <i>n, a, b, ua, ub, k, f, x, u</i> )
</pre>
where
<ul>
<li>
<i>n</i> is the number of spatial points.
</li>
<li>
<i>a, b</i> are the left and right endpoints.
</li>
<li>
<i>ua, ub</i> are the temperature values at the left and right endpoints.
</li>
<li>
<i>k</i> is the name of the function which evaluates K(X);
</li>
<li>
<i>f</i> is the name of the function which evaluates the right hand side F(X).
</li>
<li>
<i>x</i> is the X coordinates of nodes (output);
</li>
<li>
<i>u</i> is the computed value of the temperature at the nodes.
</li>
</ul>
</p>
<h3 align = "center">
Licensing:
</h3>
<p>
The computer code and data files described and made available on this web page
are distributed under
<a href = "../../txt/gnu_lgpl.txt">the GNU LGPL license.</a>
</p>
<h3 align = "center">
Languages:
</h3>
<p>
<b>FD1D_HEAT_STEADY</b> is available in
<a href = "../../c_src/fd1d_heat_steady/fd1d_heat_steady.html">a C version</a> and
<a href = "../../cpp_src/fd1d_heat_steady/fd1d_heat_steady.html">a C++ version</a> and
<a href = "../../f77_src/fd1d_heat_steady/fd1d_heat_steady.html">a FORTRAN77 version</a> and
<a href = "../../f_src/fd1d_heat_steady/fd1d_heat_steady.html">a FORTRAN90 version</a> and
<a href = "../../m_src/fd1d_heat_steady/fd1d_heat_steady.html">a MATLAB version</a>
</p>
<h3 align = "center">
Related Data and Programs:
</h3>
<p>
<a href = "../../data/fd1d/fd1d.html">
FD1D</a>,
a data directory which
contains examples of 1D FD files, two text files that can be used to describe
many finite difference models with one space variable, and either no time dependence
or a snapshot at a given time;
</p>
<p>
<a href = "../../f_src/fd1d_burgers_lax/fd1d_burgers_lax.html">
FD1D_BURGERS_LAX</a>,
a FORTRAN90 program which
applies the finite difference method and the Lax-Wendroff method
to solve the non-viscous time-dependent Burgers equation
in one spatial dimension.
</p>
<p>
<a href = "../../f_src/fd1d_burgers_leap/fd1d_burgers_leap.html">
FD1D_BURGERS_LEAP</a>,
a FORTRAN90 program which
applies the finite difference method and the leapfrog approach
to solve the non-viscous time-dependent Burgers equation in one spatial dimension.
</p>
<p>
<a href = "../../f_src/fd1d_bvp/fd1d_bvp.html">
FD1D_BVP</a>,
a FORTRAN90 program which
applies the finite difference method
to a two point boundary value problem in one spatial dimension.
</p>
<p>
<a href = "../../m_src/fd1d_display/fd1d_display.html">
FD1D_DISPLAY</a>,
a MATLAB program which
reads a pair of files defining a 1D finite difference model, and plots the data.
</p>
<p>
<a href = "../../f_src/fd1d_heat_explicit/fd1d_heat_explicit.html">
FD1D_HEAT_EXPLICIT</a>,
a FORTRAN90 program which
uses the finite difference method and explicit time stepping
to solve the time dependent heat equation in 1D.
</p>
<p>
<a href = "../../f_src/fd1d_heat_implicit/fd1d_heat_implicit.html">
FD1D_HEAT_IMPLICIT</a>,
a FORTRAN90 program which
uses the finite difference method and implicit time stepping
to solve the time dependent heat equation in 1D.
</p>
<p>
<a href = "../../f_src/fd1d_predator_prey/fd1d_predator_prey.html">
FD1D_PREDATOR_PREY</a>,
a FORTRAN90 program which
implements a finite difference algorithm for a predator-prey system
with spatial variation in 1D.
</p>
<p>
<a href = "../../f_src/fd1d_wave/fd1d_wave.html">
FD1D_WAVE</a>,
a FORTRAN90 program which
applies the finite difference method to solve the time-dependent
wave equation utt = c * uxx in one spatial dimension.
</p>
<p>
<a href = "../../m_src/fem_50_heat/fem_50_heat.html">
FEM_50_HEAT</a>,
a MATLAB program which
implements a finite element calculation for the heat equation.
</p>
<p>
<a href = "../../f_src/fem1d_heat_steady/fem1d_heat_steady.html">
FEM1D_HEAT_STEADY</a>,
a FORTRAN90 program which
uses the finite element method to solve the steady (time independent)
heat equation in 1D.
</p>
<p>
<a href = "../../f_src/fem2d_heat/fem2d_heat.html">
FEM2D_HEAT</a>,
a FORTRAN90 program which
solves the 2D time dependent heat equation on the unit square.
</p>
<h3 align = "center">
Reference:
</h3>
<p>
<ol>
<li>
George Lindfield, John Penny,<br>
Numerical Methods Using MATLAB,<br>
Second Edition,<br>
Prentice Hall, 1999,<br>
ISBN: 0-13-012641-1,<br>
LC: QA297.P45.
</li>
</ol>
</p>
<h3 align = "center">
Source Code:
</h3>
<p>
<ul>
<li>
<a href = "fd1d_heat_steady.f90">fd1d_heat_steady.f90</a>,
the library.
</li>
<li>
<a href = "fd1d_heat_steady.sh">fd1d_heat_steady.sh</a>,
commands to compile the library.
</li>
</ul>
</p>
<h3 align = "center">
Examples and Tests:
</h3>
<p>
<ul>
<li>
<a href = "problem1.f90">problem1.f90</a>,
uses K(X) = 1, F(X) = 0, so the solution should be the straight
line that connects the boundary values.
</li>
<li>
<a href = "problem1.sh">problem1.sh</a>,
commands to compile the problem and run it with the library.
</li>
<li>
<a href = "problem1_nodes.txt">problem1_nodes.txt</a>,
the coordinates of the nodes.
</li>
<li>
<a href = "problem1_values.txt">problem1_values.txt</a>,
the computed temperatures at the nodes.
</li>
<li>
<a href = "problem1.png">problem1.png</a>,
a PNG image of the solution.
</li>
</ul>
</p>
<p>
<ul>
<li>
<a href = "problem2.f90">problem2.f90</a>,
uses K(X) which is set to different constants over three subregions,
and F(X) = 0.0, so the solution will be a piecewise linear function
that connects the boundary values.
</li>
<li>
<a href = "problem2.sh">problem2.sh</a>,
commands to compile the problem and run it with the library.
</li>
<li>
<a href = "problem2_nodes.txt">problem2_nodes.txt</a>,
the coordinates of the nodes.
</li>
<li>
<a href = "problem2_values.txt">problem2_values.txt</a>,
the computed temperatures at the nodes.
</li>
<li>
<a href = "problem2.png">problem2.png</a>,
a PNG image of the solution.
</li>
</ul>
</p>
<p>
<ul>
<li>
<a href = "problem3.f90">problem3.f90</a>,
uses K(X) = 1, F(X) defines a heat source, so the solution can
rise above the boundary values.
</li>
<li>
<a href = "problem3.sh">problem3.sh</a>,
commands to compile the problem and run it with the library.
</li>
<li>
<a href = "problem3_nodes.txt">problem3_nodes.txt</a>,
the coordinates of the nodes.
</li>
<li>
<a href = "problem3_values.txt">problem3_values.txt</a>,
the computed temperatures at the nodes.
</li>
<li>
<a href = "problem3.png">problem3.png</a>,
a PNG image of the solution.
</li>
</ul>
</p>
<p>
<ul>
<li>
<a href = "problem4.f90">problem4.f90</a>,
uses K(X) = 1, F(X) defines a heat source and a heat sink, so the
solution can go above and below the boundary values.
</li>
<li>
<a href = "problem4.sh">problem4.sh</a>,
commands to compile the problem and run it with the library.
</li>
<li>
<a href = "problem4_nodes.txt">problem4_nodes.txt</a>,
the coordinates of the nodes.
</li>
<li>
<a href = "problem4_values.txt">problem4_values.txt</a>,
the computed temperatures at the nodes.
</li>
<li>
<a href = "problem4.png">problem4.png</a>,
a PNG image of the solution.
</li>
</ul>
</p>
<h3 align = "center">
List of Routines:
</h3>
<p>
<ul>
<li>
<b>FD1D_HEAT_STEADY</b> solves the steady 1D heat equation.
</li>
<li>
<b>GET_UNIT</b> returns a free FORTRAN unit number.
</li>
<li>
<b>R83NP_FS</b> factors and solves an R83NP system.
</li>
<li>
<b>R8MAT_WRITE</b> writes an R8MAT file.
</li>
<li>
<b>TIMESTAMP</b> prints the current YMDHMS date as a time stamp.
</li>
</ul>
</p>
<p>
You can go up one level to <a href = "../f_src.html">
the FORTRAN90 source codes</a>.
</p>
<hr>
<i>
Last revised on 13 September 2010.
</i>
<!-- John Burkardt -->
</body>
</html>