The second law of thermodynamics is used as a variational statement to derive a numerical procedure to satisfy the governing equations of motion. The procedure, based on numerical experimentation, appears to be stable provided the CFL condition is satisfied. This stability is manifested no matter how severe the gradients (compression or expansion) are in the flow field. For reasons of simplicity only one dimensional inviscid compressible unsteady flow is discussed here; however, the concepts and techniques are not restricted to one dimension nor are they restricted to inviscid non-reacting flow. The solution here is explicit in time. Further study is required to determine the impact of the variational principle on implicit algorithms

Prozan, R. J.

English

NASA Technical Reports Server

NASA Contractor Report 3526 
A Variational Principle for 
Compressible Fluid Mechanics 
Discussion of the One-Dimensional Theory 
Robert Joel Prozan 
CONTRACT NAS l- 16646 
APRIL 1982 
NASA 
CR 
3526 
c. 1 
https://ntrs.nasa.gov/search.jsp?R=19820014344 2020-03-21T08:16:41+00:00Z
TECH LIBRARY KAFB, NM 
I##III%#I#IwIIIIIl~ll~~~l 
OCIb22l,4 - 
NASA Contractor Report 3526 
A Variational Principle for 
Compressible Fluid Mechanics 
Discussion of the One-Dimensional Theory 
Robert Joel Prozan 
Continuum, Inc. 
Huntsville, Alabama 
Prepared for 
Langley Research Center 
under Contract NASl-16646 
National Aeronautics 
and Space Administration 
Scientific and Technial 
Information Branch 
1982 
1.0 SUMMARY AND INTRODUCTION 
Conventional numerical analyses of fluid dynamic systems 
variously use finite element or finite difference analogs 
to satisfy the governing equations of motion. To date, 
numerical solutions are more of an art than a science. 
In this report it is proposed that the second law of. 
thermodynamics, a neglected principle of physics, is of 
paramount importance to a successful solution to the 
governing equations of motion. The second law is not new, 
nor is the post priori use of the second law in ascertain- 
ing the physical correctness of the solution. What is 
new is the incorporation of the second law as a fundamental 
statement of motion which is of equal importance to the 
solution as are the familiar conservation equations. It 
is hypothesized that the second law of thermodynamics is 
actually a variational principle and, in point of fact, is 
used as such in the subsequent discussion. 
Many years of experience with the frustrations of numerical 
modeling have convinced the author that something is miss- 
ing in the conventional numerical approaches. Even if the 
conservation equations are satisfied the solution may fail. 
There must be other governing criterion that have hereto- 
fore been lacking. 
The ensuing discussion will use the second law of thermo- 
dynamics as a variational statement to derive a numerical 
procedure which provides insight to some fundamental 
questions not previously resolved. The procedure, based 
on numerical experimentation, appears to be stable pro- 
vided the CFL condition is satisfied. This stability is 
manifested no matter how severe the gradients (compression 
or expansion) are in the flow field. 
For reasons of simplicity only one-dimensional inviscid 
compressible unsteady flow will be discussed here: however, 
the concepts and techniques are not restricted to one 
dimension nor are they restricted to inviscid non-reacting 
flow. The solution here is explicit in time. Further 
study is required to determine the impact of the varia- 
tional principle on implicit algorithms. 
2.0 DISCUSSION 
2.1 The Variational Principle 
In the subsequent discussion vector notation will be used 
rather than the more frequently used Einstein notation. 
The equations of motion for compressible, inviscid, 
unsteady flow of a perfect gas are; 
- g+==o G (1) 
where T = (p,pu,pE); 4 = (Pu,P+Pu2 r PUH) 
and where p = pZT; E = ZT +Lu2; H = E+E 
(y-1) 2 P 
so that p = (y-l) (pE-ipu2) 
In text definitions will be used throughout. The symbols 
introduced above are: p(density); u(axia1 velocity); 
T(temperature); p(pressure); t(time); x(space coordinate); 
Z(gas constant): y(ratio of specific heats): E(specific 
internal + kinetic energy): H(tota1 enthalpy). 
The Second Law of Thermodynamics is written; 
aps + apus -- - 
at ax Lo 
The integral of the above equation over the region of 
interest and over an interval of time is; 
,t"" j, ($ + p)dvdt 2 0 
1 
Equation (3) may be further integrated to yield; 
t2 
$1 I, +v + ~,,p-'dA) dt L 0 
(2) 
(3) 
(4) 
Thus the Second Law states that the entropy generated 
internal to the region under investigation must be station- 
ary or increase. Another view of the Second Law results 
when equation (3) is expanded; 
I, (i@?$+$+ +~(k$+;~+pu2) +@sp) ) dv 
> 0 (5) - 
Equation (5) states that if the conservation equations are 
satisfied exactly everywhere then the system entropy 
generation is zero. In a numerical analysis, however, the 
equations of motion are only satisfied approximately so 
that in general the equality will not be satisfied. The 
Second Law is therefore interpreted as stating that the 
"numericalm entropy created by the inexact scheme must 
be positive. There are many choices made during the con- 
struction of a numerical analog. Whether by accident or 
design, these choices must not (in the long run at least) 
violate the Second Law. 
The equations of motion which must be satisfied can be 
viewed as constraints to the entropy formation. Therefore, 
the constrained entropy function (S) can be written as; 
t2 
- - 
S = jtl( ,,($ + X++$)dv + /,,pusdA)dt (6) 
where the x=(Xp,Xpu,XpE) are the Lagrange multipliers. 
Note that the constraint term looks very much like the 
familiar Method of Weighted Residuals (MWR). Now let the 
time interval be small and let At=t,-t,. The one-dimen- 
sionai domain will consist of N nodes which, for simplicity, 
will be subdivided into (N-l) elements of length Ax. 
(See figure 1 below.) 
e=l e=2 _ 
n=l n=2 n=3 
e=N-_l x 
n=N-1 n=N 
Figure 1 Schematic of region 
3 
The numbering system is such that the element 
bounded by nodes (n) and (n+l) or that element 
bounded by nodes (n) and (n+l). 
We will first develop the constraint term (@). 
and x-xn=Axc. The constraint term is then; 
(e) is 
(n) is 
Let t=AtT 
@ = AtAx/ ==' f /'=l x. (i& 2 + k ?ii)d<& (7) 
T=O e=l c=O 
For the element shown in figure 2 below 
<=O <=l 
Ax 0 n n+Y 
Figure 2 Schematic of element 
the functions 7 and 4 are assumed to vary linearly. That is 
F = (1-E)~n+5~n+,; 4 = (1-s)s,+5qn+1. Then 
(8) 
We will also assume a linear functional for 5; viz; 
x = xl;xn + T2n.xn+l (9) 
In the above expression the superscript (=) is used to 
denote a dyad. Only diagonal dyads will occur during 
this development so that a vector pre or post multiplying 
the dyad will give the same result. There are also some 
occasions when it is necessary to redefine a vector as a 
dyad. It is felt that no undue interpretation problems 
to the reader should result. 
= 
Now Aln =x1,(<) and x2n=x2n(<). These functions are assoc- 
iated with the element (e) (left node (n) and right node (n+l)). 
The constraint term becomes; 
4 
N-1 5=1 = _ = _ . . 
@=AxAt 1 1 
- - v 
n=l E=O 
(A,; ~I,+A~~- xn+l) l ( (1-E) fn+Sfn+l+ “;;-% de 
where (the time integration has been performed), 
en = E+l-z; +n+l = E$$-T$+i 
At At 
and where the superscripts denote the time level. 
Since the stated objective was to create an explicit 
scheme (this has nothing to do with the variational state- 
ment but is of great practical importance), the weight 
functions (to use MWR terminology) are chosen such that: 
! 
E.=l= 
AlnWS=O; f 
5=lT 
C=O c=o 
2n(l-S)dS = 0. 
= = 
If 'lnf '2n are linear we find that; 
= 
A 
=I 
In = bln(+ A2n = g&-$) (J-0) 
=I =I 
where bin, b2n are arbitrary constants. Now ; 
I 
5=1 = = = 
<=o 
(A,,U-S) :A2nS;Alnt x2JdS = (-bln,~2n,-~ln,~2n) 
=I = = 1 = 
where bin = 6bln; b2n = 6b2n. 
Since the g" were arbitrary constants so are the b. The 
constraint term is now written as; 
where, in order go achieve a mgre recpgnizable form a zzz= 
substitution of b2n = (l-(3,); bin = -8, was made. 
5 
The variational statement of equation (6) becomes 
1 $dvdr+@+l 
-c=l 
p&*dxAtdT (12) 
-r=o v -c=o 
Differentiation of equation (12) with respect to the nodal 
Lagrange multipliers yields the equations of constraint; 
as 
Bx,= 8,. $+ =q = 0 Ax 
aF 
- = (f-FNml) l (+N+ 
ghT-gN- ‘- 
a xN 
Ax )=O 
as= . 
axn 
(T-Fn-l+zn) l Tn +Bn’ ( 4n+ii4n) +(T-i5;,-1) l ( S-4 A:-‘) = 0 (13c) 
for 2<n<N-1 -- 
The above equations are an explicit numerical analog to the 
governing,(or constraint) equations of motion. =To deter- 
mine the 7 however, the grb&trary parameters B must be 
found. Le?: for example, fin=0 for all n. Then the interior 
nodal analog becomes; 
while if F,=T we get; 
(f4a) 
(14b) 
which are backward qnd forward differencing schemes 
respectively. The Bn therefore control the typeof 
analog generated. 
control the 
It is also apparent that the 8, 
portion of each element that is associated 
with the node. Figure 3 illustrates the point. 
Figure 3 Effective control volume 
The control volume is seen to be Av =Ax(l-f3 
Returning to the constrained equatik (12) % a 
+B ). 
p&forming 
the integration on the same control volume basis yields: 
-- 
S=AxAt ( 
N-l - 
+ c aps, n=2 af 
’ (T-Fn-l+Bn) l +n) + @ + (~USN-~US~) At (15) 
where the flux term was expanded first order in time and 
integrated with respect to time. Differentiation of the 
above functional with respect to the fn at each node yields: 
ail 
AxA&%.; --!I- apusl +F 
1 2Ax af = 0 
as= 
ai, 
AxAt(e*(%;N-l)+z af 
as= - == 2<n<N-1 
ai, 
AxAt(l-&l+E,,*(e + xn) = 0 __ 
Equations (16) therefore define the Lagrange multipliers 
so that, for an interior node: 
(16a) 
(16b) 
(16~) 
Tin = apsn af (17) 
Finally, using the above results, the constrained functional 
is differentiated with respect to the Bn. The result is; 
(18) 
Equation (18) indicates that for this low order analysis 
it is not possible to find a stationary value of the 
functional. This is of course due to the fact that the 
functional is linear in Bn. What can be done is to determine 
the value of Bn which is in the direction of increasing 
system entropy. If Pn is viewed as an interpolant then it 
must have a value between zero and unity. The approach 
taken for the purposes of the following computations is to 
choose Bn equal to unity if the term in equations (18) is 
positive and equal to zero if the term is negative. In this 
fashion the differences become one-sided in the direction 
of maximum increase (or minimum decrease) of system entropy. 
It is also interesting to note that the boundary conditions 
naturally evolve from the analysis. Suppose that the left- 
most node remains at a fixed condition for all variables. 
In this case B- is not a variable and may not be determined 
by equation (lb), rather it must be zero so that (13a) is 
removed from the equation list. In the diaphragm burst 
problem to be discussed later, the axial velocity must be 
fixed (zero) while density and energy are variable. 
Therefore the value of 6, for the momentum equation must be 
zero while the other two are determined by equation (18). 
It can be seen, therefore, that the analysis determines a 
differencing scheme for each equation at each node for each 
time step which simultaneously satisfies boundary conditions 
and the condition of maximum entropy increase. As will be 
seen, depending on the values of (3n any particular equation 
at any particular node may be differenced forward or back- 
ward or centered or unchanged. 
. 
Equations (13) may be solved for the F at each node. This 
results in a nonconservative analog. If, however, 
equations (13) are solved in the following fashion; 
= 
Ax(~-~,$+~ = -(l-&,-l) l ($&-1) (19b) 
8 
then the result is absolutely conservative. To illustrate 
this, note that; 
I + dV = &AVn = - &-cJl) 
V n=l 
(20) 
so that all the interior flux terms cancel and just the 
entering and exiting flux terms remain. 
To summarize the computational procedure, having known 
values of Fn at a given time level, then 
(1) compute sn=cn<Fn) for all nodes. 
(2) compute aps/af for all equations at all nodes. 
(3) use equations 118) (and boundary conditions) to 
determine the Bn for each element. 
. 
(4) solve equations (19) for (FnAVn) at all nodes. 
(5) integrate for next time level: 
Fk+l 
n = $+&AV,)At/AV 
(6) repeat (l)-(5). 
2.2 Results 
The one-dimensional equations restrict the problems which 
can be solved to demonstrate the method. Two relatively 
difficult problems may be investigated, however. These 
are the shock tube and diaphragm burst. In the former 
case the boundary conditions dictate that B =(O,O,O) while 
B 
t e'(*) E 
are free. In the latter case Bl=EN-l= t *,O,*) where 
indicates a free parameter. In all cases y=1.4 
and Z=l. 
Figure 4 illustrates the behavior of the method for the 
shock tube conditions. The initial disturbance of the 
square wave has not yet washed out. This impulsive start 
typically gives rise to wavelets which migrate back and 
forth between the inlet and the traveling primary wave. 
If the solution domain were long enough a better descrip- 
tion of the wave characteristics would result. In a higher 
spatial dimension analysis, however, economics prohibit 
the use of considerably more than the 40 points in a single 
direction. The traveling wave description shown in figure 4 
9 
is therefore typical of what would occur in a two-dimen- 
sional or three-dimensional analysis. The final velocity 
and pressure distributions are, of course, correct. 
Figure 5 illustrates typical results on a diaphragm burst 
problem. The results shown occur just before the shock 
wave impinges on the downstream wall and just before the 
expansion wave reaches the upstream wall. Calculations 
have been run for many thousands of time steps in which 
the waves reflect many times from both ends. These results 
are not particularly interesting and in the interest of 
brevity are not shown here. 
A typical time step is shown in figure 6. The rightmost 
three columns contain the values of beta for each equation 
at each node for this time step. These values of beta 
correspond to differencing schemes which are identified in 
the column headed by DIFF. The three letters at each node 
identify the scheme constructed for the three equations of 
motion. U means unchanged, B means backward, F means for- 
ward and C means centered. Time steps exceeding the CFL 
condition, as well as negative pressures and/or densities, 
are physically unrealistic. In the context of the varia- 
tional approach they must be viewed as inequality constraints. 
Due to the difficulty of formally entering an inequality 
constraint into the functional definition, these constraints 
were handled in the following fashion: The time step was 
never allowed to exceed some predetermined fraction of the 
least local CFL and should a violation of the pressure or 
density inequality constraint occur then the time step was 
halved. Subsequent time steps were increased by 10 percent 
each step until the maximum time step is reached. It seems 
reasonable to expect that, particularly during initial 
transients, time steps may be attempted which would lead 
to negative pressures or densities. This should only be a 
temporary condition. Instability would be manifested by 
repeated attempts to achieve unrealistic values; fortun- 
ately this never occurred. 
10 
SHOCK TUBE ANALYSIS 
Initial Conditions: p=l,p=l,u=O (field); p=5,p=2.5,~=2 (inlet) 
Boundary Conditions: fixed inlet(left end); free outlet 
/ 
final velocity distribution- --------------. ------- - 
--*4-CC .-C.~-.-.~-..C.-rC*-C~=-.-.v-.~.-~~ 
velocity .- /c 
2 
?i 
1 -2 
2 
if 
0 
1 5 9 13 17 21 25 29 33 37 
Axial node # 
Figure 4 Typical shock tube solution 
40 
DIAPHRAGM BURST ANALYSIS 
Initial Conditions: p=l,p=l,u=O; x>20 p=5,p=.2.5,u=o; xc20 
Boundary Conditions: ends closed 
i 
velocity b-- v-- 
.*/hd 
/ .d 
initial pressure 
Idist.ribution -.. 
1 5 9 13 17 21 25 29 33 37 40 
Axial node 11 
Figure 5 Typical diaphragm burst solution 
Fwo 
. . . . &Cm z, 
2.57916 
2.70773 
2.72572 
2.74249 
2.80334 
2.31434 
2 .32569 
2.32779 
2.82644 
2.3097 1 
2.32753 
2. 318$-J3 
2.84409 
2.33937 
2.36746 
3.13534 
3.2&&S& 
3.10797 
2.72104 
2.83974 
2.73827 
1 .I?1331 
1 .39(-J&5 
1 . 1 GC)34 
1.01275 
1.00’4 1 4 
1 , 00033 
1 . CKm 15 
1 . 0000 1 
u 
2 
1. ‘34482 
1. S&l&l 
1. %4509 
1 .:3378 1 
1.79843 
1.79491 
1 . p-309 
1.. 76735 
1 . 3Qtj*33 
1.81449 
1 .31106 
I.81693 
1 . 8(:)77:3 
1.791i37 
1.78805 
1.63467 
1. &4’@zj 4 
1.7 5955 
1.9472 
1 ‘Fi349 . . c* 
1 n 70407 
1.4TJC~l 4 
.830265 
.2647 13 
7.&2318E-Z 
4.62629E-3 
2.553&&E-3 
1.33Q 1 GE-4 
6. 17lOlE-5 
P 
5 
5.22644 
5. &<I833 
5. b7329 
5.73433 
5.35845 
5.83434 
5.90675 
5.89313 
5.88075 
5.30511 
5.3320 1 
5.7813 
5.94022 
5.87319 
5 99&.37 
&: 4(:)5&2 
B. 4%209 
d. 14304 
5.54245 
5.30401 
5.41724 
2.53593 
2.07443 
1.2473 
1.08433 
1 . cm74 
1. 001 14 
1 .00047 
1. QQQQ3 
.32643 1 
.3272%7 
.327797 
.331944 
.334054 
.324955 
.324039 
.321343 
.3 13495 
.31&342 
.312425 
.303 1 Q5 
.3034;7 
.318407 
. 3 10232 
. 3lke359 
.257113 
.242198 
.227744 
.311042 
.207265 
. 27933 
3 33 153E-2 A. c- 
.268744 
1.27515E-2 
4.322&E-2 
2.793 14E-3 
B.77r5.99E-4 
2.40028E-4 
1. :3824E-5 
Dll-F ---, E BFTA 
__ 2 3 
UUU 
BEB 
CBC 
ucu 
BUB 
BCB 
BFC 
BUF 
CBF 
FBF 
UBU 
Em 
BCH 
CFB 
UFE 
BUB 
Em 
CBC 
UEU 
BPB 
CEC 
FCU 
FUC 
FBIJ 
UBB 
cm 
FPB 
FPB 
FBB 
FBB 
1 0 0 0 
2 Q 0 0 
3 1 0 1 
4 0 1 0 
5 0 Q Q 
G 0 1 0 
7 0 1 1 
3 0 Q 1 
9 1 Q 1 
10 1 0 1 
11 0 Cl 0 
12 Cl 0 0 
13 0 1 0 
14 1 1 0 
15 0 1 0 
14 0 0 0 
17 0 0 0 
13 1 0 1 
19 0 0 0 
20 Q 0 Q 
21 1 0 1 
22 1 1 0 
-2 A-.. 1 0 1 
24 1 0 (:I 
25 0 0 0 
33 1 0 0 
-27 1 0 0 
33 1 0 Q 
Yq &. 1 0 r:) 
3Q 1 Q 0 
Figure 6 Illustration of differencing 
3.0 CONCLUSIONS 
In the preceding discussion it was proposed that the 
Second Law of Thermodynamics is a variational principle 
applicable to compressible gas dynamics. It was shown 
that the principle states that the system entropy must 
not decrease and that the equations of motion act as 
constraints to the entropy formation. It was further 
shown that, at least in an explicit formulation, arbi- 
trary parameters exist which control or dictate the 
type of differencing scheme. The analysis then discusses 
how to select these arbitrary parameters such that the 
variational principle is satisfied while maintaining 
exact conservation. 
The concept is demonstrated numerically for the shock 
tube and diaphragm burst situations. The success of these 
calculations demonstrates the validity of the original 
proposition, i.e., that the Second Law is indeed an 
applicable variational principle and furthermore demon- 
strates the validity of the calculational procedure 
derived from the stated principle. 
A fundamental advance in the state of the art of comp- 
utational fluid mechanics can be expected as a result of 
these findings. 
14 
1. Report No. 2. Government Accession No. 
.-.I 
3. Recipient’s Catalcq No. 
NASA CR-3526 
4. Title and Subtitle 5. Report Date 
A VARIATIONAL PRINCIPLE FOR COMPRESSIBLE FLUID MECHANICS - April 1982 
DISCUSSION OF THE ONE-DIMENSIONAL THEORY 6. Performing Organization Code 
7. Author(s) 6. Performing Organization Report No. 
Robert Joel Prozan CI-TR-0042 
10. Work Unit No. 
9. Performing Organization Name and Address 
Continuum, Inc. 
4717 University Dr. (Suite 96) 11. Contract or Grant No, 
Huntsville, Alabama 35805 NASl-16646 
13. Type of Report and Period Covered 
12. Sponsoring Agency Name and Address Contractor report 
National Aeronautics and Space Administration 
Washington, DC 20546 14. Sponsoring Agency Code 
15. Supplementary Notes 
Langley Technical Monitor: James L. Hunt 
Interim Report 
1 
-1 
1 
6. Abstract 
A variational principle governing compressible fluid mechanics 
is introduced. A numerical analog based on this principle was 
devised and demonstrated for one-dimensional compressible flow 
in shock tube and diaphragm burst problems. The numerical 
behavior of the solution is excellent. The success of the 
numerical solution demonstrates the validity of the variational 
principle. 
r. Key Words (Suggested by Author(s) ) 
Variational principle 
Numerical analysis 
Fluid mechanics 
18. Distribution Statement 
Unclassified - Unlimited 
Subject Category 02 
3. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No. of Pages 22. Price 
Unclassified Unclassified 14 A02 
For sale by the National Technical Information Service, Springfield, Virginia 22161 
NASA-Langley, 1982 


A variational principle for compressible fluid mechanics. Discussion of the one-dimensional theory

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server